// ------------------------------------------------------------------------- #include "iostream.h" #include "TClonesArray.h" #include "TGeoMCGeometry.h" #include "TGeoManager.h" #include "TLorentzVector.h" #include "TParticle.h" #include "TVirtualMC.h" #include "FairGeoInterface.h" #include "FairGeoLoader.h" #include "FairGeoNode.h" #include "CbmGeoPlane.h" #include "FairGeoRootBuilder.h" #include "PndStack.h" #include "CbmPlane.h" #include "CbmPlanePoint.h" #include "FairRootManager.h" #include "FairVolume.h" // add on for debug //#include "FairGeoG3Builder.h" #include "FairRuntimeDb.h" #include "CbmGeoPlanePar.h" #include "TObjArray.h" #include "FairRun.h" #include "TGeant3.h" // ----- Default constructor ------------------------------------------- CbmPlane::CbmPlane() { fPlane1Collection = new TClonesArray("CbmPlanePoint"); fPlane2Collection = new TClonesArray("CbmPlanePoint"); fPlane3Collection = new TClonesArray("CbmPlanePoint"); // fPlaneImPlaneCollection = new TClonesArray("CbmPlanePoint"); fPosIndex = 0; fpreflag = 0; fpostflag = 0; volDetector = 0; } // ------------------------------------------------------------------------- // ----- Standard constructor ------------------------------------------ CbmPlane::CbmPlane(const char* name, Bool_t active) : FairDetector(name, active) { fPlane1Collection = new TClonesArray("CbmPlanePoint"); fPlane2Collection = new TClonesArray("CbmPlanePoint"); fPlane3Collection = new TClonesArray("CbmPlanePoint"); // fPlaneImPlaneCollection = new TClonesArray("CbmPlanePoint"); fPosIndex = 0; fpreflag = 0; fpostflag = 0; volDetector = 0; fevent=-1; } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- CbmPlane::~CbmPlane() { if (fPlane1Collection) { fPlane1Collection->Delete(); delete fPlane1Collection; } } // ------------------------------------------------------------------------- // ----- Public method Intialize --------------------------------------- void CbmPlane::Initialize() { // controlla-- forse non tutto e' necessario. FairDetector::Initialize(); FairRun* sim = FairRun::Instance(); FairRuntimeDb* rtdb=sim->GetRuntimeDb(); CbmGeoPlanePar *par=(CbmGeoPlanePar*)(rtdb->getContainer("CbmGeoPlanePar")); TObjArray *fSensNodes = par->GetGeoSensitiveNodes(); // FairGeoNode *fm1= (FairGeoNode *) fSensNodes->FindObject("Plane0000"); // cout << "GET COPY NO: " << fm1->getCopyNo() << endl; // volDetector = fm1->getMCid(); } // ------------------------------------------------------------------------- void CbmPlane::BeginEvent(){ } // ----- Public method ProcessHits -------------------------------------- Bool_t CbmPlane::ProcessHits(FairVolume* vol) { Int_t pdgCode = gMC->TrackPid(); if(pdgCode==13){ fmass = gMC->TrackMass(); // mass (GeV) fcharge = gMC->TrackCharge(); //charge? fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); // trk ID Int_t copyNo; gMC->CurrentVolID(copyNo); // if entering if ( gMC->IsTrackEntering() ) { fpreflag = 1; fELoss = 0.; gMC->TrackPosition(fPosIn); // cm gMC->TrackMomentum(fMomIn); // GeV } // Sum energy loss for all steps in the active volume fELoss += gMC->Edep(); // Set additional parameters at exit of active volume. Create CbmPlanePoint. if ( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared() ) { fpostflag = 1; gMC->TrackPosition(fPosOut); // cm gMC->TrackMomentum(fMomOut); // GeV fLength = gMC->TrackLength(); fTime = gMC->TrackTime() * 1e9; // time // if (fELoss == 0. ) return kFALSE; Bool_t fmuch = kFALSE; // cout << vol->getMCid() << " " << gMC->CurrentVolPath() << endl; if(fmuch == kTRUE) { if(vol->getMCid()==8){ AddHit2(fTrackID, fVolumeID, TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()), TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()), TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()), TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()), fTime, fLength, fELoss, fcharge, fmass, fpreflag, fpostflag, "P1",fevent); ResetParameters(); }else if(vol->getMCid() ==3){ AddHit1(fTrackID, fVolumeID, TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()), TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()), TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()), TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()), fTime, fLength, fELoss, fcharge, fmass, fpreflag, fpostflag, "P2",fevent); ResetParameters(); } } else if(fmuch == kFALSE) { // 45 if(vol->getMCid()==3){ //36 //37 //2814 //1959 AddHit1(fTrackID, fVolumeID, TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()), TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()), TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()), TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()), fTime, fLength, fELoss, fcharge, fmass, fpreflag, fpostflag, "P1",fevent); ResetParameters(); } else if(vol->getMCid() ==4){ //37 // 38 //2815 //1960 AddHit2(fTrackID, fVolumeID, TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()), TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()), TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()), TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()), fTime, fLength, fELoss, fcharge, fmass, fpreflag, fpostflag, "P2",fevent); ResetParameters(); } else if(vol->getMCid() == 2){ //35 // 36 //2813 //1958 AddHit3(fTrackID, fVolumeID, TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()), TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()), TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()), TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()), fTime, fLength, fELoss, fcharge, fmass, fpreflag, fpostflag, "P3",fevent); ResetParameters(); } } } } return kTRUE; } // ---------------------------------------------------------------------------- // ----- Public method EndOfEvent ----------------------------------------- void CbmPlane::EndOfEvent() { if (fVerboseLevel) Print(); Reset(); if (fevent==10000-1) plane(); } // ---------------------------------------------------------------------------- void CbmPlane::plane() { // /* // GEANE stuff // cout<<"Calling GEANE..."<Gpcxyz(); //a simple test Gctmed_t *afGmed = gMC3->Gctmed(); //a simple test //cout<epsil<Ertrio(); Eropts_t *afEropts=gMC3->Eropts(); for(Int_t i=0;i<15;i++) afErtrio->errin[i]=0.; afErtrio->nepred=1; for(Int_t i=0;i<6;i++) afEropts->erpli[i]=0.; afEropts->erpli[1]=1.; afEropts->erpli[5]=1.; for(Int_t i=0;i<6;i++) afEropts->erplo[i]=afEropts->erpli[i]; afEropts->erplo[6]=1000.; afEropts->erplo[7]=0.; afEropts->erplo[8]=0.; TVector3 v1(afEropts->erplo[0],afEropts->erplo[1],afEropts->erplo[2]); TVector3 v2(afEropts->erplo[3],afEropts->erplo[4],afEropts->erplo[5]); TVector3 v3=v1.Cross(v2); afEropts->erplo[9]=v3(0); afEropts->erplo[10]=v3(1); afEropts->erplo[11]=v3(2); //ERTRAK Float_t x1[3]={0.,0.,0.}; Float_t p1[3]={1.0,0.,0.}; Float_t x2[3]; Float_t p2[3]; gMC3->Ertrak(x1,p1,x2,p2,5,"P"); ofstream f("plane.out"); f<erpout[0]<<"\t"<erpout[1]<<"\t"<erpout[2]<erplo[9]<<"\t"<erplo[10]<<"\t"<erplo[11]<leexac<leleng<leonly<leplan<lepoin<levolu<errout[i]<Register("Plane1Point","Plane1", fPlane1Collection, kTRUE); FairRootManager::Instance()->Register("Plane2Point","Plane2", fPlane2Collection, kTRUE); FairRootManager::Instance()->Register("Plane3Point","Plane3", fPlane3Collection, kTRUE); } // ---------------------------------------------------------------------------- // ----- Public method GetCollection -------------------------------------- TClonesArray* CbmPlane::GetCollection(Int_t iColl) const { if (iColl == 0) return fPlane1Collection; return NULL; } // ---------------------------------------------------------------------------- // ----- Public method Print ---------------------------------------------- void CbmPlane::Print() const { } // ---------------------------------------------------------------------------- // ----- Public method Reset ---------------------------------------------- void CbmPlane::Reset() { fPlane1Collection->Clear(); fPlane2Collection->Clear(); fPlane3Collection->Clear(); fPosIndex = 0; } // ---------------------------------------------------------------------------- // guarda in FairRootManager::CopyClones // ----- Public method CopyClones ----------------------------------------- void CbmPlane::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset ) { } // ---------------------------------------------------------------------------- // ----- Public method ConstructGeometry ---------------------------------- void CbmPlane::ConstructGeometry() { FairGeoLoader* PlanegeoLoad = FairGeoLoader::Instance(); FairGeoInterface* PlanegeoFace = PlanegeoLoad->getGeoInterface(); CbmGeoPlane* PlaneGeo = new CbmGeoPlane(); PlaneGeo->setGeomFile(GetGeometryFileName()); PlanegeoFace->addGeoModule(PlaneGeo); Bool_t rc = PlanegeoFace->readSet(PlaneGeo); if (rc) PlaneGeo->create(PlanegeoLoad->getGeoBuilder()); TList* volList = PlaneGeo->getListOfVolumes(); // store geo parameter FairRun *fRun = FairRun::Instance(); FairRuntimeDb *rtdb= FairRun::Instance()->GetRuntimeDb(); CbmGeoPlanePar* par=(CbmGeoPlanePar*)(rtdb->getContainer("CbmGeoPlanePar")); TObjArray *fSensNodes = par->GetGeoSensitiveNodes(); TObjArray *fPassNodes = par->GetGeoPassiveNodes(); TListIter iter(volList); FairGeoNode* node = NULL; FairGeoVolume *aVol=NULL; while( (node = (FairGeoNode*)iter.Next()) ) { aVol = dynamic_cast ( node ); if ( node->isSensitive() ) { fSensNodes->AddLast( aVol ); }else{ fPassNodes->AddLast( aVol ); } } ProcessNodes ( volList ); par->setChanged(); par->setInputVersion(fRun->GetRunId(),1); } // ---------------------------------------------------------------------------- CbmPlanePoint* CbmPlane::AddHit1(Int_t trackID, Int_t detID, TVector3 posIn, TVector3 posOut, TVector3 momIn, TVector3 momOut, Double_t time, Double_t length, Double_t eLoss, Double_t charge, Double_t mass, Bool_t preflag, Bool_t postflag, TString nam,Int_t evt) { TClonesArray& clref = *fPlane1Collection; Int_t size = clref.GetEntriesFast(); return new(clref[size]) CbmPlanePoint(trackID, detID, posIn, posOut, momIn, momOut, time, length, eLoss, charge, mass, fpreflag, fpostflag, nam,evt); } CbmPlanePoint* CbmPlane::AddHit2(Int_t trackID, Int_t detID, TVector3 posIn, TVector3 posOut, TVector3 momIn, TVector3 momOut, Double_t time, Double_t length, Double_t eLoss, Double_t charge, Double_t mass, Bool_t preflag, Bool_t postflag, TString nam,Int_t evt) { TClonesArray& clref = *fPlane2Collection; Int_t size = clref.GetEntriesFast(); return new(clref[size]) CbmPlanePoint(trackID, detID, posIn, posOut, momIn, momOut, time, length, eLoss, charge, mass, fpreflag, fpostflag, nam,evt); } CbmPlanePoint* CbmPlane::AddHit3(Int_t trackID, Int_t detID, TVector3 posIn, TVector3 posOut, TVector3 momIn, TVector3 momOut, Double_t time, Double_t length, Double_t eLoss, Double_t charge, Double_t mass, Bool_t preflag, Bool_t postflag, TString nam,Int_t evt) { TClonesArray& clref = *fPlane3Collection; Int_t size = clref.GetEntriesFast(); return new(clref[size]) CbmPlanePoint(trackID, detID, posIn, posOut, momIn, momOut, time, length, eLoss, charge, mass, fpreflag, fpostflag, nam,evt); } ClassImp(CbmPlane)