#include "PndRadMapBoxMesh.h" #include #include #include #include #include #include #include #include #include "PndRadMapPlane.h" #define verbose true bool momentumfilter(TVector3 mom, TVector3 diff){ bool back = true; for(int i = 0; i < 3; i++) if(mom(i)*diff(i) < 0) back = false; return back; } // PndRadMapPlane::PndRadMapPlane(){ // SetNormal(TVector3(0, 0, 0)); // SetDistance(0); // corner1 = TVector3(0, 0, 0); // corner2 = TVector3(0, 0, 0); // corner3 = TVector3(0, 0, 0); // } // PndRadMapPlane::PndRadMapPlane(TVector3 _corner1, TVector3 _corner2, TVector3 _corner3, // double dist){ // corner1 = _corner1; // corner2 = _corner2; // corner3 = _corner3; // TVector3 dpv1 = corner2-corner1; // TVector3 dpv2 = corner3-corner1; // TVector3 n = (dpv1.Cross(dpv2)).Unit(); // SetNormal(n); // SetDistance(dist); // } // TVector3 PndRadMapPlane::GetCorner(int i){ // switch(i){ // case 1: // return Corner1(); // case 2: // return Corner2(); // case 3: // return Corner3(); // default: // return TVector3(0, 0, 0); // } // } // TVector3 intersectLine(TVector3 l0, TVector3 l1,// two points of the line // TVector3 p0, TVector3 p1, TVector3 p2, // float& eff) {//three points of the plane // //axes of the plane // TVector3 dpv1 = p1-p0; // TVector3 dpv2 = p2-p0; // TVector3 n = (dpv1.Cross(dpv2)).Unit(); // // Rotate the points (of the line) to the plane // TVector3 vRotRay1 = TVector3 (dpv1*(l0-p1), dpv2*(l0-p1), n*(l0-p1)); // TVector3 vRotRay2 = TVector3 (dpv1*(l1-p1), dpv2*(l1-p1), n*(l1-p1)); // float fPercent = vRotRay1.Z() / (vRotRay1.Z()-vRotRay2.Z()); // eff = fPercent; // TVector3 ixp = l0 + fPercent*(l1-l0); // } // TVector3 PndRadMapPlane::LineIntersection(TVector3 begline, TVector3 endline){ // TVector3 diffl = endline-begline; // float nDotBeg = normal*(corner1-begline); // float nDotDiff = normal*diffl; // float tt = nDotBeg/nDotDiff; // TVector3 ixp = begline + tt*diffl; // return ixp; // } PndRadMapBoxMesh::PndRadMapBoxMesh(PndRadMapBoxMesh& m){ _MeshHisto = (TH2D*)m._MeshHisto->Clone(); _StatHisto = (TH2I*)m._StatHisto->Clone(); _orientation = m._orientation; _rotate = m._rotate; _volume = m._volume; _Name = m._Name; _plane = NULL; _isSurfaceQuantity = false; pdg = TDatabasePDG::Instance(); pdg -> ReadPDGTable(); } PndRadMapBoxMesh::PndRadMapBoxMesh(const char* name, int Xbins, Double_t Xlow, Double_t Xhigh, int Ybins, Double_t Ylow, Double_t Yhigh, int Zbins, Double_t Zlow, Double_t Zhigh){//in cm!!!! _verbose = 0; _Name = TString(name); _Xbins = Xbins; _Xlow = Xlow; _Xhigh = Xhigh; _Ybins = Ybins; _Ylow = Ylow; _Yhigh = Yhigh; _Zbins = Zbins; _Zlow = Zlow; _Zhigh = Zhigh; _filter = TFormula("formula", "1"); if((_Xlow == _Xhigh) || (_Ylow == _Yhigh) || (_Zlow == _Zhigh)) _isSurfaceQuantity = true; pdg = TDatabasePDG::Instance(); pdg -> ReadPDGTable(); // if(_verbose) std::cout << "PndRadMapBoxMesh::PndRadMapBoxMesh()\n" << _Xlow << ' ' << _Xhigh << ' ' << _Ylow << ' ' << _Yhigh << ' ' << _Zlow << ' ' << _Zhigh << ' ' << _isSurfaceQuantity << std::endl; } PndRadMapBoxMesh::PndRadMapBoxMesh(const char* name, int xbins, Double_t xlow, Double_t xhigh, int ybins, Double_t ylow, Double_t yhigh, Double_t zlow, Double_t zhigh, orientation plane, quantity Quantity){//in cm!!!! _verbose = 0; _Name = TString(name); _orientation = plane; SetQuantity(Quantity); switch (_orientation){ case XY: _Xbins = xbins; _Xlow = xlow; _Xhigh = xhigh; _Ybins = ybins; _Ylow = ylow; _Yhigh = yhigh; _Zlow = zlow; _Zhigh = zhigh; break; case YX: _Xbins = ybins; _Xlow = ylow; _Xhigh = yhigh; _Ybins = xbins; _Ylow = xlow; _Yhigh = xhigh; _Zlow = zlow; _Zhigh = zhigh; break; case XZ: _Xbins = xbins; _Xlow = xlow; _Xhigh = xhigh; _Ylow = zlow; _Yhigh = zhigh; _Zbins = ybins; _Zlow = ylow; _Zhigh = yhigh; break; case ZX: _Xbins = ybins; _Xlow = ylow; _Xhigh = yhigh; _Zbins = xbins; _Zlow = xlow; _Zhigh = xhigh; _Ylow = zlow; _Yhigh = zhigh; break; case YZ: _Xlow = zlow; _Xhigh = zhigh; _Ybins = xbins; _Ylow = xlow; _Yhigh = xhigh; _Zbins = ybins; _Zlow = ylow; _Zhigh = yhigh; break; case ZY: _Xlow = zlow; _Xhigh = zhigh; _Ybins = ybins; _Ylow = ylow; _Yhigh = yhigh; _Zbins = xbins; _Zlow = xlow; _Zhigh = xhigh; break; }; if((_Xlow == _Xhigh) || (_Ylow == _Yhigh) || (_Zlow == _Zhigh)) _isSurfaceQuantity = true; pdg = TDatabasePDG::Instance(); pdg -> ReadPDGTable(); } PndRadMapBoxMesh::~PndRadMapBoxMesh(){ if(_verbose) std::cout << "PndRadMapBoxMesh::~PndRadMapBoxMesh()\n"; delete _StatHisto; delete _MeshHisto; } void PndRadMapBoxMesh::SetFilter(const char* filter){ //Pid=x, Ch=y, Mom=z, if(_verbose) std::cout << "PndRadMapBoxMesh::SetFilter(" << filter << ")\n"; TString temps(filter); temps.ReplaceAll("Pid", "x"); temps.ReplaceAll("Ch" , "y"); temps.ReplaceAll("Mom", "z"); _filter = TFormula("formula", temps.Data()); } void PndRadMapBoxMesh::SetQuantity(quantity Quantity){ if(_verbose) std::cout << "PndRadMapBoxMesh::SetQuantity(" << Quantity << ")\n"; _quantity = Quantity; } void PndRadMapBoxMesh::SetOrientation(orientation plane, Double_t rotate, axis Ax){ if(_verbose) std::cout << "PndRadMapBoxMesh::SetOrientation(" << plane << ", " << rotate << ", " << Ax << ")\n"; _orientation = plane; SetOrientation(rotate,Ax); } void PndRadMapBoxMesh::SetOrientation(Double_t rotate, axis Ax){ if(_verbose) std::cout << "PndRadMapBoxMesh::SetOrientation(" << rotate << ", " << Ax << ")\n"; _axis = Ax; _rotate = (rotate == 99999) ? 0 : rotate; if(_verbose) std::cout << "Volume[" << _Name.Data() << "][(" << _Xhigh << ", " << _Xlow << ", " << _Xbins << "),(" << _Yhigh << ", " << _Ylow << ", " << _Ybins << "),(" << _Zhigh << ", " << _Zlow << ", " << _Zbins << ")]:\n"; // char C[32]; switch (_orientation){ case XY: makeHisto("XY", _rotate, _Xbins, _Xlow, _Xhigh, _Ybins, _Ylow, _Yhigh, _Zlow, _Zhigh); break; case YX: makeHisto("YX", _rotate, _Ybins, _Ylow, _Yhigh, _Xbins, _Xlow, _Xhigh, _Zlow, _Zhigh); break; case XZ: makeHisto("XZ", _rotate, _Xbins, _Xlow, _Xhigh, _Zbins, _Zlow, _Zhigh, _Ylow, _Yhigh); break; case ZX: makeHisto("ZX", _rotate, _Zbins, _Zlow, _Zhigh, _Xbins, _Xlow, _Xhigh, _Ylow, _Yhigh); break; case YZ: makeHisto("YZ", _rotate, _Ybins, _Ylow, _Yhigh, _Zbins, _Zlow, _Zhigh, _Xlow, _Xhigh); break; case ZY: makeHisto("ZY", _rotate, _Zbins, _Zlow, _Zhigh, _Ybins, _Ylow, _Yhigh, _Xlow, _Xhigh); break; }; if(!_isSurfaceQuantity) _volume = ((_Xhigh-_Xlow)/_Xbins) * ((_Yhigh-_Ylow)/_Ybins)* ((_Zhigh-_Zlow)/_Zbins); else _volume = 1; if(verbose) std::cout << "Volume: " << (_Xhigh-_Xlow)/_Xbins << ", " << (_Yhigh-_Ylow)/_Ybins << ", " << (_Zhigh-_Zlow)/_Zbins << " -> " << _volume << "\n"; } void PndRadMapBoxMesh::SetVerbosityLevel(int Verbose){ _verbose = Verbose; } // void PndRadMapBoxMesh::transform(Double_t Vh, Double_t Vv){ void PndRadMapBoxMesh::Transform(Double_t X, Double_t Y, Double_t Z){ if(_verbose) std::cout << "PndRadMapBoxMesh::Transform(" << X << ", " << Y << ", " << Z << ")\n"; _X = X; _Y = Y; _Z = Z; Transform( X, Y, Z, _X, _Y, _Z); } void PndRadMapBoxMesh::Transform(Double_t X, Double_t Y, Double_t Z, Double_t& X0, Double_t& Y0, Double_t& Z0){ if(_verbose){ std::cout << "PndRadMapBoxMesh::Transform(" << X << ", " << Y << ", " << Z << ")\n"; std::cout << " " << X0 << ", " << Y0 << ", " << Z0 << ")\n"; } if(_rotate == 99999) _rotate = 0; if((_rotate == 0)){ X0 = X; Y0 = Y; Z0 = Z; }else{ Double_t V1 = 0; Double_t V2 = 0; X0 = X; Y0 = Y; Z0 = Z; switch (_axis){ case Xx: V1 = Y; V2 = Z; X0 = X; break; case Yy: V1 = X; V2 = Z; Y0 = Y; break; case Zz: V1 = X; V2 = Y; Z0 = Z; break; }; Double_t _Vh, _Vv; if(_rotate != 0){ Double_t R = TMath::Sqrt(V1*V1 + V2*V2); Double_t rot0 = TMath::ATan(V1/V2)*TMath::RadToDeg(); if(V2 < 0) rot0 += 180; Double_t drot = rot0-_rotate; _Vh = R*TMath::Sin(1*drot*TMath::DegToRad()); _Vv = R*TMath::Cos(1*drot*TMath::DegToRad()); }else{ _Vh = V1; _Vv = V2; } switch (_axis){ case Xx: switch (_orientation){ case XY: case YX: case XZ: case ZX: Y0 = _Vv; Z0 = _Vh; break; case YZ: Y0 = _Vh; Z0 = _Vv; break; case ZY: Y0 = _Vv; Z0 = _Vh; break; }; break; case Yy: switch (_orientation){ case XY: case YX: case YZ: case ZY: X0 = _Vh; Z0 = _Vv; break; case XZ: X0 = _Vh; Z0 = _Vv; break; case ZX: X0 = _Vv; Z0 = _Vh; break; }; break; case Zz: switch (_orientation){ case XZ: case ZX: case YZ: case ZY: X0 = _Vh; Y0 = _Vv; break; case XY: X0 = _Vh; Y0 = _Vv; break; case YX: X0 = _Vv; Y0 = _Vh; break; }; break; }; } if(_verbose){ std::cout << " -> PndRadMapBoxMesh::Transform(" << X << ", " << Y << ", " << Z << ")\n"; std::cout << " " << X0 << ", " << Y0 << ", " << Z0 << ")\n"; } } void PndRadMapBoxMesh::Fill(FairRadMapPoint *p){ if(_verbose) std::cout << "void PndRadMapBoxMesh::Fill(FairRadMapPoint *p)\n"; Int_t pid = p->GetPdg(); Double_t mom = TMath::Sqrt((p->GetPx()*p->GetPx()) + (p->GetPy()*p->GetPy()) + (p->GetPz()*p->GetPz())); pdgpart = pdg->GetParticle(pid); Int_t cha = -1000;// to be implemented Double_t pdgmass = -1; Double_t val = 0; Double_t edep = 0; Double_t dens = 0; Double_t mass = p->GetMass(); TVector3 poststepv, prestepv; Double_t lX = 0, lY = 0, lZ = 0; if(pdgpart){ pdgmass = pdgpart->Mass(); cha = pdgpart->Charge(); } if(_filter.Eval(pid, cha, mom)){ val = 0; edep = p->GetEnergyLoss()*1.60217657*1e-10;//GeV->Joule dens = p->GetDensity(); dens *= 0.001;//g/cm3 -> kg/cm3 mass = p->GetMass(); poststepv = TVector3(p->GetXOut(), p->GetYOut(), p->GetZOut()); prestepv = TVector3(p->GetX(), p->GetY(), p->GetZ()); lX = 0; lY = 0; lZ = 0; if(_verbose) std::cout << "Density: " << dens << ", mass: " << mass << " (" << (dens*_volume) << ", " << _volume << ") " << dens*mass << ' ' << "charge: " << cha << "\n"; lX = p->GetXOut(); lY = p->GetYOut(); lZ = p->GetZOut(); bool OK = false; switch (_quantity){ case Edep: val = edep; if(!isnan(val) && val !=0) // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val); Fill(lX, lY, lZ, val); break; case Dose: val = edep/(dens*_volume);//mass=density*vol,dens[kg/cm^3],vol[cm^3]! if(!isnan(val) && val !=0){ // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val); Fill(lX, lY, lZ, val); } break; case Fluence: InterSection = _plane->LineIntersection(poststepv, prestepv); switch (_orientation){ case XY: case YX: if(((prestepv.Z() <= InterSection.Z()) && (InterSection.Z() <= poststepv.Z())) || ((poststepv.Z() <= InterSection.Z()) && (InterSection.Z() <= prestepv.Z()))) OK = true; break; case XZ: case ZX: if(((prestepv.Y() <= InterSection.Y()) && (InterSection.Y() <= poststepv.Y())) || ((poststepv.Y() <= InterSection.Y()) && (InterSection.Y() <= prestepv.Y()))) OK = true; break; case YZ: case ZY: if(((prestepv.X() <= InterSection.X()) && (InterSection.X() <= poststepv.X())) || ((poststepv.X() <= InterSection.X()) && (InterSection.X() <= prestepv.X()))) OK = true; break; }; // std::cout << "test0\n" // << "[" << prestepv.X() << "<=" << InterSection.X() << "<=" << poststepv.X() << "] - " << OK << "\n" // << "[" << poststepv.X() << "<=" << InterSection.X() << "<=" << prestepv.X() << "] - " << OK << "\n" // << "[" << prestepv.Y() << "<=" << InterSection.Y() << "<=" << poststepv.Y() << "] - " << OK << "\n" // << "[" << poststepv.Y() << "<=" << InterSection.Y() << "<=" << prestepv.Y() << "] - " << OK << "\n" // << "[" << prestepv.Z() << "<=" << InterSection.Z() << "<=" << poststepv.Z() << "] - " << OK << "\n" // << "[" << poststepv.Z() << "<=" << InterSection.Z() << "<=" << prestepv.Z() << "] - " << OK << "\n\n"; // // if(((prestepv.Z() <= InterSection.Z()) && (InterSection.Z() <= poststepv.Z())) || // // ((poststepv.Z() <= InterSection.Z()) && (InterSection.Z() <= prestepv.Z()))){ if(OK){ // std::cout << "testing 0...\n"; // std::cout << "testing 0... " // << _Xlow << "<=" << InterSection.X() << "<=" << _Xhigh << " -> " // << ((_Xlow <= InterSection.X()) && (InterSection.X() <= _Xhigh)) << ' ' << (InterSection.X() <= _Xhigh) << ' ' << (_Xlow <= InterSection.X()) << ' ' // << _Ylow << "<=" << InterSection.Y() << "<=" << _Yhigh << " -> " // << ((_Ylow <= InterSection.Y()) && (InterSection.Y() <= _Yhigh)) << ' ' // << _Zlow << "<=" << InterSection.Z() << "<=" << _Zhigh << " -> " // << ((_Zlow <= InterSection.Z()) && (InterSection.Z() <= _Zhigh)) << std::endl; Fill(InterSection.X(), InterSection.Y(), InterSection.Z()); //normalisation to the area ... now it is OK Double_t theta = InterSection.Theta()*TMath::RadToDeg(); Double_t K = TMath::Sqrt(mom*mom + pdgmass*pdgmass) - pdgmass; _EnergyHisto->Fill(theta, K); //uncommented by MP, Sep 8 // val = CalcFluence(p); // if(!isnan(val) && val !=0){ // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val); } break; case Density: val = dens; if(!isnan(val) && val !=0){ Fill(lX, lY, lZ, val); // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val); } break; case Mass: val = (dens*_volume); if(!isnan(val) && val !=0){ Fill(lX, lY, lZ, val); // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val); } break; case SimpleFluence:{ poststepv = TVector3(p->GetXOut(), p->GetYOut(), p->GetZOut()); prestepv = TVector3(p->GetX(), p->GetY(), p->GetZ()); val = (poststepv-prestepv).Mag()/1e2;// cm->m val /= (_volume/1e6);//cm^3->m^3 -> 1/m^2 if(!isnan(val) && val!=0){ if(_verbose > 1){ std::cout << "Pre: " << prestepv.X() << ", " << prestepv.Y() << ", " << prestepv.Z() << "\n"; std::cout << "Pre: " << poststepv.X() << ", " << poststepv.Y() << ", " << poststepv.Z() << "\n"; std::cout << "diff: " << (poststepv-prestepv).Mag() << "; val: " << val << std::endl; std::cout << "volume: " << _volume << "\n\n"; } Fill(lX, lY, lZ, val); // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val); } }break; case Twos: val = 2; if(!isnan(val) && val !=0){ Fill(lX, lY, lZ, val); // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val); } break; case EnergyFluence: break; case Flux: break; case Kerma: break; } } } // void PndRadMapBoxMesh::Fill(FairRadMapPoint *p, Double_t we){ void PndRadMapBoxMesh::Fill(Int_t gBin, Double_t val){ if(_verbose) std::cout << "PndRadMapBoxMesh::Fill(" << gBin << ", " << val << ")\n"; Int_t binX, binY, binZ; _MeshHisto->GetBinXYZ(gBin, binX, binY, binZ); Fill(binX, binY, binZ, val); } void PndRadMapBoxMesh::Fill(Double_t X, Double_t Y, Double_t Z, Double_t we){ if(_verbose) std::cout << "PndRadMapBoxMesh::Fill(" << X << ", " << Y << ", " << Z << ", " << we << ")\n"; Transform(X, Y, Z); // std::cout << "0 - PndRadMapBoxMesh::Fill(" << X << ", " << Y << ", " << Z << ", " << we << ")\n"; // std::cout << "1 - PndRadMapBoxMesh::Fill(" << _X << ", " << _Y << ", " << _Z << ", " << we << ")\n"; // if(IsInside(_X, _Y, _Z)){ if(IsInside()){ // std::cout << "2 - PndRadMapBoxMesh::Fill(" << _X << ", " << _Y << ", " << _Z << ", " << we << ")\n"; switch (_orientation){ case XY: _MeshHisto->Fill(_X, _Y, we); _StatHisto->Fill(_X, _Y); break; case YX: _MeshHisto->Fill(_Y, _X, we); _StatHisto->Fill(_Y, _X); break; case XZ: _MeshHisto->Fill(_X, _Z, we); _StatHisto->Fill(_X, _Z); break; case ZX: _MeshHisto->Fill(_Z, _X, we); _StatHisto->Fill(_Z, _X); break; case YZ: _MeshHisto->Fill(_Y, _Z, we); _StatHisto->Fill(_Y, _Z); break; case ZY: _MeshHisto->Fill(_Z, _Y, we); _StatHisto->Fill(_Z, _Y); break; }; } } void PndRadMapBoxMesh::Scale(Double_t sca){ if(_verbose) std::cout << "PndRadMapBoxMesh::Scale(" << sca << ")" << std::endl; _MeshHisto->Scale(sca); } TH2D* PndRadMapBoxMesh::GetHisto(){ if(_verbose) std::cout << "PndRadMapBoxMesh::GetHisto()\n"; return _MeshHisto; } PndRadMapPlane* PndRadMapBoxMesh::GetPlane(){ if(_verbose) std::cout << "PndRadMapBoxMesh::GetPlane()\n"; return _plane; } void PndRadMapBoxMesh::Save(TFile* fout){ if(_verbose) std::cout << "PndRadMapBoxMesh::Save(TFile* fout)\n"; TObjArray* toa = TString(_MeshHisto->GetName()).Tokenize('_'); TString dirname = ((TObjString*)toa->At(0))->GetString(); dirname += '_'; dirname += ((TObjString*)toa->At(2))->GetString().Data(); dirname += '_'; dirname += ((TObjString*)toa->At(3))->GetString().Data(); std::cout << dirname.Data() << std::endl; // fout->mkdir(_MeshHisto->GetName()); // fout->cd (_MeshHisto->GetName()); fout->mkdir(dirname.Data()); fout->cd (dirname.Data()); _MeshHisto->Clone()->Write(); _StatHisto->Clone()->Write(); if(_isSurfaceQuantity) _EnergyHisto->Clone()->Write(); fout->cd(); } void PndRadMapBoxMesh::Save(){ if(_verbose) std::cout << "PndRadMapBoxMesh::Save()\n"; _MeshHisto->Write(); _StatHisto->Write(); if(_isSurfaceQuantity) _EnergyHisto->Clone()->Write(); } bool PndRadMapBoxMesh::IsInside(Double_t X, Double_t Y, Double_t Z){ if(_verbose) std::cout << "PndRadMapBoxMesh::IsInside(" << X << ", " << Y << ", " << Z << ")\n"; bool ok = false; if(((_Xlow <= X) && (X <= _Xhigh)) && ((_Ylow <= Y) && (Y <= _Yhigh)) && ((_Zlow <= Z) && (Z <= _Zhigh))) ok = true; return ok; } bool PndRadMapBoxMesh::IsInside(){ if(_verbose) std::cout << "PndRadMapBoxMesh::IsInside()\n"; bool ok = false; switch (_orientation){ case XY: case YX: if(((_Xlow <= _X) && (_X <= _Xhigh)) && ((_Ylow <= _Y) && (_Y <= _Yhigh))) ok = true; break; case XZ: case ZX: if(((_Xlow <= _X) && (_X <= _Xhigh)) && ((_Zlow <= _Z) && (_Z <= _Zhigh))) ok = true; break; case YZ: case ZY: if(((_Ylow <= _Y) && (_Y <= _Yhigh)) && ((_Zlow <= _Z) && (_Z <= _Zhigh))) ok = true; break; }; // if(((_Xlow <= _X) && (_X <= _Xhigh)) && // ((_Ylow <= _Y) && (_Y <= _Yhigh)) && // ((_Zlow <= _Z) && (_Z <= _Zhigh))) // ok = true; return ok; } bool PndRadMapBoxMesh::IsInside(FairRadMapPoint *p){ if(_verbose) std::cout << "PndRadMapBoxMesh::IsInside(FairRadMapPoint *p)\n"; bool ok = false; double xin, yin, zin, xout, yout, zout; Transform(p->GetX(), p->GetY(), p->GetZ(), xin, yin, zin); Transform(p->GetXOut(), p->GetYOut(), p->GetZOut(), xout, yout, zout); if((((xin <= _Xlow) && (_Xhigh <= xout)) || ((xout <= _Xlow) && (_Xhigh <= xin )) || ((_Xlow <= xout) && (xout <= _Xhigh))) && (((yin <= _Ylow) && (_Yhigh <= yout)) || ((yout <= _Ylow) && (_Yhigh <= yin )) || ((_Ylow <= yout) && (yout <= _Yhigh))) && (((zin <= _Zlow) && (_Zhigh <= zout)) || ((zout <= _Zlow) && (_Zhigh <= zin )) || ((_Zlow <= zout) && (zout <= _Zhigh)))) ok = true; return ok; } void PndRadMapBoxMesh::makeHisto(const char* Orient, Double_t rotate, int Hbins, Double_t Hlow, Double_t Hhigh, int Vbins, Double_t Vlow, Double_t Vhigh, Double_t dlow, Double_t dhigh){ if(_verbose){ std::cout << "PndRadMapBoxMesh::makeHisto(" << Orient << ", " << rotate << ",\n"; std::cout << " " << Hbins << ", " << Hlow << ", " << Hhigh << ",\n"; std::cout << " " << Vbins << ", " << Vlow << ", " << Vhigh << ")\n"; } char C[32]; sprintf(C, "%s_Histo_%s_%g", _Name.Data(), Orient, _rotate); _MeshHisto = new TH2D(C, C, Hbins, Hlow, Hhigh, Vbins, Vlow, Vhigh); sprintf(C, "%s_StatHisto_%s_%g", _Name.Data(), Orient, _rotate); _StatHisto = new TH2I(C, C, Hbins, Hlow, Hhigh, Vbins, Vlow, Vhigh); if(_isSurfaceQuantity){ sprintf(C, "%s_EnergyHisto_%s_%g", _Name.Data(), Orient, _rotate); _EnergyHisto = new TH2D(C, C, 180, 0, 180, 10000, 0, 10); _plane = new PndRadMapPlane(dlow, _rotate, _orientation, _axis); } } Double_t PndRadMapBoxMesh::CalcFluence(FairRadMapPoint *p){ if(_verbose) std::cout << "PndRadMapBoxMesh::CalcFluence(FairRadMapPoint *p)" << std::endl; Double_t xx, yy, zz, dx = 0, dy = 0, dz =0; // float eff; std::vector CV; std::vector PV; //The Z should be between zlow and zhigh!! Transform(p->GetXOut(), p->GetYOut(), p->GetZOut(), xx, yy, zz); std::pair minp; bool OK, OK2; if(IsInside(p)){ OK = false; if(IsInside(xx, yy, zz)) OK = true; switch (_orientation){ case XY: case YX: dz = _Zlow; break; case XZ: case ZX: dy = _Ylow; break; case ZY: case YZ: dx = _Xlow; break; } TVector3 poststepv = TVector3(p->GetXOut(), p->GetYOut(), p->GetZOut()); TVector3 prestepv = TVector3(p->GetX(), p->GetY(), p->GetZ()); unsigned int prebin = _MeshHisto->FindBin(prestepv.X(), prestepv.Y()); unsigned int postbin = _MeshHisto->FindBin(poststepv.X(), poststepv.Y()); PndRadMapPlane plane; int prebinX = 0, prebinY = 0, prebinZ = 0; int postbinX = 0, postbinY = 0, postbinZ = 0; _MeshHisto->GetBinXYZ(postbin, postbinX, postbinY, postbinZ); _MeshHisto->GetBinXYZ(prebin, prebinX, prebinY, prebinZ); TVector3 tn(0, 0, 0); if(verbose) std::cout << "PndRadMapBoxMesh::CalcFluence - 1" << std::endl; // 1, only two planes: the narrowest bin lower and higher plane // 2, if the resulting two points in the range of the other two axis-limits // 3, start the stepping from one of those points // -------- 1 ------- TVector3 isl, lbeg, lend; double min_x, min_y, min_z, max_x, max_y, max_z; if(!OK){//if the 'out' is not in the given range //the 'in' and the 'out' is outside of the given range, so it can be that it will cross it if(verbose) std::cout << "PndRadMapBoxMesh::CalcFluence - 2" << std::endl; OK2 = false; min_x = _Xlow; min_y = _Ylow; min_z = _Zlow; max_x = _Xhigh; max_y = _Yhigh; max_z = _Zhigh; if(verbose) std::cout << "min_z: " << min_z << ' ' << max_z << ' ' << dz << std::endl; //the 'front' and 'back' XY planes plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, min_z+dz), TVector3(max_x+dx, min_y+dy, min_z+dz), TVector3(min_x+dx, max_y+dy, min_z+dz));//XY @ z=z_min isl = plane.LineIntersection(prestepv, poststepv); if(((min_x < isl.X()) && (isl.X() < max_x)) && ((min_y < isl.Y()) && (isl.Y() < max_y))){ OK2 = true; } if(verbose){ std::cout << "Xcross1: " << min_x << ' ' << isl.X() << ' ' << max_x << std::endl; std::cout << "Ycross1: " << min_y << ' ' << isl.Y() << ' ' << max_y << std::endl; std::cout << "=> " << OK2 << std::endl; } lbeg = isl; plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, max_z+dz), TVector3(max_x+dx, min_y+dy, max_z+dz), TVector3(min_x+dx, max_y+dy, max_z+dz));//XY @ z-z_max isl = plane.LineIntersection(prestepv, poststepv); if(((min_x < isl.X()) && (isl.X() < max_x)) && ((min_y < isl.Y()) && (isl.Y() < max_y))){ OK2 = OK2 && true; } if(verbose){ std::cout << "Xcross2: " << min_x << ' ' << isl.X() << ' ' << max_x << std::endl; std::cout << "Ycross2: " << min_y << ' ' << isl.Y() << ' ' << max_y << std::endl; std::cout << "=> " << OK2 << std::endl; } lend = isl; if(OK2){//the intersections are on the XY planes withing boarders at both z_min and z_max if(verbose) std::cout << "PndRadMapBoxMesh::CalcFluence - 3" << std::endl; prestepv = lbeg; poststepv = lend; prebin = _MeshHisto->FindBin(prestepv.X(), prestepv.Y()); postbin = _MeshHisto->FindBin(poststepv.X(), poststepv.Y()); _MeshHisto->GetBinXYZ(prebin, prebinX, prebinY, prebinZ); _MeshHisto->GetBinXYZ(postbin, postbinX, postbinY, postbinZ); }else{ if(verbose) std::cout << "PndRadMapBoxMesh::CalcFluence - 4" << std::endl; } }else{ // if(prebin != postbin){ int counter = 0; while((prebin != postbin) && !tn.Z()){ counter++; PV.clear(); // max // /---/ y| // /---/| | /z // || || | / // |/--|/ | / // ----- |/ // min +---------- // x // . corner1 corner2 corner3 //plane1 [(min, min, min), (min, min, max), (max, min, min)] ZX //plane2 [(min, min, min), (max, min, min), (min, max, min)] XY //plane3 [(min, min, min), (min, min, max), (min, max, min)] ZY //plane4 [(min, max, min), (min, max, max), (max, max, min)] ZX //plane5 [(min, min, max), (max, min, max), (min, max, max)] XY //plane6 [(max, min, min), (max, min, max), (max, max, min)] ZY // //plane4 [(max, max, max), (max, max, min), (min, max, max)] XZ // //plane5 [(max, max, max), (min, max, max), (max, min, max)] // //plane6 [(max, max, max), (max, max, min), (max, min, max)] min_x = (_MeshHisto->GetXaxis()->GetBinLowEdge(prebinX) < 0) ? 0 : _MeshHisto->GetXaxis()->GetBinLowEdge(prebinX); min_y = (_MeshHisto->GetYaxis()->GetBinLowEdge(prebinY) < 0) ? 0 : _MeshHisto->GetYaxis()->GetBinLowEdge(prebinY); min_z = (_MeshHisto->GetZaxis()->GetBinLowEdge(prebinZ) < 0) ? 0 : _MeshHisto->GetZaxis()->GetBinLowEdge(prebinZ); max_x = (_MeshHisto->GetXaxis()->GetBinUpEdge(prebinX) < 0) ? 0 : _MeshHisto->GetXaxis()->GetBinUpEdge(prebinX); max_y = (_MeshHisto->GetYaxis()->GetBinUpEdge(prebinY) < 0) ? 0 : _MeshHisto->GetYaxis()->GetBinUpEdge(prebinY); max_z = (_MeshHisto->GetZaxis()->GetBinUpEdge(prebinZ) < 0) ? 0 : _MeshHisto->GetZaxis()->GetBinUpEdge(prebinZ); max_z += _MeshHisto->GetZaxis()->GetBinWidth(1); if(verbose){ std::cout << counter << " - Corners [[" << min_x << ", " << min_y << ", " << min_z << "] - "; std::cout << "[" << max_x << ", " << max_y << ", " << max_z << "]]\n"; } //plane1 plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, min_z+dz), TVector3(min_x+dx, min_y+dy, max_z+dz), TVector3(max_x+dx, min_y+dy, min_z+dz)); PV.push_back(plane); //plane2 plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, min_z+dz), TVector3(max_x+dx, min_y+dy, min_z+dz), TVector3(min_x+dx, max_y+dy, min_z+dz)); PV.push_back(plane); //plane3 plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, min_z+dz), TVector3(min_x+dx, min_y+dy, max_z+dz), TVector3(min_x+dx, max_y+dy, min_z+dz)); PV.push_back(plane); //plane4 plane = PndRadMapPlane(TVector3(min_x+dx, max_y+dy, min_z+dz), TVector3(min_x+dx, max_y+dy, max_z+dz), TVector3(max_x+dx, max_y+dy, min_z+dz)); PV.push_back(plane); //plane5 plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, max_z+dz), TVector3(max_x+dx, min_y+dy, max_z+dz), TVector3(min_x+dx, max_y+dy, max_z+dz)); PV.push_back(plane); //plane6 plane = PndRadMapPlane(TVector3(max_x+dx, min_y+dy, min_z+dz), TVector3(max_x+dx, min_y+dy, max_z+dz), TVector3(max_x+dx, max_y+dy, min_z+dz)); PV.push_back(plane); if(verbose) std::cout << "Start\n[" << prebin << " {" << prestepv.X() << ' ' << prestepv.Y() << ' ' << prestepv.Z() << "}(" << prebinX << ' ' << prebinY << ' ' << prebinZ << ")\n" << ' ' << postbin << " {" << poststepv.X() << ' ' << poststepv.Y() << ' ' << poststepv.Z() << "}(" << postbinX << ' ' << postbinY << ' ' << postbinZ << ")]\n\n"; minp = std::make_pair(-1, 111111);//(which plane, what is the disrance) //calculate the intersection with every each planes isl = PV.at(0).LineIntersection(prestepv, poststepv); if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv))) if(minp.second > (isl-prestepv).Mag()) minp = std::make_pair(0, (isl-prestepv).Mag()); if(verbose) std::cout << "XZ (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") " << (poststepv-prestepv).Mag() << ", " << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") " << (isl-prestepv).Mag() << ' ' << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' ' << (minp.second > (isl-prestepv).Mag()) << ' ' << "[" << minp.first << ", " << minp.second << std::endl; if(!OK2){ isl = PV.at(1).LineIntersection(prestepv, poststepv); if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv))) if(minp.second > (isl-prestepv).Mag()) minp = std::make_pair(1, (isl-prestepv).Mag()); if(verbose) std::cout << "XY (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") " << (poststepv-prestepv).Mag() << ", " << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") " << (isl-prestepv).Mag() << ' ' << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' ' << (minp.second > (isl-prestepv).Mag()) << ' ' << "[" << minp.first << ", " << minp.second << "]" << std::endl; } isl = PV.at(2).LineIntersection(prestepv, poststepv); if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv))) if(minp.second > (isl-prestepv).Mag()) minp = std::make_pair(2, (isl-prestepv).Mag()); if(verbose) std::cout << "YZ (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") " << (poststepv-prestepv).Mag() << ", " << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") " << (isl-prestepv).Mag() << ' ' << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' ' << (minp.second > (isl-prestepv).Mag()) << ' ' << "[" << minp.first << ", " << minp.second << "]" << std::endl; isl = PV.at(3).LineIntersection(prestepv, poststepv); if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv))) if(minp.second > (isl-prestepv).Mag()) minp = std::make_pair(3, (isl-prestepv).Mag()); if(verbose) std::cout << "XZ (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") " << (poststepv-prestepv).Mag() << ", " << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") " << (isl-prestepv).Mag() << ' ' << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' ' << (minp.second > (isl-prestepv).Mag()) << ' ' << "[" << minp.first << ", " << minp.second << "]" << std::endl; if(!OK2){ isl = PV.at(4).LineIntersection(prestepv, poststepv); if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv))) if(minp.second > (isl-prestepv).Mag()) minp = std::make_pair(4, (isl-prestepv).Mag()); if(verbose) std::cout << "XY (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") " << (poststepv-prestepv).Mag() << ", " << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") " << (isl-prestepv).Mag() << ' ' << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' ' << (minp.second > (isl-prestepv).Mag()) << ' ' << "[" << minp.first << ", " << minp.second << "]" << std::endl; } isl = PV.at(5).LineIntersection(prestepv, poststepv); if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv))) if(minp.second > (isl-prestepv).Mag()) minp = std::make_pair(5, (isl-prestepv).Mag()); if(verbose) std::cout << "YZ (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") " << (poststepv-prestepv).Mag() << ", " << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") " << (isl-prestepv).Mag() << ' ' << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' ' << (minp.second > (isl-prestepv).Mag()) << ' ' << "[" << minp.first << ", " << minp.second << "]" << std::endl; if(verbose) std::cout << "\nmomentum: " << p->GetPx() << ' ' << p->GetPy() << ' ' << p->GetPz() << std::endl; //minp stores the smallest possible step if(0 <= minp.first){ if(verbose){ std::cout << "\nmomentum: " << p->GetPx() << ' ' << p->GetPy() << ' ' << p->GetPz() << std::endl; std::cout << "[(----)] " << minp.first << ' ' << minp.second << std::endl; std::cout << "normal: " << PV.at(minp.first).Normal().X() << ' ' << PV.at(minp.first).Normal().Y() << ' ' << PV.at(minp.first).Normal().Z() << std::endl; } isl = PV.at(minp.first).LineIntersection(prestepv, poststepv); Double_t step = (isl-prestepv).Mag(); prestepv = isl; _MeshHisto->Fill(prestepv.X(), prestepv.Y(), step); _StatHisto->Fill(prestepv.X(), prestepv.Y()); tn = PV.at(minp.first).Normal(); float ds = 0; if(tn.X()){ ds = _MeshHisto->GetXaxis()->GetBinWidth(1)/10.; if(verbose) std::cout << "Xprebin: " << prebin << " ---> "; if(p->GetPx() > 0){ prebin = _MeshHisto->FindBin(isl.X() + ds, isl.Y()); if(verbose) std::cout << "(" << isl.X() << " ---> " << isl.X() + ds << ' ' << p->GetPx() << " ) "; }else{ prebin = _MeshHisto->FindBin(isl.X() - ds, isl.Y()); if(verbose) std::cout << "(" << isl.X() << " ---> " << isl.X() - ds << ' ' << p->GetPx() << " ) "; } if(verbose) std::cout << prebin << std::endl; } if(tn.Y()){ ds = _MeshHisto->GetYaxis()->GetBinWidth(1)/10.; if(verbose) std::cout << "Yprebin: " << prebin << " ---> "; if(p->GetPy() > 0){ prebin = _MeshHisto->FindBin(isl.X(), isl.Y() + ds); if(verbose) std::cout << "1(" << isl.Y() << " ---> " << isl.Y() + ds << " ) "; }else{ prebin = _MeshHisto->FindBin(isl.X(), isl.Y()-ds); if(verbose) std::cout << "1(" << isl.Y() << " ---> " << isl.Y() - ds << " ) "; } if(verbose) std::cout << prebin << std::endl; } _MeshHisto->GetBinXYZ(prebin, prebinX, prebinY, prebinZ); if(verbose) std::cout << "Stop\n[" << prebin << " {" << prestepv.X() << ' ' << prestepv.Y() << ' ' << prestepv.Z() << "}(" << prebinX << ' ' << prebinY << ' ' << prebinZ << ")\n" << ' ' << postbin << " {" << poststepv.X() << ' ' << poststepv.Y() << ' ' << poststepv.Z() << "}(" << postbinX << ' ' << postbinY << ' ' << postbinZ << ")]\n\n"; if(IsInside(p->GetX(), p->GetY(), p->GetZ())){ if(verbose) std::cout << " Bingo\n"; } } } } } return minp.second; } void InvMatVecProd(TMatrixD mat, TVector3 vec, TVector3& res){ TMatrixD imat = mat.Invert(); double val; double resl[3]; for(int i = 0; i < 3; i++){ val = 0; for(int j = 0; j < 3; j++){ val+=imat[i][j]*vec[j]; } resl[i] = val; } res = TVector3(resl[0], resl[1], resl[2]); }