//---------------------------------------------------------------------- // File and Version Information: // $Id: // // Description: // Class PndEmcMultiWaveformToCalibratedDigi. Module to take the ADC waveforms and produces digi. // // Software developed for the BaBar Detector at the SLAC B-Factory. // Adapted for the PANDA experiment at GSI // // Author List: // Phil Strother Original Author // Dima Melnichuk - adaption for PANDA // Copyright Information: // Copyright (C) 1996 Imperial College // //---------------------------------------------------------------------- #include "PndEmcMultiWaveformToCalibratedDigi.h" #include "PndEmcMultiWaveform.h" #include "PndEmcWaveform.h" #include "PndEmcDigi.h" #include "PndEmcDigiPar.h" #include "PndEmcRecoPar.h" #include "PndEmcAsicPulseshape.h" #include "PndEmcPSAParabolic.h" #include "PndEmcPSAParabolicBaseline.h" #include "PndEmcPSAMatchedDigiFilter.h" #include "FairRootManager.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "TClonesArray.h" #include "TStopwatch.h" #include #include //#include using std::cout; using std::endl; using std::fstream; PndEmcMultiWaveformToCalibratedDigi::PndEmcMultiWaveformToCalibratedDigi(Int_t verbose, Bool_t storedigis) { fVerbose=verbose; fDigiPosMethod="depth";// "surface" or "depth" fEmcDigiRescaleFactor=1.08; fStoreDigis=storedigis; fCalibrationFileName = ""; //fPndEmcDigiPositionDepth=6.2; } //-------------- // Destructor -- //-------------- PndEmcMultiWaveformToCalibratedDigi::~PndEmcMultiWaveformToCalibratedDigi() { } InitStatus PndEmcMultiWaveformToCalibratedDigi::Init() { // Get RootManager FairRootManager* ioman = FairRootManager::Instance(); if ( ! ioman ) { cout << "-E- PndEmcMultiWaveformToCalibratedDigi::Init: " << "RootManager not instantiated!" << endl; return kFATAL; } // Get input array fWaveformArray = (TClonesArray*) ioman->GetObject("EmcMultiWaveform"); if ( ! fWaveformArray ) { cout << "-W- PndEmcMultiWaveformToCalibratedDigi::Init: " << "No PndEmcMultiWaveform array!" << endl; return kERROR; } // Create and register output array fDigiArray = new TClonesArray("PndEmcDigi"); ioman->Register("EmcDigi","Emc",fDigiArray,fStoreDigis); fSampleRate=fDigiPar->GetSampleRate(); fSampleRate_PMT=fDigiPar->GetSampleRate_PMT(); fASIC_Shaping_int_time=fDigiPar->GetASIC_Shaping_int_time(); //s fPMT_Shaping_int_time=fDigiPar->GetPMT_Shaping_int_time(); //s fPMT_Shaping_diff_time=fDigiPar->GetPMT_Shaping_diff_time(); //s fCrystal_time_constant=fDigiPar->GetCrystal_time_constant(); //s fShashlyk_time_constant=fDigiPar->GetShashlyk_time_constant(); //s fNumber_of_samples_in_waveform=fDigiPar->GetNumber_of_samples_in_waveform(); fNumber_of_samples_in_waveform_pmt=fDigiPar->GetNumber_of_samples_in_waveform_pmt(); fEnergyDigiThreshold=fDigiPar->GetEnergyDigiThreshold(); fEmcDigiPositionDepthPWO=fRecoPar->GetEmcDigiPositionDepthPWO(); fEmcDigiPositionDepthShashlyk=fRecoPar->GetEmcDigiPositionDepthShashlyk(); cout<<"fEmcDigiPositionDepthPWO: "< params; params.push_back(30); // width params.push_back(fSampleRate); // Sample rate fpsaAlgorithm = new PndEmcPSAMatchedDigiFilter(params,fPulseshape); } if(fpsaAlgorithm_pmt == NULL){ // Pulse shape analysis algorithm. // Simple parabolic fit. fpsaAlgorithm_pmt = new PndEmcPSAParabolic(); } // Determine normalisation constant for PndEmcMultiWaveform PndEmcWaveform *tmpwaveform=new PndEmcWaveform(0,101010001, fNumber_of_samples_in_waveform); PndEmcWaveform *tmpwaveform2=new PndEmcWaveform(0,101010001, fNumber_of_samples_in_waveform_pmt); PndEmcHit *gevHit=new PndEmcHit(); gevHit->SetEnergy(1.0); gevHit->SetTime(0.); tmpwaveform->UpdateWaveform(gevHit, 0, false, 1., 0., fSampleRate, fPulseshape); tmpwaveform2->UpdateWaveform(gevHit, 0, false, 1., 0., fSampleRate_PMT, fPulseshape_pmt); Double_t tmpPeakPosition; Double_t tmpPeakPosition2; fpsaAlgorithm->Process(tmpwaveform,fWfNormalisation,tmpPeakPosition); fpsaAlgorithm_pmt->Process(tmpwaveform2,fWfNormalisation_pmt,tmpPeakPosition2); fWfNormalisation_proto192 = 1; if(fCalibrationFileName !="") ReadCalibrationFile(); cout << "-I- PndEmcMultiWaveformToCalibratedDigi: Read "<< fCalibrationMap.size() << " Calibration Entries" << endl; cout << "-I- PndEmcMultiWaveformToCalibratedDigi: Intialization successfull" << endl; return kSUCCESS; } void PndEmcMultiWaveformToCalibratedDigi::Exec(Option_t*) { TStopwatch timer; if (fVerbose>2){ timer.Start(); } // Reset output array if ( ! fDigiArray ) Fatal("Exec", "No Digi Array"); fDigiArray->Delete(); Double_t peakPosition; Double_t energy; Double_t digi_time; Int_t i_digi=0; //index of digi in TClonesArray Int_t hitIndex; Int_t detId; Int_t trackId; Int_t module; Int_t nWaveforms = fWaveformArray->GetEntriesFast(); Int_t nSignal; Bool_t highgain; PndEmcAbsPSA *thePSA; Double_t theEnergyNorm; Double_t theSampleRate; Int_t nHits[4]; Double_t hittimes[4][32]; Double_t hitenergies[4][32]; std::map::iterator it; //cout<<"PndEmcMultiWaveformToCalibratedDigi: "<At(iWaveform); hitIndex=theWaveform->GetHitIndex(); detId=theWaveform->GetDetectorId(); trackId=theWaveform->GetTrackId(); module=theWaveform->GetModule(); nSignal=theWaveform->GetNumberOfWaveforms(); if(module == 5){ thePSA = fpsaAlgorithm_pmt; theEnergyNorm = fWfNormalisation_pmt; theSampleRate = fSampleRate_PMT; } else { thePSA = fpsaAlgorithm; theEnergyNorm = fWfNormalisation_proto192;; theSampleRate = fSampleRate; } for(Int_t iSignal=0; iSignal < nSignal; iSignal++){ theWaveform->SetActiveWaveform(iSignal); nHits[iSignal] = thePSA->Process(theWaveform); // std::cout << "Signal " << iSignal << ": " << nHits[iSignal] << " Hits" << std::endl; for(Int_t iHit = 0; iHit < nHits[iSignal]; iHit ++){ thePSA->GetHit(iHit,hitenergies[iSignal][iHit],hittimes[iSignal][iHit]); // std::cout << "Signal " << iSignal << ", Hit " << iHit << ": Energy " << hitenergies[iSignal][iHit] < 1 && energy > 1500.){ //todo make highgain limit a config option //look for hit in lowgain // std::cout << "Hit " << iHit << " is too high (" << energy << ") checking lowgain." << std::endl; // std::cout << "highgain time: " << hittimes[0][iHit] << std::endl; for(Int_t iHit1 = 0; iHit1second << " for hitIndex " << hitIndex << std::endl; energy*=it->second; }else{ // if(hitIndex == 45) std::cout << "no calibration value for hitIndex " << hitIndex << std::endl; } if(!highgain){ // if(hitIndex == 45) std::cout << "highgain" << std::endl; it=fGainMap.find(detId); if(it!=fGainMap.end()){ if(hitIndex == 45) std::cout << "found highgain value of " << it->second << " for hitIndex " << hitIndex << std::endl; energy*=it->second; } } // if(energy>1000){ // if(hitIndex == 45) std::cout << "energy: " << energy << " threshold: "<< fEnergyDigiThreshold << endl; // if(hitIndex == 45) std::cout << "creating digi for detid: " << detId << "hitIndex: " << hitIndex <fEnergyDigiThreshold) { // if(hitIndex == 45) std::cout << "energy: " << energy << endl; PndEmcDigi* myDigi = new((*fDigiArray)[i_digi]) PndEmcDigi(trackId,detId, energy, digi_time, hitIndex); myDigi->AddLink(FairLink("EmcMultiWaveform", iWaveform)); i_digi++; } } } if (fVerbose>2){ timer.Stop(); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); cout << "PndEmcMultiWaveformToCalibratedDigi, Real time " << rtime << " s, CPU time " << ctime << " s" << endl; } } void PndEmcMultiWaveformToCalibratedDigi::SetParContainers() { // Get run and runtime database FairRun* run = FairRun::Instance(); if ( ! run ) Fatal("SetParContainers", "No analysis run"); FairRuntimeDb* db = run->GetRuntimeDb(); if ( ! db ) Fatal("SetParContainers", "No runtime database"); // Get Emc digitisation parameter container fDigiPar = (PndEmcDigiPar*) db->getContainer("PndEmcDigiPar"); // Get Emc reconstruction parameter container fRecoPar = (PndEmcRecoPar*) db->getContainer("PndEmcRecoPar"); } void PndEmcMultiWaveformToCalibratedDigi::SetStorageOfData(Bool_t val) { fStoreDigis = val; return; } void PndEmcMultiWaveformToCalibratedDigi::ReadCalibrationFile(){ std::ifstream in(fCalibrationFileName,std::ifstream::in); // in.open(fCalibrationFileName); if (!in.good()) { std::cerr << "Cannot open calibration file!" << endl; return; } char buf[255]; while (in.getline(buf, 255)) { TString tmp = buf; TObjArray *tokens = tmp.Tokenize(" \t;"); if(tokens->GetEntries()<2){ continue; } tmp=tokens->UncheckedAt(0)->GetName(); TString tmp2 = tokens->UncheckedAt(1)->GetName(); if(tmp.IsDigit() && tmp2.IsFloat()){ fCalibrationMap.insert(std::pair(tmp.Atoi(),tmp2.Atof())); } if(tokens->GetEntries()<3){ continue; } tmp2 = tokens->UncheckedAt(2)->GetName(); if(tmp.IsDigit() && tmp2.IsFloat()){ fGainMap.insert(std::pair(tmp.Atoi(),tmp2.Atof())); } delete tokens; } if(fVerbose >1){ std::cout << "calibration values:" << std::endl; std::map::iterator it; std::map::iterator it2; Int_t detid; for ( it=fCalibrationMap.begin() ; it != fCalibrationMap.end(); it++ ){ detid = (*it).first; std::cout << detid << "\t" << (*it).second << "\t"; it2=fGainMap.find(detid); if(it2!=fGainMap.end()){ std::cout << (*it2).second; } std::cout << endl; } } } ClassImp(PndEmcMultiWaveformToCalibratedDigi)