// ------------------------------------------------------------------------- // ----- PndDrcDigiTask source file ----- // ----- Created 15/04/13 by M.Patsyuk ----- // ----- HARPHOOL KUMAWAT h.kumawat@gsi.de ----- // ----- ----- // ------------------------------------------------------------------------- #include #include #include #include "stdio.h" #include "TGeoManager.h" #include "FairRootManager.h" #include "PndDrcDigiTask.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 "FairBaseParSet.h" #include "FairGeoTransform.h" #include "PndGeoDrcPar.h" using std::endl; using std::cout; // ----- Default constructor ------------------------------------------- PndDrcDigiTask::PndDrcDigiTask() :FairTask("PndDrcDigiTask"){ fGeo = new PndGeoDrc(); fGeoH = NULL; fPDPointArray = NULL; SetParameters(); Reset(); } // ----- Standard constructor with verbosity level ------------------------------------------- PndDrcDigiTask::PndDrcDigiTask(Int_t verbose) :FairTask("PndDrcDigiTask",verbose){ fVerbose = verbose; fGeo = new PndGeoDrc(); fGeoH = NULL; fPDPointArray = NULL; SetParameters(); Reset(); } void PndDrcDigiTask::SetParameters(){ fTimeOrderedDigi = kFALSE; fChargeSharing = kFALSE; fDeadTime = 5; fDetType=1; // Detector Type =1 fSigmat=0.1; // Time Resolution is 100 ps ############################ fTimeGranularity = 0.001;// [ns] = 98 ps granularity of the time signal 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(); fThreshold = 0.20; //threshold to detect charge shared hits } // ----- Destructor ---------------------------------------------------- PndDrcDigiTask::~PndDrcDigiTask(){ Reset(); if (fGeo) delete fGeo; } // ----- Initialization ----------------------------------------------- InitStatus PndDrcDigiTask::Init(){ cout << " ---------- INITIALIZATION ------------" << endl; nevents = 0; // Get RootManager FairRootManager* ioman = FairRootManager::Instance(); if ( ! ioman ) { cout << "-E- PndDrcDigiTask::Init: " << "RootManager not instantiated!" << endl; return kFATAL; } // PndGeoHandling if ( fGeoH == NULL ) fGeoH = PndGeoHandling::Instance(); fGeoH->SetParContainers(); fGeoH->SetVerbose(fVerbose); 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- PndDrcDigiTask::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; } // Create and register output buffer fDataBuffer = new PndDrcDigiWriteoutBuffer("DrcDigi","PndDrc", kTRUE); fDataBuffer = (PndDrcDigiWriteoutBuffer*)ioman->RegisterWriteoutBuffer("DrcDigi", fDataBuffer); fDataBuffer->SetVerbose(fVerbose); fDataBuffer->ActivateBuffering(fTimeOrderedDigi); cout << "-I- PndDrcDigiTask: Intialization successfull" << endl; return kSUCCESS; } // ----- Execution of Task --------------------------------------------- void PndDrcDigiTask::Exec(Option_t*){ // ion //[R.K.03/2017] unused variable(s) Reset(); if (fVerbose > 0) cout << "-I- PndDrcDigiTask: Event # " << nevents << " has " << fPDPointArray->GetEntriesFast() << " PD points" << endl; ProcessPhotonPoint(); nevents++; } //-------- Photon Detector Hit production with efficiency------------- void PndDrcDigiTask::ProcessPhotonPoint(){ for(Int_t k=0; k < fPDPointArray->GetEntriesFast(); k++) { fPpt = (PndDrcPDPoint*)fPDPointArray->At(k); // fMCtrk = (PndMCTrack*)fMCArray->At(fPpt->GetTrackID()); Int_t bpId=fPpt->GetBarPointID(); if(bpId<0){ std::cout<<" Error: bpId = "<At(bpId); // transform to local sensor system TVector3 PptPosition; fPpt->Position(PptPosition); TVector3 posL = fGeoH->MasterToLocalShortId(PptPosition,fPpt->GetDetectorID()); 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; fTime=fPpt->GetTime(); if(fSigmat>0) Smear(fTime,fSigmat); ActivatePixel(sensorId, k, 0); if(fChargeSharing){ Double_t distance = 999.; TVector3 corner, point; // left pixel distance = posLshifted.X() - TMath::Floor(posLshifted.X()/fPixelStep)*fPixelStep; if(exp(- distance/fPixelSigma) > gRandom->Uniform(0.,1.) && exp(- distance/fPixelSigma) > fThreshold){ if((Ncol-1) >= 0){ ActivatePixel(mcpId * 100 +Ncol-1+Nrow*fNpix,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(mcpId * 100 +Ncol+1+Nrow*fNpix, 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(mcpId * 100 +Ncol+(Nrow-1)*fNpix, 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(mcpId * 100 +Ncol+(Nrow+1)*fNpix, 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(mcpId * 100 +Ncol-1+(Nrow+1)*fNpix, 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(mcpId * 100 +Ncol-1+(Nrow-1)*fNpix, 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(mcpId * 100 +Ncol+1+(Nrow-1)*fNpix, 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(mcpId * 100 +Ncol+1+(Nrow+1)*fNpix, k, 1); } } } } } // ----- Private method ActivatePixel ------------------------------- void PndDrcDigiTask::ActivatePixel(Int_t sensorId, 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 PndDrcDigi* digi = new PndDrcDigi(k, fDetectorID, sensorId, 0., fTime, csflag, 0.); Double_t timeStamp = fTime + FairRootManager::Instance()->GetEventTime(); digi->SetPdgCode(fBarPoint->GetPdgCode()); digi->SetTimeStamp(timeStamp); digi->SetTimeAtBar(fBarPoint->GetTime()); if(fSigmat>0) digi->SetTimeStampError(fSigmat/TMath::Sqrt(fSigmat)); else digi->SetTimeStampError(0); FairEventHeader* evtHeader = (FairEventHeader*)FairRootManager::Instance()->GetObject("EventHeader."); digi->SetLink(fPpt->GetLink(0)); // MCTrack digi->AddLink(FairLink(evtHeader->GetInputFileId(), evtHeader->GetMCEntryNumber(), "DrcPDPoint", k)); fDataBuffer->FillNewData(digi, timeStamp, timeStamp + fDeadTime); // if ( fPixelMap.find(sensorId) == fPixelMap.end() ){ // // pixel not yet active, create new digi // fDataBuffer->FillNewData(digi, fTime, fTime + fDeadTime); // fPixelMap[sensorId] = fNDigis; // fNDigis++; // }else{ // // // if there was already hit in this pixel... // // PndDrcDigi* ddigi = dynamic_cast(fDrcDigiArray->At(fPixelMap[sensorId])); // // // ... check the time difference with another one in the same pixel // // if(fabs(ddigi->GetTimeStamp() - fTime) > fDeadTime){ // fDataBuffer->FillNewData(digi, fTime, fTime + fDeadTime); // fPixelMap[sensorId] = fNDigis; // fNDigis++; // // } // } } //-------------------------------------------------------------- TVector3 PndDrcDigiTask::GetSensorDimensions(Int_t sensorID){ gGeoManager->cd(fGeoH->GetPath(sensorID)); TGeoVolume* actVolume = gGeoManager->GetCurrentVolume(); TGeoBBox* actBox = (TGeoBBox*)(actVolume->GetShape()); TVector3 result(actBox->GetDX(),actBox->GetDY(),actBox->GetDZ()); return result; } //-------------Smear Time--------------------------------------- void PndDrcDigiTask::Smear(Double_t& time, Double_t sigt){ time += gRandom->Gaus(0,sigt); } // ----- Private method Reset ------------------------------ void PndDrcDigiTask::Reset() { fNDigis = 0; fPixelMap.clear(); } // ----- Finish Task --------------------------------------- void PndDrcDigiTask::Finish(){ cout << "-I- PndDrcDigiTask: Finish" << endl; } ClassImp(PndDrcDigiTask)