/* $Id: CbmRichProjectionProducer.cxx,v 1.9 2006/06/22 08:23:40 hoehne Exp $ */ /* History of cvs commits: * * $Log: CbmRichProjectionProducer.cxx,v $ * Revision 1.9 2006/06/22 08:23:40 hoehne * bug found in filling some variables * * Revision 1.8 2006/06/21 13:54:23 hoehne * flag added in order to choose between projection from imaginary z-plane or mirror * * Revision 1.7 2006/01/26 09:56:36 hoehne * updates due to changes in extrapolation routines: acts now on "RichTrackZ" * * Revision 1.6 2005/12/19 19:04:31 friese * New CbmTask design * * Revision 1.5 2005/12/08 15:13:04 turany * change GetRegistered, ActivateBranch and CheckActivatedBranch * to GetObject, also add a flag to the arguments of Register which control saving * to file(kTRUE) or only in memory(kFALSE) * * Revision 1.4 2005/11/30 15:38:17 hoehne * hardcoded parameters removed * * Revision 1.3 2005/07/18 15:03:09 kharlov * Correct projection counter * * Revision 1.2 2005/07/18 10:11:29 kharlov * Typos are corrected * * Revision 1.1 2005/07/18 09:53:44 kharlov * Track projection from imaginary plane onto a photodetector * */ // ------------------------------------------------------------------------- // ----- CbmRichProjectionProducer source file ----- // ----- Created 23/06/05 by P.Stolpovsky(P.Stolpovsky at gsi.de) ----- // ----- Project track by strate line from imagiry plane to the mirror ----- // ----- and reflect it to the photodetector plane ----- // ------------------------------------------------------------------------- #include "CbmRichProjectionProducer.h" #include "CbmRichRing.h" #include "CbmRootManager.h" #include "CbmMCTrack.h" #include "CbmRichPoint.h" #include "CbmTrackParam.h" #include "CbmRichRing.h" #include "CbmRuntimeDb.h" #include "CbmRun.h" #include "CbmGeoNode.h" #include "CbmGeoRichPar.h" #include "CbmGeoTransform.h" #include "CbmGeoVector.h" #include "TObjArray.h" #include "TVector3.h" #include "TLorentzVector.h" #include "TRandom.h" #include "TFormula.h" #include "iostream.h" #include "CbmRunAna.h" #include "TMath.h" #include "fstream.h" #include "iostream.h" #include "stdio.h" // ----- Default constructor ------------------------------------------- CbmRichProjectionProducer::CbmRichProjectionProducer() : CbmTask("RichProjectionProducer") { SetDefaultParameters(); } // ------------------------------------------------------------------------- // ----- Standard constructor ------------------------------------------- CbmRichProjectionProducer::CbmRichProjectionProducer(Int_t verbose, Int_t zflag) :CbmTask("RichProjectionProducer") { fVerbose = verbose; fZflag = zflag; // 1 - use imaginary plane; 2 - use mirror point for extrapolation } // ------------------------------------------------------------------------- // ------------------------------------------------------------------------- /** Constructor with name and title */ CbmRichProjectionProducer::CbmRichProjectionProducer(const char *name, const char *title) : CbmTask(name) {} // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- CbmRichProjectionProducer::~CbmRichProjectionProducer() { CbmRootManager *fManager =CbmRootManager::Instance(); fManager->Write(); } // ------------------------------------------------------------------------- // ------------------------------------------------------------------------- /** Private method SetDefaultParameters */ void CbmRichProjectionProducer::SetDefaultParameters() { fVerbose = 1; fZflag = 1; } // ------------------------------------------------------------------------- //----- Initialization of Parameter Containers ------------------------- void CbmRichProjectionProducer::SetParContainers() { CbmRunAna* sim = CbmRunAna::Instance(); CbmRuntimeDb* rtdb=sim->GetRuntimeDb(); par=(CbmGeoRichPar*)(rtdb->getContainer("CbmGeoRichPar")); } // ----- Initialization ------------------------------------------------ InitStatus CbmRichProjectionProducer::Init() { CbmRootManager* fManager = CbmRootManager::Instance(); fSensNodes = par->GetGeoSensitiveNodes(); fPassNodes = par->GetGeoPassiveNodes(); // par->printParams(); // fSensNodes->Dump(); // fPassNodes->Dump(); // get detector position: CbmGeoNode *det= (CbmGeoNode *) fSensNodes->FindObject("rich1d#1"); CbmGeoTransform* detTr=det->getLabTransform(); // detector position in labsystem CbmGeoVector detPosLab=detTr->getTranslation(); // ... in cm CbmGeoTransform detCen=det->getCenterPosition(); // center in Detector system CbmGeoVector detPosCen=detCen.getTranslation(); fDetZ = detPosLab.Z() + detPosCen.Z(); /** z coordinate of photodetector (Labsystem, cm) */ fDetY = detPosLab.Y() + detPosCen.Y(); /** y coordinate of photodetector (Labsystem, cm) */ fDetX = detPosLab.X() + detPosCen.X(); /** x coordinate of photodetector (Labsystem, cm) */ TArrayD *fdetA=det->getParameters(); // get other geometry parameters: width in x, width in y, thickness fDetWidthX = fdetA->At(0); fDetWidthY = fdetA->At(1); // for(Int_t i=0;iGetSize();i++) cout << "Array detector " << fdetA->At(i)<< endl; if (fVerbose) { cout << "---------------------- RICH Projection Producer ---------------------------------------" << endl; cout << " detector position in (x,y,z): " << fDetX << " " << fDetY << " " << fDetZ << endl; cout << " detector size in x and y: " << fDetWidthX << " " << fDetWidthY << endl; } // get mirror position: //CbmGeoNode *mir= (CbmGeoNode *) fPassNodes->FindObject("rich1mgl#1"); CbmGeoNode *mir= (CbmGeoNode *) fSensNodes->FindObject("rich1mgl#1"); CbmGeoTransform* mirTr=mir->getLabTransform(); // position of mirror center in labsystem CbmGeoVector mirPosLab=mirTr->getTranslation(); // ... in cm fZm = mirPosLab.Z(); fYm = mirPosLab.Y(); fXm = mirPosLab.X(); TArrayD *fmirA=mir->getParameters(); // get other geometry parameters: radius, fR = fmirA->At(0); // for(Int_t i=0;iGetSize();i++) cout << "Array mirror " << fmirA->At(i)<< endl; if (fVerbose){ cout << " mirror center (x,y,z): " << fXm << " " << fYm << " " << fZm << endl; cout << " mirror radius: " << fR << endl; } fEvent = 0; fNHits = 0; fListRICHImPlanePoint = (TClonesArray*)fManager->GetObject("RichTrackParamZ"); // imaginary plane in Rich detector for track extrapolation if (fVerbose && fZflag == 1) cout << " use tracks in imaginary plane for projection to photodetector plane" << endl; if (fVerbose && fZflag == 2) cout << " use tracks in RICH mirror for projection to photodetector plane" << endl; cout << "--------------------------------------------------------------------------------" << endl; if ( ! fListRICHImPlanePoint) { cout << "-W- CbmRichProjectionProducer::Init: No Rich Z-Point array!" << endl; } fListStack = (TClonesArray*)fManager->GetObject("MCTrack"); // for checking workability fProjectionTrackParam = new TClonesArray("CbmTrackParam"); // coordinates of projected tracks to photodetectors plane fManager->Register("RichProjection","RICH",fProjectionTrackParam, kTRUE); return kSUCCESS; } // ------------------------------------------------------------------------- // ----- Execution of Task --------------------------------------------- void CbmRichProjectionProducer::Exec(Option_t* option) { fEvent++; printf("\n\n=====> Event no. %d\n",fEvent); fProjectionTrackParam->Clear(); fNHits = 0; CbmTrackParam* point=NULL; // some default variables TMatrixFSym covMat(5); for(Int_t i=0;i<5;i++) for(Int_t j=0; j<=i; j++) covMat(i,j) = 0; covMat(0,0) = covMat(1,1) = covMat(2,2) = covMat(3,3) = covMat(4,4) = 1.e-4; Double_t rho1,rho2; Double_t xX; Double_t yY; Double_t zZ; Double_t normP2,RxP,dist; Double_t crossPx,crossPy,crossPz; if (fVerbose > 1) cout << "Number of tracks in Imaginary z-Plane: " << fListRICHImPlanePoint->GetEntriesFast() << endl; for(Int_t j=0; jGetEntriesFast(); j++) { point = (CbmTrackParam*)fListRICHImPlanePoint->At(j); new((*fProjectionTrackParam)[j]) CbmTrackParam(0.,0.,0.,0.,0.,0.,covMat); if (fVerbose > 1) cout << "Data in Imaginary z-Plane: " << j << " " << point->GetX() << " " << point->GetY() << " " << point->GetZ() << " " << point->GetTx() << " " << point->GetTy() << " " << point->GetQp() << endl; // check if Array was filled if (point->GetX() == 0 && point->GetY() == 0 && point->GetZ() == 0 && point->GetTx() == 0 && point->GetTy() ==0) continue; if (point->GetQp()==0) continue; rho1 = rho2 = 0.; xX=9999.; yY=9999.; zZ=0.; normP2 = RxP = dist = 0.; TVector3 startP, momP, crossP, centerP; // operate on ImPlane point //------------------------------------------------------------------------------------------------- if (fZflag ==1) { Double_t p = 1./TMath::Abs(point->GetQp()); Double_t pz = p/TMath::Sqrt(1+point->GetTx()*point->GetTx()+point->GetTy()*point->GetTy()); Double_t px = pz*point->GetTx(); Double_t py = pz*point->GetTy(); momP.SetXYZ(px,py,pz); point->Position(startP); if ((fYm*startP.y())<0) fYm = -fYm; // check that mirror center and startP are in same hemisphere if (fVerbose > 1) cout << " --> position: " << startP.x() << " " << startP.y() << " " << startP.z() << endl; if (fVerbose > 1) cout << " --> momentum: " << momP.x() << " " << momP.y() << " " << momP.z() << endl; // calculation of intersection of track with selected mirror // corresponds to calculation of intersection between a straight line and a sphere: // vector: r = startP - mirrorCenter // RxP = r*momP // normP2 = momP^2 // dist = r^2 - fR^2 // -> rho1 = (-RxP+sqrt(RxP^2-normP2*dist))/normP2 extrapolation factor for: // intersection point crossP = startP + rho1 * momP RxP=(momP.x()*(startP.x()-fXm)+momP.y()*(startP.y()-fYm)+momP.z()*(startP.z()-fZm)); normP2=(momP.x()*momP.x()+momP.y()*momP.y()+momP.z()*momP.z()); dist=(startP.x()*startP.x()+fXm*fXm+startP.y()*startP.y()+fYm*fYm+startP.z()*startP.z()+fZm*fZm-2*startP.x()*fXm-2*startP.y()*fYm-2*startP.z()*fZm-fR*fR); if (normP2!=0.) rho1=(-RxP+TMath::Sqrt(RxP*RxP-normP2*dist))/normP2; if (normP2 == 0) cout << " Error in track extrapolation: momentum = 0 " << endl; crossPx=startP.x()+rho1*momP.x(); crossPy=startP.y()+rho1*momP.y(); crossPz=startP.z()+rho1*momP.z(); crossP.SetXYZ(crossPx,crossPy,crossPz); // check if crosspoint with mirror and chosen mirrorcenter (y) are in same hemisphere // if not recalculate crossing point if ((fYm*crossP.y())<0) { fYm = -fYm; RxP=(momP.x()*(startP.x()-fXm)+momP.y()*(startP.y()-fYm)+momP.z()*(startP.z()-fZm)); normP2=(momP.x()*momP.x()+momP.y()*momP.y()+momP.z()*momP.z()); dist=(startP.x()*startP.x()+fXm*fXm+startP.y()*startP.y()+fYm*fYm+startP.z()*startP.z()+fZm*fZm-2*startP.x()*fXm-2*startP.y()*fYm-2*startP.z()*fZm-fR*fR); if (normP2!=0.) rho1=(-RxP+TMath::Sqrt(RxP*RxP-normP2*dist))/normP2; if (normP2 == 0) cout << " Error in track extrapolation: momentum = 0 " << endl; crossPx=startP.x()+rho1*momP.x(); crossPy=startP.y()+rho1*momP.y(); crossPz=startP.z()+rho1*momP.z(); crossP.SetXYZ(crossPx,crossPy,crossPz); } centerP.SetXYZ(fXm,fYm,fZm); // mirror center } // operate on Rich Mirror point //--------------------------------------------------------------------------------------------------------- if (fZflag ==2) { Double_t p = 1./TMath::Abs(point->GetQp()); Double_t pz = p/TMath::Sqrt(1+point->GetTx()*point->GetTx()+point->GetTy()*point->GetTy()); Double_t px = pz*point->GetTx(); Double_t py = pz*point->GetTy(); momP.SetXYZ(px,py,pz); point->Position(crossP); if ((fYm*crossP.y())<0) fYm = -fYm; // check that mirror center and crossP are in same hemisphere centerP.SetXYZ(fXm,fYm,fZm); // mirror center } // calculate normal on crosspoint with mirror TVector3 normP(crossP.x()-centerP.x(),crossP.y()-centerP.y(),crossP.z()-centerP.z()); normP=normP.Unit(); // check that normal has same z-direction as momentum if ((normP.z()*momP.z())<0.) normP = TVector3(-1.*normP.x(),-1.*normP.y(),-1.*normP.z()); // reflect track Double_t np=normP.x()*momP.x()+normP.y()*momP.y()+normP.z()*momP.z(); Double_t refX = 2*np*normP.x()-momP.x(); Double_t refY = 2*np*normP.y()-momP.y(); Double_t refZ = 2*np*normP.z()-momP.z(); // crosspoint whith photodetector plane: // calculate intersection between straight line and plane: // intersection point = crossP + rho2 * refl_track if (refZ!=0.) { rho2 = -1*(crossP.z() - fDetZ)/refZ; xX = crossP.x() + refX * rho2; yY = crossP.y() + refY * rho2; zZ = crossP.z() + refZ * rho2; if( xX > (fDetX-fDetWidthX) && xX < (fDetX+fDetWidthX)) { //check that crosspoint inside the plane if(TMath::Abs(yY) > (fDetY-fDetWidthY) && TMath::Abs(yY) < (fDetY+fDetWidthY)) { CbmTrackParam richtrack(xX,yY,zZ,0.,0.,0.,covMat); * (CbmTrackParam*)(fProjectionTrackParam->At(j)) = richtrack; if (fVerbose > 1) cout << "-I- RichProjectionProducer: track extrapolated to photodetector plane: iTrack " << j << endl; } } } } cout<<"projected hits:"<Clear(); fListRICHImPlanePoint->Clear(); fListStack->Clear(); } // ------------------------------------------------------------------------- ClassImp(CbmRichProjectionProducer)