///////////////////////////////////////////////////////////// // // CbmEmc // // Filler of CbmEmcPoint // // Created 14/08/06 by S.Spataro // /////////////////////////////////////////////////////////////// #include "iostream" #include "TClonesArray.h" #include "TGeoMCGeometry.h" #include "TGeoManager.h" #include "TLorentzVector.h" #include "TParticle.h" #include "TVirtualMC.h" #include "TGeoArb8.h" #include "CbmGeoInterface.h" #include "CbmGeoLoader.h" #include "CbmGeoRootBuilder.h" #include "CbmEmc.h" #include "CbmEmcPoint.h" #include "CbmEmcReader.h" #include "CbmRootManager.h" #include "CbmVolume.h" #include "CbmGeoMedia.h" // add on for debug #include "CbmGeoG3Builder.h" #include "CbmRuntimeDb.h" #include "TObjArray.h" #include "CbmRun.h" #include "EmcStructure.h" // ----- Default constructor ------------------------------------------- CbmEmc::CbmEmc() { fEmcCollection = new TClonesArray("CbmEmcPoint"); fPosIndex = 0; fEventID=-1; bIsFastFsc = kFALSE; } // ------------------------------------------------------------------------- // ----- Standard constructor ------------------------------------------ CbmEmc::CbmEmc(const char* name, Bool_t active, Bool_t fast) : CbmDetector(name, active) { fEmcCollection = new TClonesArray("CbmEmcPoint"); fPosIndex = 0; fEventID=-1; bIsFastFsc = fast; } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- CbmEmc::~CbmEmc() { if (fEmcCollection) { fEmcCollection->Delete(); delete fEmcCollection; } } // ------------------------------------------------------------------------- // ----- Public method Intialize --------------------------------------- void CbmEmc::Initialize() { // Init function CbmDetector::Initialize(); CbmRun* sim = CbmRun::Instance(); CbmRuntimeDb* rtdb=sim->GetRuntimeDb(); par=(CbmGeoEmcPar*)(rtdb->getContainer("CbmGeoEmcPar")); } // ------------------------------------------------------------------------- void CbmEmc::BeginEvent(){ // Begin of the event } // ----- Public method ProcessHits -------------------------------------- Bool_t CbmEmc::ProcessHits(CbmVolume* vol) { // Increment number of emc points for TParticle if (gMC->IsTrackEntering()) { Int_t points = gMC->GetStack()->GetCurrentTrack()->GetMother(1); Int_t nEmcPoints = (points & (15<<16)) >> 16; nEmcPoints ++; if (nEmcPoints <= 15) { points = ( points & ( ~ (15<<16) ) ) | (nEmcPoints << 16); gMC->GetStack()->GetCurrentTrack()->SetMother(1,points); } } if (gMC->Edep()<=0) return kTRUE; // skip all the points which have no energy loss (i.e. Entering) TString nam = gMC->CurrentVolName(); fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); // trk ID fEventID = gMC->CurrentEvent(); fELoss = gMC->Edep(); fLength = gMC->TrackLength(); fTime = gMC->TrackTime(); gMC->TrackPosition(fPos); // cm gMC->TrackMomentum(fMom); // GeV Int_t copyNo = -1, id = -1; Int_t nMod = -1, nRow = -1, nCrys = -1; if (nam.BeginsWith("emc")) { sscanf(nam,"emc%dr%dc%d", &nMod, &nRow, &nCrys); // Crys 1-5000; copyNo 1-20; nRow 1-100, nMod 1-3 if ((nMod==1) || (nMod==2)) id = gMC->CurrentVolOffID(2,copyNo); if ((nMod==3) || (nMod==4)) id = gMC->CurrentVolOffID(1,copyNo); // 1 -because the pad stays inside flayer4 (Emc4), so only "1" as inheritance. // In barrel part one pad stays into Emc1 which stays inside Emc12 (and after Emc12 is // copied and rotated -> the inheritance level is "2" } else if (bIsFastFsc) { nMod = 5; nRow = abs((Int_t)(fPos.X()/11.))+1; nCrys = abs((Int_t)(fPos.Y()/11.))+1; if ((fPos.X()<0.) && (fPos.Y()>0.)) copyNo = 1; if ((fPos.X()<0.) && (fPos.Y()<0.)) copyNo = 2; if ((fPos.X()>0.) && (fPos.Y()<0.)) copyNo = 3; if ((fPos.X()>0.) && (fPos.Y()>0.)) copyNo = 4; } else if (!bIsFastFsc) { nam = gMC->CurrentVolOffName(2); sscanf(nam,"emc%dr%dc%d", &nMod, &nRow, &nCrys); if (nMod==5) id = gMC->CurrentVolOffID(3,copyNo); } fVolumeID = nMod*100000000 + nRow*1000000 + copyNo*10000 + nCrys; TVector3 pos(fPos.X(), fPos.Y(), fPos.Z()); AddHit(fTrackID, fVolumeID, fEventID, TVector3(fPos.X(), fPos.Y(), fPos.Z()), TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), fTime, fLength, fELoss, nMod, nRow, nCrys, copyNo); ResetParameters(); return kTRUE; } // ---------------------------------------------------------------------------- // ----- Public method EndOfEvent ----------------------------------------- void CbmEmc::EndOfEvent() { if (fVerboseLevel) Print(); Reset(); } // ---------------------------------------------------------------------------- // ----- Public method Register ------------------------------------------- void CbmEmc::Register() { CbmRootManager::Instance()->Register("EmcPoint","Emc", fEmcCollection, kTRUE); } // ---------------------------------------------------------------------------- // ----- Public method GetCollection -------------------------------------- TClonesArray* CbmEmc::GetCollection(Int_t iColl) const { if (iColl == 0) return fEmcCollection; return NULL; } // ---------------------------------------------------------------------------- // ----- Public method Print ---------------------------------------------- void CbmEmc::Print() const { Int_t nHits = fEmcCollection->GetEntriesFast(); cout << "-I- CbmEmc: " << nHits << " points registered in this event." << endl; if (fVerboseLevel>1) for (Int_t i=0; iPrint(); } // ---------------------------------------------------------------------------- // ----- Public method Reset ---------------------------------------------- void CbmEmc::Reset() { fEmcCollection->Clear(); fPosIndex = 0; } // ---------------------------------------------------------------------------- // guarda in CbmRootManager::CopyClones // ----- Public method CopyClones ----------------------------------------- void CbmEmc::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset ) { Int_t nEntries = cl1->GetEntriesFast(); //cout << "-I- CbmEmc: " << nEntries << " entries to add." << endl; TClonesArray& clref = *cl2; CbmEmcPoint* oldpoint = NULL; for (Int_t i=0; iAt(i); Int_t index = oldpoint->GetTrackID() + offset; oldpoint->SetTrackID(index); new (clref[fPosIndex]) CbmEmcPoint(*oldpoint); fPosIndex++; } cout << " -I- CbmEmc: " << cl2->GetEntriesFast() << " merged entries." << endl; } // ---------------------------------------------------------------------------- // ----- Public method ConstructGeometry ---------------------------------- void CbmEmc::ConstructGeometry() { // Definition of materials CbmGeoLoader*geoLoad = CbmGeoLoader::Instance(); CbmGeoInterface *geoFace = geoLoad->getGeoInterface(); CbmGeoMedia *Media = geoFace->getMedia(); CbmGeoBuilder *geobuild=geoLoad->getGeoBuilder(); CbmGeoMedium *CbmMediumPb = Media->getMedium("lead"); CbmGeoMedium *CbmMediumPWO = Media->getMedium("PWO"); CbmGeoMedium *CbmMediumFsc = Media->getMedium("FscScint"); Int_t nmedPb=geobuild->createMedium(CbmMediumPb); Int_t nmedPWO=geobuild->createMedium(CbmMediumPWO); Int_t nmedFsc=geobuild->createMedium(CbmMediumFsc); TGeoVolume *flayer1 = new TGeoVolumeAssembly("EmcLayer1"); TGeoVolume *flayer2 = new TGeoVolumeAssembly("EmcLayer2"); TGeoVolume *flayer3 = new TGeoVolumeAssembly("Emc3"); TGeoVolume *flayer4 = new TGeoVolumeAssembly("Emc4"); TGeoVolume *flayer5 = new TGeoVolumeAssembly("Fsc"); Bool_t bIsModuleOn[5] = {kFALSE, kFALSE, kFALSE, kFALSE, kFALSE}; Bool_t isFirst = kTRUE; CbmEmcReader read(GetGeometryFileName() ); for(Int_t module=read.GetMinModules(); module<=read.GetMaxModules(); module++) { cout << "Emc module = " << module; if (module==5 && bIsFastFsc) { cout << "\t FAST" << endl << "******** " << endl; continue; } cout << endl << "******** " << endl; for(Int_t row=read.GetMinRows(module); row<=read.GetMaxRows(module); row++){ for(Int_t crystal=read.GetMinCrystals(module,row); crystal<=read.GetMaxCrystals(module,row); crystal++) { Text_t buffer[30]; sprintf(buffer,"emc0%dr%dc%d",module, row, crystal); DataG4 data = read.GetData(module,row,crystal); if (data.module==-1) continue; //if the pad is not present, do not create geometry if (module<5) { // Construction of target spectrometer geometry TGeoTrap *trap = new TGeoTrap(data.pDz/10., data.pTheta, data.pPhi, data.pDy1/10., data.pDx1/10., data.pDx2/10., data.pAlp1, data.pDy2/10., data.pDx3/10., data.pDx4/10., data.pAlp2); TGeoVolume *volume; volume = new TGeoVolume(buffer, trap, gGeoManager->GetMedium("PWO")); TGeoRotation rot; rot.RotateZ(data.tau); rot.RotateY(data.theta); rot.RotateZ(data.phi); if(module ==1) flayer1->AddNode(volume,0, new TGeoCombiTrans(data.posX/10., data.posY/10., data.posZ/10.+3.7, new TGeoRotation (rot))); // shift of 37 mm with respect to interaction point if(module ==2) flayer2->AddNode(volume,0, new TGeoCombiTrans(data.posX/10., data.posY/10., data.posZ/10.+3.7, new TGeoRotation (rot))); // shift of 37 mm with respect to interaction point if(module ==3) flayer3->AddNode(volume,0, new TGeoCombiTrans(data.posX/10., data.posY/10., data.posZ/10., new TGeoRotation (rot))); if(module ==4) flayer4->AddNode(volume,0, new TGeoCombiTrans(data.posX/10., data.posY/10., data.posZ/10., new TGeoRotation (rot))); bIsModuleOn[module-1] = kTRUE; AddSensitiveVolume(volume); } if (module==5 && !bIsFastFsc) { // Construction of forward spectrometer geometry // For the forward calorimeter the box geometry was written giving different meaning to the DataG4 members: // Dx1 Dy1 : widths in XY of the absorber // Dx2 Dy2 : widths of XYthe scintillator // pAlp1 pAlp2 : withs on Y of absorber/scintillator // pDz : space in between absorber and scintillator TGeoVolume *volAbs, *volSci; TGeoRotation rot; TGeoVolume *volume = new TGeoVolumeAssembly(buffer); TGeoBBox *absorber = new TGeoBBox(data.pDx1/10., data.pDy1/10., data.pAlp1/10.); TGeoBBox *scintillator = new TGeoBBox(data.pDx2/10., data.pDy2/10., data.pAlp2/10.); volAbs = new TGeoVolume("FscAbsorber", absorber, gGeoManager->GetMedium("lead")); volSci = new TGeoVolume("FscScintillator", scintillator, gGeoManager->GetMedium("FscScint")); TGeoVolume *ffsclay = new TGeoVolumeAssembly("FscLayer"); ffsclay->AddNode(volAbs, 0, new TGeoCombiTrans(0., 0., data.pAlp1/10., new TGeoRotation (rot))); ffsclay->AddNode(volSci, 0, new TGeoCombiTrans(0., 0., 2*data.pAlp1/10+data.pDz/10.+data.pAlp2/10., new TGeoRotation (rot))); AddSensitiveVolume(volSci); for (Int_t ll=1; ll<=378; ll++) { volume->AddNode(ffsclay,ll, new TGeoCombiTrans(0., 0., ll*(2.*data.pDz+2.*data.pAlp1+2.*data.pAlp2)/10., new TGeoRotation (rot))); } flayer5->AddNode(volume, 0, new TGeoCombiTrans(data.posX/10., data.posY/10., data.posZ/10., new TGeoRotation (rot))); bIsModuleOn[module-1] = kTRUE; } } } } // Fast construction of the fsc geometry if (bIsFastFsc) { TGeoVolume *volAbs1, *volSci1, *volAbs2, *volSci2, *volAbs3, *volSci3; TGeoRotation rot; TGeoVolume *volume = new TGeoVolumeAssembly("FscBox"); Int_t padX = 28, padY = 14, holeX = 4, holeY = 2, padX1 = 11, padX2 = 13; Float_t pDx1 = 11., pDx2 = 11., pDy1 = 11., pDy2 = 11., pAlp1 = 0.0275, pAlp2 = 0.15, pDz = 0.00375, posZ = 690.; TGeoBBox *absorber1 = new TGeoBBox(pDx1*padX/2. , pDy1*(padY-holeY)/4., pAlp1/2.); TGeoBBox *scintillator1 = new TGeoBBox(pDx2*padX/2. , pDy2*(padY-holeY)/4., pAlp2/2.); TGeoBBox *absorber2 = new TGeoBBox(pDx1*padX1/2., pDy1*holeY/4. , pAlp1/2.); TGeoBBox *scintillator2 = new TGeoBBox(pDx2*padX1/2., pDy2*holeY/4. , pAlp2/2.); TGeoBBox *absorber3 = new TGeoBBox(pDx1*padX2/2., pDy1*holeY/4. , pAlp1/2.); TGeoBBox *scintillator3 = new TGeoBBox(pDx2*padX2/2., pDy2*holeY/4. , pAlp2/2.); volAbs1 = new TGeoVolume("FscAbsorber1", absorber1, gGeoManager->GetMedium("lead")); volSci1 = new TGeoVolume("FscScintillator1", scintillator1, gGeoManager->GetMedium("FscScint")); volAbs2 = new TGeoVolume("FscAbsorber2", absorber2, gGeoManager->GetMedium("lead")); volSci2 = new TGeoVolume("FscScintillator2", scintillator2, gGeoManager->GetMedium("FscScint")); volAbs3 = new TGeoVolume("FscAbsorber3", absorber3, gGeoManager->GetMedium("lead")); volSci3 = new TGeoVolume("FscScintillator3", scintillator3, gGeoManager->GetMedium("FscScint")); TGeoVolume *volSci1X = gGeoManager->Division("FscScintillator1X" ,"FscScintillator1" , 1, padX , -pDx2*padX/2. , pDx2); TGeoVolume *volSci1XY = gGeoManager->Division("FscScintillator1XY","FscScintillator1X", 2, (padY-holeY)/2, -pDy2*(padY-holeY)/4., pDy2); TGeoVolume *volSci2X = gGeoManager->Division("FscScintillator2X" ,"FscScintillator2" , 1, padX1 , -pDx2*padX1/2. , pDx2); TGeoVolume *volSci2XY = gGeoManager->Division("FscScintillator2XY","FscScintillator2X", 2, holeY/2 , -pDy2*holeY/4. , pDy2); TGeoVolume *volSci3X = gGeoManager->Division("FscScintillator3X" ,"FscScintillator3" , 1, padX2 , -pDx2*padX2/2. , pDx2); TGeoVolume *volSci3XY = gGeoManager->Division("FscScintillator3XY","FscScintillator3X", 2, holeY/2 , -pDy2*holeY/4. , pDy2); TGeoVolume *ffsclay = new TGeoVolumeAssembly("FscLayer"); ffsclay->AddNode(volAbs1, 0, new TGeoCombiTrans(0. , pDy1*(padY+holeY)/4., pAlp1/2. , new TGeoRotation (rot))); ffsclay->AddNode(volSci1, 0, new TGeoCombiTrans(0. , pDy2*(padY+holeY)/4., pAlp1+pDz+pAlp2/2., new TGeoRotation (rot))); ffsclay->AddNode(volAbs2, 0, new TGeoCombiTrans(pDx1*(-padX+padX1)/2., pDy1*holeY/4. , pAlp1/2. , new TGeoRotation (rot))); ffsclay->AddNode(volSci2, 0, new TGeoCombiTrans(pDx2*(-padX+padX1)/2., pDy2*holeY/4 , pAlp1+pDz+pAlp2/2., new TGeoRotation (rot))); ffsclay->AddNode(volAbs3, 0, new TGeoCombiTrans(pDx1*(padX-padX2)/2. , pDy1*holeY/4. , pAlp1/2. , new TGeoRotation (rot))); ffsclay->AddNode(volSci3, 0, new TGeoCombiTrans(pDx2*(padX-padX2)/2. , pDy2*holeY/4. , pAlp1+pDz+pAlp2/2., new TGeoRotation (rot))); AddSensitiveVolume(volSci1XY); AddSensitiveVolume(volSci2XY); AddSensitiveVolume(volSci3XY); for (Int_t ll=1; ll<=300; ll++) { volume->AddNode(ffsclay,ll, new TGeoCombiTrans(0., 0., ll*(2.*pDz+pAlp1+pAlp2), new TGeoRotation (rot))); } flayer5->AddNode(volume, 0, new TGeoCombiTrans(0., 0., posZ, new TGeoRotation (rot))); bIsModuleOn[4] = kTRUE; } TGeoVolume *flayer12 = new TGeoVolumeAssembly("Emc12"); if (bIsModuleOn[0]) flayer12->AddNode(flayer1,0, new TGeoCombiTrans(0., 0., 0., new TGeoRotation(0))); if (bIsModuleOn[1]) flayer12->AddNode(flayer2,0, new TGeoCombiTrans(0., 0., 0., new TGeoRotation(0))); TString vname = "cave"; vname = vname.Strip(); TGeoVolume* vcave = gGeoManager->FindVolumeFast(vname.Data()); if (bIsModuleOn[0] || bIsModuleOn[1]) vcave->AddNode(flayer12, 1); if (bIsModuleOn[2]) vcave->AddNode(flayer3, 1); if (bIsModuleOn[3]) vcave->AddNode(flayer4, 1); if (bIsModuleOn[4]) vcave->AddNode(flayer5, 1); // 15 copies for barrel part of EMC (1st copy exists in emc_module12345.dat) if (bIsModuleOn[0] || bIsModuleOn[1]) for (Int_t n=1;n<=15;n++){ TGeoRotation rot1; rot1.RotateZ(22.5*n); vcave->AddNode(flayer12, n+1,new TGeoCombiTrans(0., 0., 0., new TGeoRotation (rot1)) ); } // 3 copies for forward endcap of EMC (1st copy exists in emc_module12345.dat) if (bIsModuleOn[2]) for (Int_t n=1;n<=3;n++){ TGeoCombiTrans reflection; if (n==1) reflection.ReflectY(1); if (n==2) {reflection.ReflectY(1); reflection.ReflectX(1);} if (n==3) reflection.ReflectX(1); vcave->AddNode(flayer3, n+1, new TGeoCombiTrans(reflection)); } // 3 copies for backward endcap of EMC (1st copy exists in emc_module12345.dat) if (bIsModuleOn[3]) for (Int_t n=1;n<=3;n++){ TGeoCombiTrans reflection; if (n==1) reflection.ReflectY(1); if (n==2) {reflection.ReflectY(1); reflection.ReflectX(1);} if (n==3) reflection.ReflectX(1); vcave->AddNode(flayer4, n+1, new TGeoCombiTrans(reflection)); } // 3 copies for forward calorimeter (1st copy exists in emc_module12345.dat) if (bIsModuleOn[4]) { if (bIsFastFsc) { TGeoCombiTrans reflection; reflection.ReflectY(1); vcave->AddNode(flayer5, 2, new TGeoCombiTrans(reflection)); } else for (Int_t n=1;n<=3;n++){ TGeoCombiTrans reflection; if (n==1) reflection.ReflectY(1); if (n==2) {reflection.ReflectY(1); reflection.ReflectX(1);} if (n==3) reflection.ReflectX(1); vcave->AddNode(flayer5, n+1, new TGeoCombiTrans(reflection)); } } } // ----- Private method AddHit -------------------------------------------- CbmEmcPoint* CbmEmc::AddHit(Int_t trackID, Int_t detID, Int_t evtID, TVector3 pos, TVector3 mom, Double_t time, Double_t length, Double_t eLoss, Short_t mod, Short_t row, Short_t crys, Short_t copy) { TClonesArray& clref = *fEmcCollection; Int_t size = clref.GetEntriesFast(); if (fVerboseLevel>1) cout << "-I- CbmEmc: Adding Point at IN (" << pos.X() << ", " << pos.Y() << ", " << pos.Z() << ") cm, detector " << detID << ", evt " << evtID << ", track " << trackID <<", energy loss " << eLoss*1e06 << " keV, module " << mod << " row " << row << " crystal " << crys << " copy " << copy << endl; return new(clref[size]) CbmEmcPoint(trackID, detID, evtID, pos, mom, time, length, eLoss, mod, row, crys, copy); } inline void CbmEmc::ResetParameters() { fTrackID = -999; fVolumeID = -999; fEventID = -999; fPos.SetXYZT(0., 0., 0., 0.); fMom.SetXYZT(0., 0., 0., 0.) ; fTime = -999; fLength = -999; fELoss = -999; } // ---- ClassImp(CbmEmc)