FairRoot/PandaRoot
PndEmcHitsToWaveform.cxx
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 // File and Version Information:
3 // $Id: //
4 // Description:
5 // Class PndEmcHitsToWaveform. Module to take the hit list for the
6 // calorimeter and make ADC waveforms from them.
7 //
8 // Software developed for the BaBar Detector at the SLAC B-Factory.
9 // Adapted for the PANDA experiment at GSI
10 //
11 // Author List:
12 // Phil Strother Original author
13 // Dima Melnichuk - adaption for PANDA
14 // Copyright Information:
15 // Copyright (C) 1996 Imperial College
16 //
17 //----------------------------------------------------------------------
18 
19 #include "PndEmcHitsToWaveform.h"
20 #include "FairEventHeader.h"
21 #include "PndEmcHit.h"
22 #include "PndEmcWaveform.h"
23 #include "PndEmcAsicPulseshape.h"
24 #include "PndEmcMapper.h"
25 #include "PndEmcStructure.h"
26 #include "PndEmcDigiPar.h"
27 #include "PndEmcGeoPar.h"
28 #include "PndEmcDataTypes.h"
29 #include "TRandom.h"
30 #include "FairRootManager.h"
31 #include "FairRunAna.h"
32 #include "FairRun.h"
33 #include "FairRuntimeDb.h"
34 #include "TGraph.h"
35 #include "TFile.h"
36 #include "TStopwatch.h"
37 #include "TROOT.h"
38 #include "TClonesArray.h"
39 
40 #include <iostream>
41 #include <cassert>
42 #include "PndEmcWaveformWriteoutBuffer.h"
43 //#include <map>
44 //#include <string>
45 
46 //typedef std::pair<Int_t, Double_t> Element;
47 //std::map<Int_t, std::vector<Element> > WaveformMap;
48 
49 using std::cout;
50 using std::endl;
51 using std::fstream;
52 //#define TIMEBASEDSIM
53 
54 
55 PndEmcHitsToWaveform::PndEmcHitsToWaveform(Int_t verbose, Bool_t storewaves):
56  fHitArray(0), fWaveformArray(0), fDataBuffer(0), fTimeOrderedWaveform(kFALSE), fOneBitResolution(0), fOneBitResolutionBW(0), fOneBitResolutionPMT(0), fNBits(0), fDetectedPhotonsPerMeV(0), fDetectedPhotonsPerMeV_PMT(0), fNPhotoElectronsPerMeVAPDBarrel(0), fNPhotoElectronsPerMeVAPDBWD(0), fNPhotoElectronsPerMeVVPT(0), fNPhotoElectronsPerMeVPMT(0), fSensitiveAreaAPD(0), fSensitiveAreaVPT(0), fQuantumEfficiencyAPD(0), fQuantumEfficiencyVPT(0), fQuantumEfficiencyPMT(0), fExcessNoiseFactorAPD(0), fExcessNoiseFactorVPT(0), fExcessNoiseFactorPMT(0), fIncoherent_elec_noise_width_GeV_APD(0), fIncoherent_elec_noise_width_GeV_VPT(0), fEnergyRange(0), fEnergyRangeBW(0), fFirstSamplePhase(0), fNumber_of_samples_in_waveform(0), fNumber_of_samples_in_waveform_pmt(0), fASIC_Shaping_int_time(0), fPMT_Shaping_int_time(0), fPMT_Shaping_diff_time(0), fCrystal_time_constant(0), fShashlyk_time_constant(0), fShashlykSamplingFactor(0), fSampleRate(0), fSampleRate_PMT(0), fUse_shaped_noise(0), fUse_photon_statistic(0), fNoiseAllChannels(0), fMapVersion(0), fFirstADCBinTime(0), fGevPeakAnalogue(0), fGevPeakAnalogue_PMT(0), fDigiPar(new PndEmcDigiPar()), fGeoPar(new PndEmcGeoPar()), fVerbose(verbose)
57 {
58  HowManyHit = 0;
59  SetPersistency(storewaves);
60 }
61 
62 //--------------
63 // Destructor --
64 //--------------
65 
66 PndEmcHitsToWaveform::~PndEmcHitsToWaveform()
67 {
68 // delete fw_endcap_noise;
69  delete pulseshape1;
70  delete pulseshape2;
71  delete pulseshape3;
72 }
73 
74 
86 InitStatus PndEmcHitsToWaveform::Init()
87 {
88  // Get RootManager
89  FairRootManager* ioman = FairRootManager::Instance();
90  if ( ! ioman )
91  {
92  cout << "-E- PndEmcHitsToWaveform::Init: "
93  << "RootManager not instantiated!" << endl;
94  return kFATAL;
95  }
96 
97  // Get input array
98  fHitArray = (TClonesArray*) ioman->GetObject("EmcHit");
99  if ( ! fHitArray ) {
100  cout << "-W- PndEmcHitsToWaveform::Init: "
101  << "No EmcHit array!" << endl;
102  return kERROR;
103  }
104  if(fTimeOrderedWaveform){
105  // Create and register output buffer
106  fDataBuffer = new PndEmcWaveformWriteoutBuffer("EmcWaveform", "Emc", GetPersistency());
107  fDataBuffer->ActivateBuffering(fTimeOrderedWaveform);
108  fDataBuffer->SetVerbose(fVerbose);
109  ioman->RegisterWriteoutBuffer("EmcWaveform", fDataBuffer);
110  fDataBuffer->SaveToTree(kTRUE);
111  }else{
112  // Create and register output array
113  fWaveformArray = new TClonesArray("PndEmcWaveform");
114  ioman->Register("EmcWaveform","Emc",fWaveformArray,GetPersistency());
115  }
116 
117  cout << "-I- PndEmcHitsToWaveform: Intialization successfull" << endl;
118 
119  fNBits=fDigiPar->GetNBits();
120  fDetectedPhotonsPerMeV=fDigiPar->GetDetectedPhotonsPerMeV();
121  fDetectedPhotonsPerMeV_PMT=fDigiPar->GetDetectedPhotonsPerMeV_PMT();
122  fSensitiveAreaAPD=fDigiPar->GetSensitiveAreaAPD();
123  fSensitiveAreaVPT=fDigiPar->GetSensitiveAreaVPT();
124  fQuantumEfficiencyAPD=fDigiPar->GetQuantumEfficiencyAPD();
125  fQuantumEfficiencyVPT=fDigiPar->GetQuantumEfficiencyVPT();
126  fQuantumEfficiencyPMT=fDigiPar->GetQuantumEfficiencyPMT();
127  fExcessNoiseFactorAPD=fDigiPar->GetExcessNoiseFactorAPD();
128  fExcessNoiseFactorVPT=fDigiPar->GetExcessNoiseFactorVPT();
129  fExcessNoiseFactorPMT = fDigiPar->GetExcessNoiseFactorPMT();
130  fIncoherent_elec_noise_width_GeV_APD=fDigiPar->GetIncoherent_elec_noise_width_GeV_APD(); //GeV
131  fIncoherent_elec_noise_width_GeV_VPT=fDigiPar->GetIncoherent_elec_noise_width_GeV_VPT(); //GeV
132  fEnergyRange=fDigiPar->GetEnergyRange(); //GeV
133  fEnergyRangeBW=fDigiPar->GetEnergyRangeBW(); //GeV
134  fFirstSamplePhase=fDigiPar->GetFirstSamplePhase();
135  fNumber_of_samples_in_waveform=fDigiPar->GetNumber_of_samples_in_waveform();
136  fNumber_of_samples_in_waveform_fwd=fDigiPar->GetNumber_of_samples_in_waveform_fwd();
137  fNumber_of_samples_in_waveform_pmt=fDigiPar->GetNumber_of_samples_in_waveform_pmt();
138  fASIC_Shaping_int_time=fDigiPar->GetASIC_Shaping_int_time(); //s
139  fPMT_Shaping_int_time=fDigiPar->GetPMT_Shaping_int_time(); //s
140  fPMT_Shaping_diff_time=fDigiPar->GetPMT_Shaping_diff_time(); //s
141  fFWD_Shaping_int_time=fDigiPar->GetFWD_Shaping_int_time(); //s
142  fFWD_time_constant=fDigiPar->GetFWD_time_constant(); //s
143  fCrystal_time_constant=fDigiPar->GetCrystal_time_constant(); //s
144  fShashlyk_time_constant=fDigiPar->GetShashlyk_time_constant(); //s
145  fShashlykSamplingFactor=fDigiPar->GetShashlykSamplingFactor();
146  fSampleRate=fDigiPar->GetSampleRate();
147  fSampleRate_PMT=fDigiPar->GetSampleRate_PMT();
148  fSampleRate_FWD=fDigiPar->GetSampleRate_FWD();
149  fUse_shaped_noise=fDigiPar->GetUse_shaped_noise();
150  fUse_photon_statistic=fDigiPar->GetUse_photon_statistic();
151  fNoiseAllChannels=fDigiPar->GetNoiseAllChannels();
152 
153  fDigiPar->printParams();
154  // Test how parameters were read from DB.
155  cout<<"EMC digitisation parameters "<<endl;
156  cout<<" nBits "<<fNBits<<endl;
157  cout<<" detectedPhotonsPerMeV "<<fDetectedPhotonsPerMeV<<endl;
158  cout<<" detectedPhotonsPerMeV_PMT "<<fDetectedPhotonsPerMeV_PMT<<endl;
159  cout<<" excessNoiseFactor APD"<<fExcessNoiseFactorAPD<<endl;
160  cout<<" excessNoiseFactor VPT"<<fExcessNoiseFactorVPT<<endl;
161  cout<<" excessNoiseFactor PMT"<<fExcessNoiseFactorPMT<<endl;
162  cout<<" incoherent_elec_noise_width_GeV_APD "<<fIncoherent_elec_noise_width_GeV_APD<<endl;
163  cout<<" incoherent_elec_noise_width_GeV_VPT "<<fIncoherent_elec_noise_width_GeV_VPT<<endl;
164  cout<<" energyRange "<<fEnergyRange<<endl;
165  cout<<" energyRangeBW "<<fEnergyRangeBW<<endl;
166  cout<<" firstSamplePhase "<<fFirstSamplePhase<<endl;
167  cout<<" number_of_samples_in_waveform "<<fNumber_of_samples_in_waveform<<endl;
168  cout<<" number_of_samples_in_waveform_pmt "<<fNumber_of_samples_in_waveform_pmt<<endl;
169  cout<<" number_of_samples_in_waveform_fwd "<<fNumber_of_samples_in_waveform_fwd<<endl;
170  cout<<" ASIC_Shaping_int_time "<<fASIC_Shaping_int_time<<endl;
171  cout<<" PMT_Shaping_int_time "<<fPMT_Shaping_int_time<<endl;
172  cout<<" PMT_Shaping_diff_time "<<fPMT_Shaping_diff_time<<endl;
173  cout<<" crystal_time_constant "<<fCrystal_time_constant<<endl;
174  cout<<" shashlyk_time_constant "<<fShashlyk_time_constant<<endl;
175  cout<<" ShashlykSamplingFactor "<<fShashlykSamplingFactor<<endl;
176  cout<<" sampleRate "<<fSampleRate<<endl;
177  cout<<" sampleRate_PMT "<<fSampleRate_PMT<<endl;
178  cout<<" use_shaped_noise "<<fUse_shaped_noise<<endl;
179  cout<<" use_photon_statistic "<<fUse_photon_statistic<<endl;
180  cout<<" EMC mapper "<<fGeoPar->GetMapperVersion()<<endl;
181 
182  fGeoPar->InitEmcMapper();
183  PndEmcStructure::Instance();
184 
185  //fUse_shaped_noise = 0;
186  //fNumber_of_samples_in_waveform = 64;
187  //fNumber_of_samples_in_waveform_FWD = 64;
188  // Calculate 1 bit resolution (in units of FADC amplitude)
189  PndEmcWaveform *tmpwaveform1=new PndEmcWaveform(0,101010001, fNumber_of_samples_in_waveform);
190  pulseshape1 =new PndEmcAsicPulseshape(fASIC_Shaping_int_time,fCrystal_time_constant);
191 
192  PndEmcWaveform *tmpwaveform2=new PndEmcWaveform(0,101010001, fNumber_of_samples_in_waveform_pmt);
193  pulseshape2 =new PndEmcCRRCPulseshape(fPMT_Shaping_int_time,fPMT_Shaping_diff_time,fShashlyk_time_constant);
194 
195  PndEmcWaveform *tmpwaveform3=new PndEmcWaveform(0,101010001, fNumber_of_samples_in_waveform_fwd);//length 10
196  pulseshape3 =new PndEmcAsicPulseshape(fFWD_Shaping_int_time, fFWD_time_constant);//LNP raw signal, decay time constant 25 microseconds
197 
198 
199  //temp parameters;
200  //fSampleRate_FWD = 100e6;//100 MHz
201 
202  fGevPeakAnalogue = tmpwaveform1->GetScale(fSampleRate, pulseshape1);
203  fGevPeakAnalogue_PMT = tmpwaveform2->GetScale(fSampleRate_PMT, pulseshape2);
204  fGevPeakAnalogue_FWD = tmpwaveform3->GetScale(fSampleRate_FWD, pulseshape3);//forward endcap,100 MHz
205 
206  // cout<<" -I- PndEmcHitsToWaveform::Init:=========== "<<endl;
207  fOneBitResolution=fEnergyRange/((double) (1<<fNBits))*fGevPeakAnalogue;
208  fOneBitResolutionBW=fEnergyRangeBW/((double) (1<<fNBits))*fGevPeakAnalogue;
209  fOneBitResolutionPMT=fEnergyRange/((double) (1<<fNBits))*fGevPeakAnalogue_PMT;
210  fOneBitResolutionFWD=fEnergyRange/((double) (1<<fNBits))*fGevPeakAnalogue_FWD;
211 
212  cout<<" fGevPeakAnalogue= "<<fGevPeakAnalogue<<endl;
213  cout<<" fGevPeakAnalogue_PMT= "<<fGevPeakAnalogue_PMT<<endl;
214  cout<<" fGevPeakAnalogue_FWD= "<<fGevPeakAnalogue_FWD<<endl;
215  cout<<" fOneBitResolution= "<<fOneBitResolution<<endl;
216  cout<<" fOneBitResolutionBW= "<<fOneBitResolutionBW<<endl;
217  cout<<" fOneBitResolutionPMT= "<<fOneBitResolutionPMT<<endl;
218 
219  fFirstADCBinTime=fFirstSamplePhase/fSampleRate;
220  // cout<<" fFirstADCBinTime= "<<fFirstADCBinTime<<endl;
221 
222  // Calculate number of photoelectrons for APD and VPT
223  // The number fDetectedPhotonsPerMeV is the measured number of photoelectrons with PM covering the whole rear surface divided by quantum efficiency of PM (18%)
224  // To estimate Number of photoelectrons in barrel the rare surface is taken equal for all the crystals 745 mm^2, which is average surface, hovewer it varies depending on the type of the crystal
225  // For forward and backward endcap rear surface is equal 26x26=676 mm^2
226  // Therefore the different number of photoelectrons are used with APD for barrel and backward endcap
227  fNPhotoElectronsPerMeVAPDBarrel=fDetectedPhotonsPerMeV*fSensitiveAreaAPD/745.*fQuantumEfficiencyAPD;
228  fNPhotoElectronsPerMeVAPDBWD=fDetectedPhotonsPerMeV*fSensitiveAreaAPD/676.*fQuantumEfficiencyAPD;
229  fNPhotoElectronsPerMeVVPT=fDetectedPhotonsPerMeV*fSensitiveAreaVPT/676.*fQuantumEfficiencyVPT;
230  fNPhotoElectronsPerMeVPMT=fDetectedPhotonsPerMeV_PMT*fQuantumEfficiencyPMT;
231  cout<<" fNPhotoElectronsPerMeVAPDBarrel= "<<fNPhotoElectronsPerMeVAPDBarrel<<endl;
232  cout<<" fNPhotoElectronsPerMeVAPDBWD= "<<fNPhotoElectronsPerMeVAPDBWD<<endl;
233  cout<<" fNPhotoElectronsPerMeVVPT= "<<fNPhotoElectronsPerMeVVPT<<endl;
234  cout<<" fNPhotoElectronsPerMeVPMT= "<<fNPhotoElectronsPerMeVPMT<<endl;
235  //delete pulseshape1;
236  delete tmpwaveform1;
237  //delete pulseshape2;
238  delete tmpwaveform2;
239  //delete pulseshape3;
240  delete tmpwaveform3;
241 
242  nWaveformProduced = 0;
243 
244  return kSUCCESS;
245 }
246 
257 void PndEmcHitsToWaveform::Exec(Option_t*)
258 {
259  TStopwatch timer;
260  if (fVerbose>2){
261  timer.Start();
262  }
263  // Reset output array
264  if(fTimeOrderedWaveform) {
265  if ( ! fDataBuffer ) Fatal("Exec", "No Waveform Data Buffer");
266  }else{
267  if ( ! fWaveformArray ) Fatal("Exec", "No Waveform Array");
268  fWaveformArray->Delete();
269  }
270 
271  Double_t EventTime = FairRootManager::Instance()->GetEventTime();//nano seconds
272  Int_t nHits = fHitArray->GetEntriesFast();
273  Int_t evtNo = FairRun::Instance()->GetEventHeader()->GetMCEntryNumber();
274  if (fVerbose>1){
275  cout<<"**************************************"<<endl;
276  cout<<"Event No. #"<<evtNo<<", EvtTime #"<<EventTime<<std::endl;
277  cout<<"PndEmcHitsToWaveform:: Hit array contains "<<nHits<< " hits"<<endl;
278  cout<<"fDataBuffer size #"<< (fTimeOrderedWaveform ? fDataBuffer->GetNData() : fWaveformArray->GetEntriesFast()) <<std::endl;
279  cout<<"**************************************"<<endl;
280  }
281 
282  // Variable declaration
283  PndEmcHit* theHit = NULL;
284  PndEmcWaveform* theWaveform = NULL;
285  PndEmcHit* tmpHit = NULL;
286  std::set<Int_t> waveformInd;
287  Int_t NumOfSamples;
288 
289  // Loop over PndEmcHits to add them to correspondent waveforms
290  // <set> fWaveformInd contains indexes of detectors for which Waveforms are created
291  //PndEmcAsicPulseshape *pulseshape= new PndEmcAsicPulseshape(fASIC_Shaping_int_time,fCrystal_time_constant);
292  //PndEmcAbsPulseshape *pulseshape2=new PndEmcCRRCPulseshape(fPMT_Shaping_int_time,fPMT_Shaping_diff_time,fShashlyk_time_constant);
293  //PndEmcAbsPulseshape *pulseshape3=new PndEmcAsicPulseshape(10.e-9,25e-6);//fCrystal_time_constant = decay time constant 25 micorseconds
294 
295  Double_t sampleRate;
296  Double_t TimeMin = 99999.;
297  Double_t TimeMax = 0.;
298  Double_t WaveformTimeStamp(0.);
299  Double_t TimeError(0.);
300  Double_t EnergyError(0.);
301  Int_t detId, module, MCTrackID;
302 
303  HowManyHit += nHits;
304  for (Int_t iHit=0; iHit<nHits; iHit++)
305  {
306  theHit = (PndEmcHit*) fHitArray->At(iHit);
307  module = theHit->GetModule();
308 
309  if(module > 5) {
310  std::cout<<" UpdateWaveform: Unknown module number "<<module<<" in EMC digitization. Detector ID = "<<detId<<std::endl;
311  continue;
312  }
313  detId=theHit->GetDetectorID();
314  //if(theHit->GetEnergy() < 0.001) continue;//1MeV
315  TimeError = gRandom->Gaus(0, 0.55 + 5.5*TMath::Exp(-27.7*theHit->GetEnergy()));//ns
316  EnergyError = sqrt(0.01*0.01+0.02*0.02/theHit->GetEnergy());
317  const std::vector<Int_t>& mcTrack = theHit->GetMcList();
318  MCTrackID = mcTrack.size() > 0 ? mcTrack[0] : -1;
319 
320 
321  waveformInd.insert(detId);
322  if(module == 5){
323  sampleRate = fSampleRate_PMT;
324  NumOfSamples = fNumber_of_samples_in_waveform_pmt;
325  }else if(module == 3){
326  sampleRate = fSampleRate_FWD;//100 MHz
327  NumOfSamples = fNumber_of_samples_in_waveform_fwd;
328  }else{
329  sampleRate = fSampleRate;
330  NumOfSamples = fNumber_of_samples_in_waveform;
331  }
332  WaveformTimeStamp = EventTime + theHit->GetTime()*1.0e9 ;//to nano seconds
333  theWaveform = AddWaveform(detId,iHit,NumOfSamples, WaveformTimeStamp, sampleRate, MCTrackID);
334  theWaveform->AddEvt(evtNo);//very important
335  theWaveform->SetTimeStampError(TimeError);
336 
337  //WaveformMap[detId].push_back(Element(evtNo, theWaveform->GetTimeStamp()));
338  if(theWaveform->GetTimeStamp() > TimeMax) TimeMax = theWaveform->GetTimeStamp();
339  if(theWaveform->GetTimeStamp() < TimeMin) TimeMin = theWaveform->GetTimeStamp();
340  // Int_t module_wf = theWaveform->GetModule();
341 
342  switch (module){
343  case 1: // Barrel
344  theWaveform->UpdateWaveform(theHit, fNPhotoElectronsPerMeVAPDBarrel, fUse_photon_statistic, fExcessNoiseFactorAPD, fFirstADCBinTime, fSampleRate, pulseshape1, EnergyError);
345  if (fUse_shaped_noise==0)
346  {
347  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution, EnergyError);
348  }
349  else {
350  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution, pulseshape1, fFirstADCBinTime, fSampleRate, EnergyError);
351  }
352  break;
353  case 2: // Barrel
354  theWaveform->UpdateWaveform(theHit, fNPhotoElectronsPerMeVAPDBarrel, fUse_photon_statistic, fExcessNoiseFactorAPD, fFirstADCBinTime, fSampleRate, pulseshape1, EnergyError);
355  if (fUse_shaped_noise==0)
356  {
357  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution, EnergyError);
358  }
359  else {
360  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution, pulseshape1, fFirstADCBinTime, fSampleRate, EnergyError);
361  }
362  break;
363  case 3: // Fwd endcap
364  theWaveform->UpdateWaveform(theHit, fNPhotoElectronsPerMeVVPT, fUse_photon_statistic, fExcessNoiseFactorVPT, fFirstADCBinTime, fSampleRate_FWD, pulseshape3, EnergyError);//use a different shape
365  if (fUse_shaped_noise==0)
366  {
367  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_VPT*fGevPeakAnalogue_FWD,fOneBitResolutionFWD, EnergyError);
368  }
369  else {
370  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_VPT*fGevPeakAnalogue_FWD,fOneBitResolutionFWD, pulseshape3, fFirstADCBinTime, fSampleRate_FWD, EnergyError);
371  }
372  break;
373  case 4: // Bwd endcap
374  theWaveform->UpdateWaveform(theHit, fNPhotoElectronsPerMeVVPT, fUse_photon_statistic, fExcessNoiseFactorVPT, fFirstADCBinTime, fSampleRate, pulseshape1, EnergyError);
375  if (fUse_shaped_noise==0)
376  {
377  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolutionBW, EnergyError);
378  }
379  else {
380  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolutionBW, pulseshape1, fFirstADCBinTime, fSampleRate, EnergyError);
381  }
382  break;
383  case 5: // Shashlyk calorimetr
384  tmpHit = theHit;
385  tmpHit->SetEnergy(tmpHit->GetEnergy()*fShashlykSamplingFactor);
386  EnergyError = sqrt(0.01*0.01+0.02*0.02/tmpHit->GetEnergy());
387  theWaveform->UpdateWaveform(theHit, fNPhotoElectronsPerMeVPMT, fUse_photon_statistic,fExcessNoiseFactorPMT, fFirstADCBinTime, fSampleRate_PMT, pulseshape2, EnergyError);
388  if (fUse_shaped_noise==0)
389  {
390  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue_PMT,fOneBitResolutionPMT, EnergyError);
391  }
392  else {
393  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue_PMT,fOneBitResolutionPMT, pulseshape2, fFirstADCBinTime, fSampleRate_PMT, EnergyError);
394  }
395  break;
396  default:
397  std::cout<<" UpdateWaveform: Unknown module number "<<module<<" in EMC digitization. Detector ID = "<<detId<<std::endl;
398  //abort();
399  }
400  if(fTimeOrderedWaveform){
401  fDataBuffer->FillNewData(theWaveform, theWaveform->GetTimeStamp(), theWaveform->GetActiveTime());
402  }
403  }
404  // Produce waveforms in all the crystals, not only where hits took place
405  // Since it is time consuming, by default it is off
406  if (fNoiseAllChannels)
407  {
408  Int_t detId_tmp, modId_tmp;
409  std::map<Int_t,PndEmcTwoCoordIndex*> intTwoCoordMap = PndEmcMapper::Instance()->GetTciMap();
410  for(std::map<Int_t,PndEmcTwoCoordIndex* >::iterator iter = intTwoCoordMap.begin();
411  iter != intTwoCoordMap.end(); ++iter){
412  detId_tmp=(*iter).first;
413  modId_tmp = detId_tmp/100000000;
414  if(modId_tmp == 5){
415  NumOfSamples = fNumber_of_samples_in_waveform_pmt;
416  sampleRate = fSampleRate_PMT;
417  }else if(modId_tmp == 3){
418  NumOfSamples = fNumber_of_samples_in_waveform_fwd;
419  sampleRate = fSampleRate_FWD;
420  }else{
421  NumOfSamples = fNumber_of_samples_in_waveform;
422  sampleRate = fSampleRate;
423  }
424 
425  if (waveformInd.insert(detId_tmp).second){
426  WaveformTimeStamp = gRandom->Uniform(TimeMin, TimeMax);
427  theWaveform = AddWaveform(detId_tmp,-1,NumOfSamples, WaveformTimeStamp, sampleRate, -1); // -1 correponds to Waveform produced not from EmcHit but from Noise
428  switch (modId_tmp){
429  case 1: // Barrel
430  if (fUse_shaped_noise==0)
431  {
432  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution);
433  }
434  else {
435  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution, pulseshape1, fFirstSamplePhase, fSampleRate);
436  }
437  break;
438  case 2: // Barrel
439  if (fUse_shaped_noise==0)
440  {
441  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution);
442  }
443  else {
444  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution, pulseshape1, fFirstSamplePhase, fSampleRate);
445  }
446  break;
447  case 3: // Fwd endcap
448  if (fUse_shaped_noise==0)
449  {
450  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_VPT*fGevPeakAnalogue_FWD,fOneBitResolutionFWD);
451  }
452  else {
453  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_VPT*fGevPeakAnalogue_FWD,fOneBitResolutionFWD, pulseshape3, fFirstSamplePhase, fSampleRate_FWD);
454  }
455  break;
456  case 4: // Bwd endcap
457  if (fUse_shaped_noise==0)
458  {
459  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolutionBW);
460  }
461  else {
462  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolutionBW, pulseshape1, fFirstSamplePhase, fSampleRate);
463  }
464  break;
465  case 5: // Shashlyk calorimetr
466  if (fUse_shaped_noise==0)
467  {
468  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue_PMT,fOneBitResolutionPMT);
469  }
470  else {
471  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue_PMT,fOneBitResolutionPMT, pulseshape2, fFirstSamplePhase, fSampleRate_PMT);
472  }
473  break;
474  default:
475  std::cout<<" UpdateWaveform: Unknown module number "<<modId_tmp<<" in EMC digitization. Detector ID = "<<detId_tmp<<std::endl;
476  //abort();
477  }
478  if(fTimeOrderedWaveform){
479  fDataBuffer->FillNewData(theWaveform, theWaveform->GetTimeStamp(), theWaveform->GetActiveTime());
480  }
481  }
482  }
483  }
484 
485  // Add electronic noise
486  // There are two options how to add noise (before and after shaping)
487 
488  /*Int_t nWf = fWaveformArray->GetEntriesFast();
489  if (fVerbose>2){
490  cout << "Number of waveforms processed= "<<nWf<<endl;
491  }
492 
493  for (Int_t iWf=0; iWf<nWf; iWf++) {
494  theWaveform = (PndEmcWaveform*) fWaveformArray->At(iWf);
495  Int_t module = theWaveform->GetModule();
496  switch (module){
497  case 1: // Barrel
498  if (fUse_shaped_noise==0)
499  {
500  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution);
501  }
502  else {
503  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution, pulseshape, fFirstSamplePhase, fSampleRate);
504  }
505  break;
506  case 2: // Barrel
507  if (fUse_shaped_noise==0)
508  {
509  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution);
510  }
511  else {
512  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolution, pulseshape, fFirstSamplePhase, fSampleRate);
513  }
514  break;
515  case 3: // FWD Endcap
516  if (fUse_shaped_noise==0)
517  {
518  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_VPT*fGevPeakAnalogue,fOneBitResolution);
519  }
520  else {
521  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_VPT*fGevPeakAnalogue,fOneBitResolution, pulseshape, fFirstSamplePhase, fSampleRate);
522  }
523  break;
524  case 4: // BWD Endcap
525  if (fUse_shaped_noise==0)
526  {
527  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolutionBW);
528  }
529  else {
530  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue,fOneBitResolutionBW, pulseshape, fFirstSamplePhase, fSampleRate);
531  }
532  break;
533  case 5: // shashlyk calorimetr
534  if (fUse_shaped_noise==0)
535  {
536  theWaveform->AddElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue_PMT,fOneBitResolutionPMT);
537  }
538  else {
539  theWaveform->AddShapedElecNoiseAndDigitise(fIncoherent_elec_noise_width_GeV_APD*fGevPeakAnalogue_PMT,fOneBitResolutionPMT, pulseshape2, fFirstSamplePhase, fSampleRate_PMT);
540  }
541  break;
542  default:
543  std::cout<<"Add Noise: Unknown module number in EMC digitization"<<std::endl;
544  //abort();
545  }
546  }*/
547 
548  //std::cout<<"nWaveformProduced = "<<nWaveformProduced<<std::endl;
549  //delete pulseshape;
550  //delete pulseshape2;
551  //delete pulseshape3;
552 
553  if (fVerbose>2){
554  timer.Stop();
555  Double_t rtime = timer.RealTime();
556  Double_t ctime = timer.CpuTime();
557  cout << "PndEmcHitsToWaveform, Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
558  }
559 }
560 
561 void PndEmcHitsToWaveform::SetParContainers() {
562 
563  // Get run and runtime database
564  FairRun* run = FairRun::Instance();
565  if ( ! run ) Fatal("SetParContainers", "No analysis run");
566 
567  FairRuntimeDb* db = run->GetRuntimeDb();
568  if ( ! db ) Fatal("SetParContainers", "No runtime database");
569 
570  // Get Emc geometry parameter container
571  fGeoPar = (PndEmcGeoPar*) db->getContainer("PndEmcGeoPar");
572  // Get Emc digitisation parameter container
573  fDigiPar = (PndEmcDigiPar*) db->getContainer("PndEmcDigiPar");
574 
575 }
576 
577 // ----- Private method AddWaveform --------------------------------------------
578 /*PndEmcWaveform* PndEmcHitsToWaveform::AddWaveform(Int_t detID, Int_t iHit,Int_t numOfSamples){
579  TClonesArray& clref = *fWaveformArray;
580  Int_t size = clref.GetEntriesFast();
581  return new(clref[size]) PndEmcWaveform(0,detID,numOfSamples,iHit);
582  }*/
594 PndEmcWaveform* PndEmcHitsToWaveform::AddWaveform(Int_t detID,
595  Int_t iHit,
596  Int_t numOfSamples,
597  Double_t timeStamp,
598  Double_t sampleRate,
599  Int_t MCTrackID)
600 {
601  nWaveformProduced ++;
602  PndEmcWaveform* thisWave(0);
603  if(fTimeOrderedWaveform){
604  thisWave = new PndEmcWaveform(MCTrackID,detID,sampleRate, numOfSamples, iHit, timeStamp);
605  }else{
606  TClonesArray& clref = *fWaveformArray;
607  Int_t size = clref.GetEntriesFast();
608  thisWave = new(clref[size]) PndEmcWaveform(MCTrackID,detID,sampleRate, numOfSamples,iHit, timeStamp);
609  }
610  return thisWave;
611 }
612 
613 void PndEmcHitsToWaveform::SetStorageOfData(Bool_t val)
614 {
615  SetPersistency(val);
616  return;
617 }
626 void PndEmcHitsToWaveform::FinishTask()
627 {
628  std::cout<<"==================================================="<<std::endl;
629  std::cout<<"PndEmcHitsToWaveform::FinishTask"<<std::endl;
630  std::cout<<"***************************************************"<<std::endl;
631  std::cout<<"Read Hits# "<<HowManyHit<<std::endl;
632  std::cout<<"Produc waveforms# "<<nWaveformProduced<<std::endl;
633  std::cout<<"***************************************************"<<std::endl;
634 
635  if(fTimeOrderedWaveform){
636  fDataBuffer->Write();
637  }
638 }
639 
640 ClassImp(PndEmcHitsToWaveform)
PndEmcHitsToWaveform
Takes list of PndEmcHits and creates PndEmcWaveform.
Definition: PndEmcHitsToWaveform.h:41
PndEmcHitsToWaveform::Init
virtual InitStatus Init()
Init Task.
Definition: PndEmcHitsToWaveform.cxx:86
fVerbose
int fVerbose
Definition: poormantracks.C:24
PndEmcDigiPar::GetQuantumEfficiencyVPT
Double_t GetQuantumEfficiencyVPT()
Definition: PndEmcDigiPar.h:21
PndEmcHitsToWaveform::PndEmcHitsToWaveform
PndEmcHitsToWaveform(Int_t verbose=0, Bool_t storewaves=kTRUE)
Definition: PndEmcHitsToWaveform.cxx:55
PndEmcDigiPar::GetNBits
Int_t GetNBits()
Definition: PndEmcDigiPar.h:43
PndEmcWaveform::GetActiveTime
Double_t GetActiveTime() const
Definition: PndEmcWaveform.h:132
PndEmcDigiPar::GetUse_photon_statistic
Int_t GetUse_photon_statistic()
Definition: PndEmcDigiPar.h:48
PndEmcWaveform::AddShapedElecNoiseAndDigitise
void AddShapedElecNoiseAndDigitise(Double_t noise_width, Double_t oneBitResolution, PndEmcAbsPulseshape *pulseshape, Double_t firstADCBinTime, Double_t sampleRate, Double_t=0)
Definition: PndEmcWaveform.cxx:291
PndEmcDigiPar::GetSampleRate
Double_t GetSampleRate()
Definition: PndEmcDigiPar.h:39
PndEmcDigiPar::GetPMT_Shaping_diff_time
Double_t GetPMT_Shaping_diff_time()
Definition: PndEmcDigiPar.h:33
PndEmcDigiPar::GetQuantumEfficiencyPMT
Double_t GetQuantumEfficiencyPMT()
Definition: PndEmcDigiPar.h:22
PndEmcHitsToWaveform::fIncoherent_elec_noise_width_GeV_VPT
Double_t fIncoherent_elec_noise_width_GeV_VPT
Definition: PndEmcHitsToWaveform.h:98
run
Int_t run
Definition: autocutx.C:47
PndEmcDigiPar::GetIncoherent_elec_noise_width_GeV_VPT
Double_t GetIncoherent_elec_noise_width_GeV_VPT()
Definition: PndEmcDigiPar.h:27
PndEmcDigiPar::GetUse_shaped_noise
Int_t GetUse_shaped_noise()
Definition: PndEmcDigiPar.h:47
PndEmcDigiPar::GetNumber_of_samples_in_waveform
Int_t GetNumber_of_samples_in_waveform()
Definition: PndEmcDigiPar.h:44
PndEmcDigiPar::GetQuantumEfficiencyAPD
Double_t GetQuantumEfficiencyAPD()
Definition: PndEmcDigiPar.h:20
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
val
Double_t val[nBoxes][nFEBox]
Definition: createCalib.C:11
PndEmcMapper::GetTciMap
const std::map< Int_t, PndEmcTwoCoordIndex * > & GetTciMap()
Definition: PndEmcMapper.cxx:441
PndEmcDigiPar::GetIncoherent_elec_noise_width_GeV_APD
Double_t GetIncoherent_elec_noise_width_GeV_APD()
Definition: PndEmcDigiPar.h:26
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
PndEmcDigiPar::GetEnergyRangeBW
Double_t GetEnergyRangeBW()
Definition: PndEmcDigiPar.h:29
PndEmcHitsToWaveform::fWaveformArray
TClonesArray * fWaveformArray
Definition: PndEmcHitsToWaveform.h:73
PndEmcDigiPar::GetCrystal_time_constant
Double_t GetCrystal_time_constant()
Definition: PndEmcDigiPar.h:36
PndEmcDigiPar::GetFirstSamplePhase
Double_t GetFirstSamplePhase()
Definition: PndEmcDigiPar.h:30
PndPersistencyTask::SetPersistency
void SetPersistency(Bool_t val=kTRUE)
Definition: PndPersistencyTask.h:31
PndEmcHitsToWaveform::pulseshape2
PndEmcAbsPulseshape * pulseshape2
Definition: PndEmcHitsToWaveform.h:140
PndEmcHitsToWaveform::fDataBuffer
PndEmcWaveformWriteoutBuffer * fDataBuffer
Definition: PndEmcHitsToWaveform.h:74
PndEmcDigiPar::GetASIC_Shaping_int_time
Double_t GetASIC_Shaping_int_time()
Definition: PndEmcDigiPar.h:31
PndEmcGeoPar::GetMapperVersion
Int_t GetMapperVersion()
Definition: PndEmcGeoPar.h:20
PndEmcCRRCPulseshape
Pulseshape from an CRRC-Shaper.
Definition: PndEmcCRRCPulseshape.h:25
PndEmcDigiPar::GetSampleRate_FWD
Double_t GetSampleRate_FWD()
Definition: PndEmcDigiPar.h:41
PndEmcGeoPar::InitEmcMapper
void InitEmcMapper()
Definition: PndEmcGeoPar.cxx:41
PndEmcDigiPar::GetShashlykSamplingFactor
Double_t GetShashlykSamplingFactor()
Definition: PndEmcDigiPar.h:38
PndEmcDigiPar::GetSensitiveAreaAPD
Double_t GetSensitiveAreaAPD()
Definition: PndEmcDigiPar.h:18
PndEmcHitsToWaveform::SetStorageOfData
void SetStorageOfData(Bool_t val)
Definition: PndEmcHitsToWaveform.cxx:613
PndEmcWaveform::GetScale
Double_t GetScale(Double_t sampleRate, PndEmcAbsPulseshape *pulseshape) const
Definition: PndEmcWaveform.cxx:143
PndEmcHit::GetMcList
const std::vector< Int_t > & GetMcList()
Definition: PndEmcHit.h:69
nHits
int nHits
Definition: RiemannTest.C:16
PndEmcAsicPulseshape
Pulseshape from an APFEL ASIC preamplifier shaper.
Definition: PndEmcAsicPulseshape.h:31
PndEmcHitsToWaveform::fNPhotoElectronsPerMeVAPDBarrel
Double_t fNPhotoElectronsPerMeVAPDBarrel
Definition: PndEmcHitsToWaveform.h:85
PndEmcHitsToWaveform::fIncoherent_elec_noise_width_GeV_APD
Double_t fIncoherent_elec_noise_width_GeV_APD
Definition: PndEmcHitsToWaveform.h:97
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndEmcHit::GetEnergy
virtual Double_t GetEnergy() const
Definition: PndEmcHit.h:54
timer
TStopwatch timer
Definition: hit_dirc.C:51
PndEmcDigiPar
parameter set of Emc digitisation
Definition: PndEmcDigiPar.h:12
PndEmcDigiPar::GetFWD_Shaping_int_time
Double_t GetFWD_Shaping_int_time()
Definition: PndEmcDigiPar.h:34
PndEmcWaveform::AddElecNoiseAndDigitise
void AddElecNoiseAndDigitise(Double_t, Double_t, Double_t=0)
Definition: PndEmcWaveform.cxx:259
PndEmcDigiPar::GetEnergyRange
Double_t GetEnergyRange()
Definition: PndEmcDigiPar.h:28
PndEmcWaveform::UpdateWaveform
void UpdateWaveform(PndEmcHit *hit, Double_t pePerMeV, Bool_t usePhotonStatistic, Double_t excessNoiseFactor, Double_t firstADCBinTime, Double_t sampleRate, PndEmcAbsPulseshape *pulseshape, Double_t=0)
Definition: PndEmcWaveform.cxx:173
PndEmcHit::SetEnergy
virtual void SetEnergy(Double32_t energy)
Definition: PndEmcHit.h:50
PndEmcDigiPar::GetDetectedPhotonsPerMeV
Double_t GetDetectedPhotonsPerMeV()
Definition: PndEmcDigiPar.h:16
PndEmcWaveform
represents a simulated waveform in an emc crystal
Definition: PndEmcWaveform.h:56
ctime
Double_t ctime
Definition: hit_dirc.C:114
PndEmcHitsToWaveform::pulseshape1
PndEmcAbsPulseshape * pulseshape1
Definition: PndEmcHitsToWaveform.h:139
PndEmcDigiPar::GetNumber_of_samples_in_waveform_fwd
Int_t GetNumber_of_samples_in_waveform_fwd()
Definition: PndEmcDigiPar.h:46
PndEmcHitsToWaveform::Exec
virtual void Exec(Option_t *opt)
Runs the task.
Definition: PndEmcHitsToWaveform.cxx:257
PndEmcHit::GetTime
virtual Double_t GetTime() const
Definition: PndEmcHit.h:55
PndEmcHit
represents the deposited energy of one emc crystal from simulation
Definition: PndEmcHit.h:26
PndEmcDigiPar::GetExcessNoiseFactorVPT
Double_t GetExcessNoiseFactorVPT()
Definition: PndEmcDigiPar.h:24
PndEmcHitsToWaveform::AddWaveform
PndEmcWaveform * AddWaveform(Int_t detId, Int_t hitIndex, Int_t numOfSamples, Double_t timeStamp, Double_t sampleRate, Int_t MCTrackID)
Create a new PndEmcWaveform from the passed parameters.
Definition: PndEmcHitsToWaveform.cxx:594
PndEmcDigiPar::GetDetectedPhotonsPerMeV_PMT
Double_t GetDetectedPhotonsPerMeV_PMT()
Definition: PndEmcDigiPar.h:17
PndEmcDigiPar::GetExcessNoiseFactorAPD
Double_t GetExcessNoiseFactorAPD()
Definition: PndEmcDigiPar.h:23
PndEmcDigiPar::GetFWD_time_constant
Double_t GetFWD_time_constant()
Definition: PndEmcDigiPar.h:35
ClassImp
ClassImp(PndAnaContFact)
PndEmcHitsToWaveform::pulseshape3
PndEmcAbsPulseshape * pulseshape3
Definition: PndEmcHitsToWaveform.h:141
PndEmcWaveform::AddEvt
void AddEvt(Int_t evtNo)
Definition: PndEmcWaveform.h:135
PndEmcDigiPar::GetSampleRate_PMT
Double_t GetSampleRate_PMT()
Definition: PndEmcDigiPar.h:40
PndEmcStructure::Instance
static PndEmcStructure * Instance()
Definition: PndEmcStructure.cxx:70
PndEmcHitsToWaveform::FinishTask
void FinishTask()
Called at end of task.
Definition: PndEmcHitsToWaveform.cxx:626
rtime
Double_t rtime
Definition: hit_dirc.C:113
PndEmcHit::GetModule
Short_t GetModule() const
Definition: PndEmcHit.h:58
PndEmcDigiPar::GetShashlyk_time_constant
Double_t GetShashlyk_time_constant()
Definition: PndEmcDigiPar.h:37
PndEmcMapper::Instance
static PndEmcMapper * Instance()
Definition: PndEmcMapper.cxx:80
PndEmcDigiPar::GetNumber_of_samples_in_waveform_pmt
Int_t GetNumber_of_samples_in_waveform_pmt()
Definition: PndEmcDigiPar.h:45
PndEmcDigiPar::GetPMT_Shaping_int_time
Double_t GetPMT_Shaping_int_time()
Definition: PndEmcDigiPar.h:32
PndEmcDigiPar::GetSensitiveAreaVPT
Double_t GetSensitiveAreaVPT()
Definition: PndEmcDigiPar.h:19
PndEmcDigiPar::GetExcessNoiseFactorPMT
Double_t GetExcessNoiseFactorPMT()
Definition: PndEmcDigiPar.h:25
PndEmcDigiPar::GetNoiseAllChannels
Int_t GetNoiseAllChannels()
Definition: PndEmcDigiPar.h:50
Generated on Wed Apr 3 2019 08:43:54 for FairRoot/PandaRoot by   doxygen 1.8.5