// ------------------------------------------------------------------------- // ----- PndDrcTimeDigiTask source file ----- // ----- HARPHOOL KUMAWAT h.kumawat@gsi.de ----- // ----- ----- // ------------------------------------------------------------------------- #include #include #include #include "stdio.h" #include "TGeoManager.h" #include "FairRootManager.h" #include "PndDrcDigiWriteoutBuffer.h" #include "PndDrcTimeDigiTask.h" #include "PndDrcDigi.h" #include "PndDrcPDPoint.h" #include "PndDrcBarPoint.h" #include "PndMCTrack.h" #include "FairHit.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "FairGeoVector.h" #include "FairGeoNode.h" #include "FairGeoVolume.h" #include "TVector3.h" #include "TRandom.h" #include "TString.h" #include "TGeoBBox.h" #include "TCanvas.h" #include "FairBaseParSet.h" #include "FairGeoTransform.h" #include "PndGeoDrcPar.h" using std::endl; using std::cout; // ----- Default constructor ------------------------------------------- PndDrcTimeDigiTask::PndDrcTimeDigiTask() :FairTask("PndDrcTimeDigiTask") { fGeo = new PndGeoDrc(); fGeoH = NULL; SetParameters(); fisDetEff= kTRUE; fisPixel= kTRUE; fChargeSharing=kTRUE; fTimeSmearing= kTRUE; fPDPointArray = NULL; fTimeOrderedDigi =kFALSE;// kTRUE; fDrcDigiArray = NULL; fDrcTimeDigiArray = NULL; Reset(); } // ------------------------------------------------------------------------- // ----- Standard constructor with verbosity level ------------------------------------------- PndDrcTimeDigiTask::PndDrcTimeDigiTask(Int_t verbose) :FairTask("PndDrcTimeDigiTask",verbose) { fVerbose = verbose; fDetType= 1; fisDetEff= kTRUE; fisPixel= kTRUE; fGeo = new PndGeoDrc(); fGeoH = NULL; SetParameters(); fPDPointArray = NULL; fTimeOrderedDigi=kFALSE;//kTRUE; fDrcDigiArray = NULL; fDrcTimeDigiArray = NULL; Reset(); } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndDrcTimeDigiTask::~PndDrcTimeDigiTask() { if (fGeo) delete fGeo; if ( fDrcDigiArray ) { fDrcDigiArray->Delete(); delete fDrcDigiArray; } Reset(); } // ----- Initialization ----------------------------------------------- // ------------------------------------------------------------------------- InitStatus PndDrcTimeDigiTask::Init() { cout << " ---------- INITIALIZATION ------------" << endl; nevents = 0; // Get RootManager FairRootManager* ioman = FairRootManager::Instance(); if ( ! ioman ) { cout << "-E- PndDrcTimeDigiTask::Init: " << "RootManager not instantiated!" << endl; return kFATAL; } // PndGeoHandling if ( fGeoH == NULL ) fGeoH = PndGeoHandling::Instance(); fGeoH->SetParContainers(); if(fVerbose>1) Info("SetParContainers","done."); // Get input array fMCArray = (TClonesArray*) ioman->GetObject("MCTrack"); if ( ! fMCArray ) { cout << "-W- PndDrcRecoLookupMap::Init: " << "No MCTrack array!" << endl; return kERROR; } // Get input array fBarPointArray = (TClonesArray*) ioman->GetObject("DrcBarPoint"); if ( ! fBarPointArray ) { cout << "-W- PndDrcTimeDigiTask::Init: " << "No DrcBarPoint array!" << endl; return kERROR; } // Get Photon point array fPDPointArray = (TClonesArray*) ioman->GetObject("DrcPDPoint"); if ( ! fPDPointArray ) { cout << "-W- PndDrcAna::Init: " << "No DrcPDPoint array!" << endl; return kERROR; } fDrcDigiArray = new TClonesArray("PndDrcDigi"); ioman->Register("DrcDigiNormal","Drc",fDrcDigiArray, kFALSE); // Create and register output array fDrcTimeDigiArray = new PndDrcDigiWriteoutBuffer("DrcDigi","PndDrc", kTRUE); fDrcTimeDigiArray = (PndDrcDigiWriteoutBuffer*)ioman->RegisterWriteoutBuffer("DrcDigi", fDrcTimeDigiArray); fDrcTimeDigiArray->ActivateBuffering(fTimeOrderedDigi); cout << "-I- PndDrcTimeDigiTask: Intialization successfull" << endl; return kSUCCESS; } // ----- Execution of Task --------------------------------------------- void PndDrcTimeDigiTask::Exec(Option_t*) // ion //[R.K.03/2017] unused variable(s) { fTimeMap.clear(); Reset(); fDrcDigiArray->Delete(); nevents++; fPileup=0; //Int_t PDHIT = fPDPointArray->GetEntriesFast(); //[R.K. 01/2017] unused variable? // if (fVerbose > 0) if(nevents%10==0) // cout<<"event time "<GetEventTime()<SetVerbose(fVerbose); if (fVerbose > 1) { cout <<" Number of Photon MC Points in Photon Detector Plane : "<GetEntries()<GetEntriesFast() <GetEntriesFast(); k++) { fPpt = (PndDrcPDPoint*)fPDPointArray->At(k); Int_t trID= fPpt->GetTrackID(); if(trID<0) continue; fMCtrk = (PndMCTrack*)fMCArray->At(trID); Int_t trMrID= fMCtrk->GetMotherID(); fMrIdPho = trMrID; fTimeAtBar = fMCtrk->GetStartTime(); PndMCTrack* fMCtrk1 = (PndMCTrack*)fMCArray->At(trMrID); fPdg = fMCtrk1->GetPdgCode(); fMrId = fMCtrk1->GetMotherID(); //Int_t NbouncesX, NbouncesY; //[R.K. 01/2017] unused variable? //Double_t angleX, angleY; //[R.K. 01/2017] unused variable? Int_t BarId=-1; fTrackId= trID; fTrackMom.SetXYZ(fMCtrk1->GetMomentum().X(),fMCtrk1->GetMomentum().Y(),fMCtrk1->GetMomentum().Z()); if(fPpt->GetBarPointID()!=-1){ fBarPoint= (PndDrcBarPoint*)fBarPointArray->At(fPpt->GetBarPointID()); BarId = fBarPoint->GetDetectorID(); fBarId = fBarPoint->GetBarId(); fBoxId = fBarPoint->GetBoxId(); // start vertex of the photon TVector3 StartVertex = fMCtrk->GetStartVertex(); fTrackIniVertex = StartVertex; // initial direction of the photon in the bar coord system TVector3 PphoInitBar = fGeoH->MasterToLocalShortId(fMCtrk->GetMomentum(), BarId)- fGeoH->MasterToLocalShortId(TVector3(0.,0.,0.),BarId); // vector } //Double_t PPx= fPpt->GetPx(); //[R.K. 01/2017] unused variable? //Double_t PPy= fPpt->GetPy(); //[R.K. 01/2017] unused variable? //Double_t PPz= fPpt->GetPz(); //[R.K. 01/2017] unused variable? //Double_t etot = sqrt(PPx*PPx + PPy*PPy +PPz*PPz);// in GeV //[R.K. 01/2017] unused variable? //Double_t lambda=197.0*2.0*fpi/(etot*1.0E9);//wavelength of photon in nm //[R.K. 01/2017] unused variable? //#################################### // transport efficiency if(BarId != -1){ fDetection=1; } if(fDetection==1){ // transform to local sensor system... (mc point has the ID not the path to the volume) TVector3 PptPosition; fPpt->Position(PptPosition); TVector3 posL = fGeoH->MasterToLocalShortId(PptPosition,fPpt->GetDetectorID()); // point TVector3 sensorDim = GetSensorDimensions(fPpt->GetDetectorID()); //usually sensors have origin in the middle, let's move it to the left lower corner: TVector3 posLshifted = posL + sensorDim; // calculate the number of the fired pixel Int_t Ncol = (Int_t)TMath::Floor(posLshifted.X()/fPixelStep); Int_t Nrow = (Int_t)TMath::Floor(posLshifted.Y()/fPixelStep); Int_t NpixelLocal = Ncol + Nrow * fNpix; fDetectorID = fPpt->GetDetectorID() * 100 + NpixelLocal;; Int_t mcpId = fPpt->GetMcpId(); Int_t sensorId = mcpId * 100 + NpixelLocal; // time is smeared and digitized = has granularity fTime=(Int_t)(fPpt->GetTime()/fTimeGranularity)*fTimeGranularity; fTimeOfFlight = fTime; fEventTime = FairRootManager::Instance()->GetEventTime(); if(fTimeSmearing == kTRUE) Smear(fTime,fSigmat); ActivatePixel(fDetectorID, sensorId, fTime, k, 0); if(fChargeSharing == kTRUE){ // find fired pixels: Double_t distance = 999.; TVector3 corner; TVector3 point; // left pixel distance = posLshifted.X() - TMath::Floor(posLshifted.X()/fPixelStep)*fPixelStep; //[cm] if(exp(- distance/fPixelSigma) > gRandom->Uniform(0.,1.) && exp(- distance/fPixelSigma) > fThreshold){ if((Ncol-1) >= 0){ ActivatePixel(fDetectorID, mcpId * 100 +Ncol-1+Nrow*fNpix, fTime, k, 1); } } // right pixel distance = TMath::Ceil(posLshifted.X()/fPixelStep)*fPixelStep - posLshifted.X(); if(exp(- distance/fPixelSigma) > gRandom->Uniform(0.,1.) && exp(- distance/fPixelSigma) > fThreshold){ if(Ncol+1 < fNpix){ ActivatePixel(fDetectorID, mcpId * 100 +Ncol+1+Nrow*fNpix, fTime, k, 1); } } // lower pixel distance = posLshifted.Y() - TMath::Floor(posLshifted.Y()/fPixelStep)*fPixelStep; if(exp(- distance/fPixelSigma) > gRandom->Uniform(0.,1.) && exp(- distance/fPixelSigma) > fThreshold){ if(Nrow-1 >= 0){ ActivatePixel(fDetectorID, mcpId * 100 +Ncol+(Nrow-1)*fNpix, fTime, k, 1); } } // upper pixel distance = TMath::Ceil(posLshifted.Y()/fPixelStep)*fPixelStep - posLshifted.Y(); if(exp(- distance/fPixelSigma) > gRandom->Uniform(0.,1.) && exp(- distance/fPixelSigma) > fThreshold){ if(Nrow+1 < fNpix){ ActivatePixel(fDetectorID, mcpId * 100 +Ncol+(Nrow+1)*fNpix, fTime, k, 1); } } point.SetXYZ(posLshifted.X(),posLshifted.Y(),0.); // upper left pixel corner.SetXYZ(TMath::Floor(posLshifted.X()/fPixelStep)*fPixelStep, TMath::Ceil(posLshifted.Y()/fPixelStep)*fPixelStep,0.); distance = (point-corner).Mag(); if(exp(-distance/fPixelSigma) > gRandom->Uniform(0.,1.) && exp(- distance/fPixelSigma) > fThreshold){ if(Ncol-1 >= 0 && Nrow+1 < fNpix){ ActivatePixel(fDetectorID, mcpId * 100 +Ncol-1+(Nrow+1)*fNpix, fTime, k, 1); } } // bottom left pixel corner.SetXYZ(TMath::Floor(posLshifted.X()/fPixelStep)*fPixelStep, TMath::Floor(posLshifted.Y()/fPixelStep)*fPixelStep,0.); distance = (point-corner).Mag(); if(exp(-distance/fPixelSigma) > gRandom->Uniform(0.,1.) && exp(- distance/fPixelSigma) > fThreshold){ if(Ncol-1 >= 0 && Nrow-1 >= 0){ ActivatePixel(fDetectorID, mcpId * 100 +Ncol-1+(Nrow-1)*fNpix, fTime, k, 1); } } // bottom right pixel corner.SetXYZ(TMath::Ceil(posLshifted.X()/fPixelStep)*fPixelStep, TMath::Floor(posLshifted.Y()/fPixelStep)*fPixelStep,0.); distance = (point-corner).Mag(); if(exp(-distance/fPixelSigma) > gRandom->Uniform(0.,1.) && exp(- distance/fPixelSigma) > fThreshold){ if(Ncol+1 < fNpix && Nrow-1 >= 0){ ActivatePixel(fDetectorID, mcpId * 100 +Ncol+1+(Nrow-1)*fNpix, fTime, k, 1); } } // upper right pixel corner.SetXYZ(TMath::Ceil(posLshifted.X()/fPixelStep)*fPixelStep, TMath::Ceil(posLshifted.Y()/fPixelStep)*fPixelStep,0.); distance = (point-corner).Mag(); if(exp(-distance/fPixelSigma) > gRandom->Uniform(0.,1.) && exp(- distance/fPixelSigma) > fThreshold){ if(Ncol+1 < fNpix && Nrow+1 < fNpix){ ActivatePixel(fDetectorID, mcpId * 100 +Ncol+1+(Nrow+1)*fNpix, fTime, k, 1); } } } // if charge sharing true } // if detection } } // ----- Private method ActivatePixel ------------------------------- void PndDrcTimeDigiTask::ActivatePixel(Int_t detectorId, Int_t sensorId, Double_t signalTime, Int_t k, Int_t csflag) { // in case when the same pixel was fired by two different photons this function takes care of which hits from that pixel to write if ( fPixelMap.find(sensorId) == fPixelMap.end() ) { // pixel not yet active, create new digi PndDrcDigi* Digi = new PndDrcDigi(k, detectorId, sensorId, 0., fTimeOfFlight, csflag, 0.); Digi->SetTimeStamp(signalTime); Digi->SetTimeStampError(fSigmat/TMath::Sqrt(fSigmat)); Digi->SetBarID(fBarId); Digi->SetBoxID(fBoxId); Digi->SetTrackID(fTrackId); Digi->SetTrackIniVertex(fTrackIniVertex); Digi->SetMotherID(fMrId); Digi->SetMotherIDPho(fMrIdPho); Digi->SetPdgCode(fPdg); Digi->SetTrackMom(fTrackMom); Digi->SetTimeAtBar(fTimeAtBar); Digi->SetEvtTim(signalTime); Digi->SetEventTim(fEventTime); Digi->SetEventNo(nevents); fDrcTimeDigiArray->FillNewData(Digi, signalTime, signalTime + fActiveTime); fPixelMap[sensorId] = fNDigis; fTimeMap[sensorId] = signalTime; fNDigis++; Digi->SetLink(FairLink("PndDrcPDPoint", k)); } else { //Int_t iDigi = fPixelMap[sensorId]; //[R.K. 01/2017] unused variable? if(fabs(fTimeMap.find(sensorId)->second-signalTime)<5.0){ fPileup++; } { PndDrcDigi* digi = new PndDrcDigi(k, detectorId, sensorId, 0., fTimeOfFlight, csflag, 0.); digi->SetTimeStamp(signalTime); digi->SetTimeStampError(fSigmat/TMath::Sqrt(fSigmat)); digi->SetBarID(fBarId); digi->SetBoxID(fBoxId); digi->SetTrackID(fTrackId); digi->SetTrackIniVertex(fTrackIniVertex); digi->SetMotherID(fMrId); digi->SetMotherIDPho(fMrIdPho); digi->SetPdgCode(fPdg); digi->SetTrackMom(fTrackMom); digi->SetTimeAtBar(fTimeAtBar); digi->SetEvtTim(signalTime); digi->SetEventTim(fEventTime); digi->SetEventNo(nevents); fDrcTimeDigiArray->FillNewData(digi, signalTime, signalTime + fActiveTime); fPixelMap[sensorId] = fNDigis; fTimeMap[sensorId] = signalTime; fNDigis++; fPileup=0.0; } } } //---------------------------------------------------------------------------------------------- Double_t PndDrcTimeDigiTask::FindPhiRot(Double_t xx, Double_t yy){ // returns [degrees] TVector3 hit; hit.SetXYZ(xx,yy,0.); Double_t startPhi = hit.Phi()/fpi*180.; // [degrees] if(startPhi < 0.){startPhi = 360. + hit.Phi()*180./fpi;} Double_t PhiRot = 0.; //[degrees] if(startPhi >= 0. && startPhi < 90.){ PhiRot = TMath::Floor(startPhi/fdphi) *fdphi + fdphi/2.; } if(startPhi >= 90. && startPhi < 270.){ PhiRot = 90. + fpipehAngle + TMath::Floor((startPhi-90.-fpipehAngle)/fdphi) *fdphi + fdphi/2.; } if(startPhi >= 270. && startPhi < 360.){ PhiRot = 270. + fpipehAngle + TMath::Floor((startPhi-270.-fpipehAngle)/fdphi) *fdphi + fdphi/2.; } //cout<<"-I- FindPhiRot: PhiRot = "<= 0. && x1 + xEn <= a){xK = x1 + xEn;} if(x0 >= 0. && x1 + xEn > a){xK = 2*a - x1 - xEn; n = 1. + n;} if(x0 < 0. && x1 + xEn >= 0.){xK = a - (x1 + xEn); n = -1. -n;} if(x0 < 0. && x1 + xEn < 0.) {xK = a + x1 + xEn; n = -n;} if(print){std::cout<<"xK = "<< xK<<", n = "<= 0. && x1 + xEn <= a){xK = a - (x1 + xEn);} if(x0 >= 0. && x1 + xEn > a){xK = x1 + xEn - a; n = 1. + n;} if(x0 < 0. && x1 + xEn >= 0.){xK = x1 + xEn; n = -1. -n;} if(x0 < 0. && x1 + xEn < 0.) {xK = - (x1 + xEn); n = -n;} if(print){std::cout<<"xK = "<< xK<<", n = "<cd(fGeoH->GetPath(sensorID)); TGeoVolume* actVolume = gGeoManager->GetCurrentVolume(); TGeoBBox* actBox = (TGeoBBox*)(actVolume->GetShape()); TVector3 result(actBox->GetDX(),actBox->GetDY(),actBox->GetDZ()); return result; } //-------------Set Parameter------------------------------------ void PndDrcTimeDigiTask::SetParameters(){ fDetType=1; // Detector Type =1 nRefrac=fGeo->nEV();//1.467; //Refractive index of SOB // fSigmat=0.1; //Time Resolution is 100 ps ############################ fCollectionEff=0.65;//Collection Efficiency fPackingFraction=0.80;//Packing Efficiency fRoughness = 0.001; // 10 A fTimeGranularity = 0.001;// [ns] = 98 ps granularity of the time signal // basic DIRC parameters: fpi = TMath::Pi(); fzup = fGeo->barBoxZUp(); fzdown = fGeo->barBoxZDown(); fradius = fGeo->radius(); //radius in the middle of the bar = 50.cm fhthick = fGeo->barHalfThick(); //half thickness of the bars=1.7/2 cm fpipehAngle = fGeo->PipehAngle(); fbbGap = fGeo->BBoxGap(); fbbnum = fGeo->BBoxNum(); fbarnum = fGeo->barNum(); fphi0 = (180.-2.*fpipehAngle)/fbbnum + fpipehAngle; fdphi = (180.-2.*fpipehAngle)/fbbnum*2.; flside = fGeo->Lside(); fbarwidth = fGeo->BarWidth(); fPixelSize = fGeo->PixelSize(); fMcpActiveArea = fGeo->McpActiveArea(); fNpix = fGeo->Npixels(); fPixelGap = (fMcpActiveArea - (Double_t)fNpix*fPixelSize) / ((Double_t)fNpix - 1.); fPixelStep = fMcpActiveArea/fNpix; //fPixelSize + 0.5*fPixelGap; fPixelSigma = fGeo->SigmaCharge(); fDeadTime = fGeo->DeadTime(); fThreshold = 0.20; //threshold to detect charge shared hits } //-------------Smear Time------------------------------------ void PndDrcTimeDigiTask::Smear(Double_t& time, Double_t sigt){ Double_t EvtTime = FairRootManager::Instance()->GetEventTime(); time += EvtTime + gRandom->Gaus(0,sigt); } // ----- Private method Reset ------------------------------------------ void PndDrcTimeDigiTask::Reset() { fNDigis = 0; fPixelMap.clear(); if ( fDrcDigiArray ) fDrcDigiArray->Clear(); } // ------------------------------------------------------------------------- // ----- Finish Task --------------------------------------------------- void PndDrcTimeDigiTask::Finish() { if ( fDrcDigiArray ) fDrcDigiArray->Clear(); cout << "-I- PndDrcTimeDigiTask: Finish" << endl; } // ------------------------------------------------------------------------- ClassImp(PndDrcTimeDigiTask)