///////////////////////////////////////////////////////////// // PndFtsHitProducerRealFull // // Class for digitalization for FTS // // authors: Pablo Genova - Pavia University // Lia Lavezzi - Pavia University // // modified for FTS by Isabella Garzia ///////////////////////////////////////////////////////////// #include "PndFtsHitProducerRealFull.h" #include "PndFtsHit.h" #include "PndFtsHitInfo.h" #include "PndFtsPoint.h" #include "PndFtsSingleStraw.h" #include "PndFtsMapCreator.h" #include "PndFtsTube.h" #include "PndFtsHitWriteoutBuffer.h" #include "FairRootManager.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "FairEventHeader.h" #include "TGeoManager.h" #include "TVector3.h" #include "TRandom.h" #include "TClonesArray.h" #include using std::cout; using std::endl; // ----- Default constructor ------------------------------------------- PndFtsHitProducerRealFull::PndFtsHitProducerRealFull() : FairTask("Real FTS Hit Producer",0), fPointArray(0), //fHitArray(0), fHitInfoArray(0), fFtsParameters(new PndGeoFtsPar()), fPersistence(kTRUE), fTimeOrderedDigi(kFALSE) { } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndFtsHitProducerRealFull::~PndFtsHitProducerRealFull() { } // ------------------------------------------------------------------------- // ----- Public method Init -------------------------------------------- InitStatus PndFtsHitProducerRealFull::Init() { // Get RootManager FairRootManager* ioman = FairRootManager::Instance(); if ( ! ioman ) { cout << "-E- PndFtsHitProducerRealFull::Init: " << "RootManager not instantiated!" << endl; return kFATAL; } // Get input array fPointArray = dynamic_cast (ioman->GetObject("FTSPoint")); if ( ! fPointArray ) { cout << "-W- PndFtsHitProducerRealFull::Init: " << "No FTSPoint array!" << endl; return kERROR; } // Create and register output array // fHitArray = new TClonesArray("PndFtsHit"); // ioman->Register("FTSHit","FTS",fHitArray, fPersistence); fDataBuffer = new PndFtsHitWriteoutBuffer("FTSHit", "FTS", fPersistence); fDataBuffer = (PndFtsHitWriteoutBuffer*)ioman->RegisterWriteoutBuffer("FTSHit", fDataBuffer); fDataBuffer->ActivateBuffering(fTimeOrderedDigi); // Create and register output array fHitInfoArray = new TClonesArray("PndFtsHitInfo"); ioman->Register("FTSHitInfo", "FTS", fHitInfoArray, kFALSE); // CHECK added PndFtsMapCreator *mapper = new PndFtsMapCreator(fFtsParameters); fTubeArray = mapper->FillTubeArray(); cout << "-I- PndFtsHitProducerRealFull: Intialization successfull" << endl; return kSUCCESS; } // ------------------------------------------------------------------------- // CHECK added void PndFtsHitProducerRealFull::SetParContainers() { FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb(); fFtsParameters = (PndGeoFtsPar*) rtdb->getContainer("PndGeoFtsPar"); } // ----- Public method Exec -------------------------------------------- void PndFtsHitProducerRealFull::Exec(Option_t*) { // Int_t evtn=0; // PndFtsPoint *ptemp =(PndFtsPoint*) fPointArray->At(0); // if(ptemp !=NULL) { // evtn=ptemp->GetEventID(); // if(evtn%50==0)cout << "Event Number "<Delete(); fHitInfoArray->Clear(); Int_t detID = 0; // detectorID TVector3 pos, dpos; // position and error vectors // Declare some variables PndFtsPoint* point = NULL; // Loop over FtsPoints Int_t nPoints = fPointArray->GetEntriesFast(); for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) { point = (PndFtsPoint*) fPointArray->At(iPoint); if (point == NULL) continue; detID = point->GetDetectorID(); // tubeID CHECK added Int_t skew=0; Int_t layerID = point->GetLayerID(); Int_t tubeID = point->GetTubeID(); Int_t chamberID = point->GetChamberID(); PndFtsTube *tube = (PndFtsTube*) fTubeArray->At(tubeID); //if skewed tube: skew==1 if(layerID>=3 && layerID<=6){skew=1;} //skewed tudes fts1 if(layerID>=11 && layerID<=14){skew=1;} //skewed tudes fts2 if(layerID>=19 && layerID<=22){skew=1;} //skewed tudes fts3 if(layerID>=27 && layerID<=30){skew=1;} //skewed tudes fts4 if(layerID>=35 && layerID<=38){skew=1;} //skewed tudes fts5 if(layerID>=43 && layerID<=46){skew=1;} //skewed tudes fts6 double InOut[6]; memset(InOut, 0, sizeof(InOut)); InOut[0] = point->GetXInLocal(); InOut[1] = point->GetYInLocal(); InOut[2] = point->GetZInLocal(); InOut[3] = point->GetXOutLocal(); InOut[4] = point->GetYOutLocal(); InOut[5] = point->GetZOutLocal(); // single straw tube simulation ----------------------- PndFtsSingleStraw fts; //setting the single straw tube simulation constants // 3 options currently available: // TConst(tube radius (cm), gas pressure (bar), Ar%, CO2%) // fts.TConst(0.4, 1, 0.9, 0.1); // fts.TConst(0.5, 1, 0.9, 0.1); fts.TConst(0.5, 2, 0.8, 0.2); // wire positioning fts.PutWireXYZ(0., 0., -75., 0., 0., 75.); // get particle momentum TVector3 momentum(point->GetPxOut(),point->GetPyOut(),point->GetPzOut()); // GeV/c Double_t GeV=1.; // position in cm (already in cm); momentum in GeV (already in GeV); mass in GeV (already in GeV) // drift time calculation Double_t pulset = fts.PartToTime(point->GetMass()/GeV, momentum.Mag()/GeV, InOut); // simulated radius (cm) double radius = fts.TimnsToDiscm(pulset); if(radius < 0.) radius = 0.; // CHECK // true radius (cm) //double true_rad = fts.TrueDist(InOut); //[R.K. 03/2017] unused variable? // dE calculation double depCharge = fts.PartToADC(); // dE/dx calculation //double dedx = -999; //[R.K. 01/2017] unused variable? // fts: detID, pos, dpos, index come from -------------- // fts (FairHit): Double_t closestDistanceError = 0.0150; // per adesso (stessa che in Ideal: // radialResolution = 0.0150) TVector3 position = tube->GetPosition(); // ---------------- // TVector3 posInLocal(point->GetXInLocal(), point->GetYInLocal(), point->GetZInLocal()); // TVector3 posOutLocal(point->GetXOutLocal(), point->GetYOutLocal(), point->GetZOutLocal()); // Double_t zpos = position.Z() + ((posOutLocal.Z() + posInLocal.Z()) / 2.); // Double_t zposError; // FoldZPosWithResolution(zpos, zposError, posInLocal, posOutLocal); // pos.SetXYZ(position.X(), position.Y(), zpos); // <--- stt2 // ---------------- pos.SetXYZ(position.X(), position.Y(), position.Z()); // <--- stt1 // dpos.SetXYZ(innerStrawDiameter / 2., innerStrawDiameter / 2., GetLongitudinalResolution(position.Z())); dpos.SetXYZ(0.5, 0.5, 3.); // per adesso (stessi che in Ideal: // innerStrawDiameter/2 = 0.5, // longitudinalResolution = 3.) // create hit AddHit(detID, tubeID, chamberID, layerID, skew, iPoint, pos, dpos, pulset, radius, closestDistanceError, depCharge, point->GetTime()); AddHitInfo(0, 0, point->GetTrackID(), iPoint, 0, kFALSE); }// Loop over MCPoints // Event summary cout << "-I- PndFtsHitProducerRealFull: " << nPoints << " FtsPoints, " << nPoints << " Hits created." << endl; } // ------------------------------------------------------------------------- void PndFtsHitProducerRealFull::FoldZPosWithResolution(Double_t &zpos, Double_t &zposError, TVector3 , TVector3 ) // localInPos localOutPos //[R.K.03/2017] unused variable(s) { //Double_t //zPosInStrawFrame = (localOutPos.Z() - localInPos.Z()) / 2.; //[R.K. 01/2017] unused variable? // zposError = gRandom->Gaus(0., GetLongitudinalResolution(zPosInStrawFrame)); zposError = gRandom->Gaus(0., 3.); // per adesso (stesso che in Ideal: // longitudinalResolution = 3.) zpos += zposError; } // ----- Private method AddHit -------------------------------------------- PndFtsHit* PndFtsHitProducerRealFull::AddHit(Int_t detID, Int_t tubeID, Int_t chamberID, Int_t layerID, Int_t skew, Int_t iPoint, TVector3& pos, TVector3& dpos, Double_t p, Double_t rsim, Double_t closestDistanceError, Double_t depcharge, Double_t timeOfFlight) { // see PndFtsHit for hit description Double_t EventTime = FairRootManager::Instance()->GetEventTime(); PndFtsHit *hitnew = new PndFtsHit(detID, tubeID, chamberID, layerID, skew, iPoint, pos, dpos, p+EventTime+timeOfFlight, rsim, closestDistanceError, depcharge); if (fTimeOrderedDigi){ hitnew->ResetLinks(); FairEventHeader* evtHeader = (FairEventHeader*)FairRootManager::Instance()->GetObject("EventHeader."); hitnew->AddLink(FairLink(evtHeader->GetInputFileId(), evtHeader->GetMCEntryNumber(), "FTSPoint", iPoint)); hitnew->AddLink(FairLink(-1, FairRootManager::Instance()->GetEntryNr(), "EventHeader.", -1)); } fDataBuffer->FillNewData(hitnew, p+EventTime+timeOfFlight, timeOfFlight+EventTime); return hitnew; } // ---- // ----- Private method AddHitInfo -------------------------------------------- PndFtsHitInfo* PndFtsHitProducerRealFull::AddHitInfo(Int_t fileNumber, Int_t eventNumber, Int_t trackID, Int_t pointID, Int_t nMerged, Bool_t isFake){ // see PndSttHitInfo for hit description TClonesArray& clref = *fHitInfoArray; Int_t size = clref.GetEntriesFast(); return new(clref[size]) PndFtsHitInfo(fileNumber, eventNumber, trackID, pointID, nMerged, isFake); } ClassImp(PndFtsHitProducerRealFull)