Macros used to produce the OCDB entries for MC. One for default (p-p) and one for...
[u/mrichter/AliRoot.git] / VZERO / AliVZERODigitizer.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15  
16 /* $Id$ */
17
18 ///_________________________________________________________________________
19 ///
20 /// This class constructs Digits out of Hits
21 ///
22 ///
23
24 // --- Standard library ---
25
26 // --- ROOT system ---
27 #include <TMath.h>
28 #include <TTree.h>
29 #include <TMap.h>
30 #include <TGeoManager.h>
31 #include <TGeoPhysicalNode.h>
32 #include <AliGeomManager.h>
33 #include <TRandom.h>
34 #include <TF1.h>
35
36 // --- AliRoot header files ---
37 #include "AliVZEROConst.h"
38 #include "AliRun.h"
39 #include "AliVZERO.h"
40 #include "AliVZEROhit.h"
41 #include "AliRunLoader.h"
42 #include "AliLoader.h"
43 #include "AliGRPObject.h"
44 #include "AliRunDigitizer.h"
45 #include "AliCDBManager.h"
46 #include "AliCDBStorage.h"
47 #include "AliCDBEntry.h"
48 #include "AliVZEROCalibData.h"
49
50 #include "AliVZEROdigit.h"
51 #include "AliVZERODigitizer.h"
52
53 ClassImp(AliVZERODigitizer)
54
55  AliVZERODigitizer::AliVZERODigitizer()
56                    :AliDigitizer(),
57                     fCalibData(GetCalibData()),
58                     fPhotoCathodeEfficiency(0.18),
59                     fPMVoltage(768.0),
60                     fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)),
61                     fNdigits(0),
62                     fDigits(0),
63                     fCollisionMode(0),
64                     fBeamEnergy(0.),
65                     fSignalShape(NULL)
66    
67 {
68   // default constructor
69
70 //    fNdigits = 0;
71 //    fDigits  = 0;
72 //   
73 //    fPhotoCathodeEfficiency =   0.18;
74 //    fPMVoltage              =  768.0;
75 //    fPMGain = TMath::Power((fPMVoltage / 112.5) ,7.04277); 
76    
77 //   fCalibData = GetCalibData();
78
79   fSignalShape = new TF1("VZEROSignalShape",this,&AliVZERODigitizer::SignalShape,0,200,6,"AliVZERODigitizer","SignalShape");
80   fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,2.9,6.40982,3.69339e-01);
81 }
82
83 //____________________________________________________________________________ 
84   AliVZERODigitizer::AliVZERODigitizer(AliRunDigitizer* manager)
85                     :AliDigitizer(manager),
86                      fCalibData(GetCalibData()),
87                      fPhotoCathodeEfficiency(0.18),
88                      fPMVoltage(768.0),
89                      fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)),
90                      fNdigits(0),
91                      fDigits(0),
92                      fCollisionMode(0),
93                      fBeamEnergy(0.),
94                      fSignalShape(NULL)
95                                         
96 {
97   // constructor
98   
99 //   fNdigits = 0;
100 //   fDigits  = 0;
101 //   
102 //   fPhotoCathodeEfficiency =   0.18;
103 //   fPMVoltage              =  768.0;
104 //   fPMGain = TMath::Power( (fPMVoltage / 112.5) ,7.04277 );
105   
106 //  fCalibData = GetCalibData();
107  
108   fSignalShape = new TF1("VZEROSignalShape",this,&AliVZERODigitizer::SignalShape,0,200,6,"AliVZERODigitizer","SignalShape");
109   fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,2.9,6.40982,3.69339e-01);
110   
111 }
112            
113 //____________________________________________________________________________ 
114   AliVZERODigitizer::~AliVZERODigitizer()
115 {
116   // destructor
117   
118   if (fDigits) {
119     fDigits->Delete();
120     delete fDigits;
121     fDigits=0; 
122   }
123
124   if (fSignalShape) {
125     delete fSignalShape;
126     fSignalShape = NULL;
127   }
128 }
129
130 //_____________________________________________________________________________
131 Bool_t AliVZERODigitizer::Init()
132 {
133   // Initialises the digitizer
134
135   // Initialises the Digit array
136   fDigits = new TClonesArray ("AliVZEROdigit", 1000);
137   
138   //  TGeoHMatrix *im = AliGeomManager::GetMatrix("VZERO/V0C");
139   //  im->Print();
140
141   GetCollisionMode();
142   return kTRUE;
143 }
144
145 //____________________________________________________________________________
146 void AliVZERODigitizer::Exec(Option_t* /*option*/) 
147 {   
148   // Creates digits from hits
149      
150   Float_t     map[80];    // 48 values on V0C + 32 on V0A
151 //  Int_t       pmNumber[80];
152   Float_t     adc[64];    // 32 PMs on V0C + 32 PMs on V0A
153   Float_t     time[80], time_ref[80], time2[64];
154   Float_t     adc_gain[80]; 
155   Float_t     adc_pedestal[64],adc_sigma[64];    
156   fNdigits     =    0;  
157   Float_t     pmGain_smeared[64];  
158   Float_t     cPM[80];
159   
160   // Smearing of the PM tubes intrinsic characteristics
161   
162   for(Int_t ii=0; ii<64; ii++) {
163       pmGain_smeared[ii] = gRandom->Gaus(fPMGain, fPMGain/5); }
164               
165   // Retrieval of ADC gain values and pedestal information from local CDB 
166   // I use only the first 64th values of the calibration array in CDB 
167   // as I have no beam burst structure - odd or even beam burst number
168   
169   // Reminder : We have 16 scintillating cells mounted on 8 PMs 
170   // on Ring 3 and Ring 4 in V0C -  added to produce  ADC outputs 
171   // on these rings... 
172    
173   for(Int_t i=0; i<16; i++) { 
174         adc_gain[i] = fCalibData->GetGain(i); 
175         cPM[i]      = fPhotoCathodeEfficiency * pmGain_smeared[i];
176   }
177   
178   for(Int_t j=16; j<48; j=j+2) { 
179         Int_t i=(j+17)/2;
180         adc_gain[j]   = fCalibData->GetGain(i);     
181         adc_gain[j+1] = fCalibData->GetGain(i); 
182         cPM[j]        = fPhotoCathodeEfficiency * pmGain_smeared[i];   
183         cPM[j+1]      = fPhotoCathodeEfficiency * pmGain_smeared[i]; 
184   }
185             
186   for(Int_t i=48; i<80; i++){ 
187         adc_gain[i] = fCalibData->GetGain(i-16); 
188         cPM[i]      = fPhotoCathodeEfficiency * pmGain_smeared[i-16];
189   };
190   
191   for(Int_t  i=0; i<64; i++){ 
192           adc_pedestal[i] = fCalibData->GetPedestal(i);
193           adc_sigma[i]    = fCalibData->GetSigma(i); 
194   }; 
195                                 
196 //  for(Int_t i=0; i<64; i++) { printf(" i = %d pedestal = %f sigma = %f \n\n", 
197 //                                       i, adc_pedestal[i], adc_sigma[i] );} 
198             
199   AliRunLoader* outRunLoader =  AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());    
200   if (!outRunLoader) {
201     Error("Exec", "Can not get output Run Loader");
202     return;
203   }
204     
205   AliLoader* outLoader = outRunLoader->GetLoader("VZEROLoader");
206
207   if (!outLoader) {
208     Error("Exec", "Can not get output VZERO Loader");
209     return;
210   }
211
212   const char* mode = "update";
213   if(outRunLoader->GetEventNumber() == 0) mode = "recreate";
214   outLoader->LoadDigits(mode);
215
216   if (!outLoader->TreeD()) outLoader->MakeTree("D");
217   outLoader->MakeDigitsContainer();
218   TTree* treeD  = outLoader->TreeD();
219   Int_t bufsize = 16000;
220   treeD->Branch("VZERODigit", &fDigits, bufsize); 
221   
222   for (Int_t iInput = 0; iInput < fManager->GetNinputs(); iInput++) {
223      AliRunLoader* runLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
224      AliLoader* loader = runLoader->GetLoader("VZEROLoader");
225      if (!loader) {
226        Error("Exec", "Can not get VZERO Loader for input %d", iInput);
227        continue;
228          }
229       
230      if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
231
232      AliVZERO* vzero = (AliVZERO*) runLoader->GetAliRun()->GetDetector("VZERO");
233      if (!vzero) {
234        Error("Exec", "No VZERO detector for input %d", iInput);
235        continue;
236          }
237       
238      loader->LoadHits();
239      TTree* treeH = loader->TreeH();
240      if (!treeH) {
241        Error("Exec", "Cannot get TreeH for input %d", iInput);
242        continue; 
243          }
244        
245      for(Int_t i=0; i<80; i++) {map[i] = 0; time[i] = 0.0;}
246      
247      TClonesArray* hits = vzero->Hits();
248              
249 //  Now makes Digits from hits
250      
251      Int_t nTracks = (Int_t) treeH->GetEntries();
252      for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
253          for (Int_t i=0; i<80; i++) {time_ref[i] = 999999.0;}   
254          vzero->ResetHits();
255          treeH->GetEvent(iTrack);
256          Int_t nHits = hits->GetEntriesFast();
257          for (Int_t iHit = 0; iHit < nHits; iHit++) {
258                          AliVZEROhit* hit = (AliVZEROhit *)hits->UncheckedAt(iHit);
259                          Int_t nPhot = hit->Nphot();
260                          Int_t cell  = hit->Cell();                          
261                          map[cell] += Float_t(nPhot);
262                          Float_t dt_scintillator = gRandom->Gaus(0,0.7);
263                          Float_t t = dt_scintillator + 1e9*hit->Tof();
264                          if (t > 0.0) {
265                                  if(t < time_ref[cell]) time_ref[cell] = t;
266                                  time[cell] = TMath::Min(t,time_ref[cell]); 
267                          }
268          }           // hit   loop      
269      }             // track loop
270
271      loader->UnloadHits();
272
273   }               // input loop
274
275 // Now builds the scintillator cell response (80 cells i.e. 80 responses)
276          
277    for (Int_t i=0; i<80; i++) {    
278         Float_t q1 = Float_t ( map[i] )* cPM[i] * kQe;
279         Float_t noise = gRandom->Gaus(10.5,3.22);
280         Float_t pmResponse  =  q1/kC*TMath::Power(ktheta/kthau,1/(1-ktheta/kthau)) 
281         + noise*1e-3;   
282         if(fCollisionMode >0) adc_gain[i] = adc_gain[i]/70.0; // reduce dynamics in Ion Collision Mode
283         map[i] =  pmResponse * adc_gain[i];
284         Float_t MIP = 1.0/fCalibData->GetMIPperADC(GetPMNumber(i));
285         if(fCollisionMode >0) MIP=2.0;
286 //      printf("cell = %d,  ADC = %d, TDC = %f \n",i,map[i], time[i]*10.0 );
287         if(map[i] > (MIP/2.) )
288                   {map[i] = gRandom->Gaus(map[i], (MIP/6.) );}
289    }
290       
291 // Now transforms 80 cell responses into 64 photomultiplier responses
292 // Also adds the ADC pedestals taken out of the calibration data base
293         
294    for (Int_t j=0; j<16; j++){
295         adc[j]  = map [j] + gRandom->Gaus(adc_pedestal[j], adc_sigma[j]);
296         time2[j]= time[j];}
297         
298    for (Int_t j=48; j<80; j++){
299         adc[j-16]  = map [j] + gRandom->Gaus(adc_pedestal[j-16],adc_sigma[j-16]);
300         time2[j-16]= time[j]; }
301         
302    for (Int_t j=0; j<16; j++){
303         adc[16+j] = map [16+2*j]+ map [16+2*j+1] + gRandom->Gaus(adc_pedestal[16+j], adc_sigma[16+j]);
304         Float_t min_time = TMath::Min(time [16+2*j],time [16+2*j+1]);
305         time2[16+j] = min_time;
306         if(min_time==0.0){time2[16+j]=TMath::Max(time[16+2*j],time[16+2*j+1]);}
307    }
308         
309
310 // Now add digits to the digit Tree 
311         
312    for (Int_t i=0; i<64; i++) {      
313       if(adc[i] > 0) {
314 //           printf(" Event, cell, adc, tof = %d %d %d %d\n", 
315 //                    outRunLoader->GetEventNumber(),i, adc[i], Int_t((time2[i]*10.0) +0.5));
316 //           multiply by 10 to have 100 ps per channel :
317                   
318                   AddDigit(i, adc[i], time2[i] ) ;
319           }      
320    }
321   treeD->Fill();
322   outLoader->WriteDigits("OVERWRITE");  
323   outLoader->UnloadDigits();     
324   ResetDigit();
325 }
326
327 //____________________________________________________________________________
328 void AliVZERODigitizer::AddDigit(Int_t PMnumber, Float_t adc, Float_t time) 
329  { 
330  
331 // Adds Digit 
332  
333   TClonesArray &ldigits = *fDigits;  
334   Bool_t integrator;
335   if (((Int_t) gRandom->Uniform(2))<1) integrator = kFALSE;
336   else integrator = kTRUE;
337          
338   new(ldigits[fNdigits++]) AliVZEROdigit(PMnumber,adc,time,0,kFALSE,kFALSE,integrator);
339          
340 }
341 //____________________________________________________________________________
342 void AliVZERODigitizer::ResetDigit()
343 {
344
345 // Clears Digits
346
347   fNdigits = 0;
348   if (fDigits) fDigits->Delete();
349 }
350
351 //____________________________________________________________________________
352 void AliVZERODigitizer::GetCollisionMode()
353 {
354 // Retrieves the collision mode from GRP data
355
356 // Initialization of the GRP entry 
357
358    Int_t run = AliCDBManager::Instance()->GetRun();
359   
360 //   printf("\n ++++++ Run Number retrieved as %d \n",run); 
361  
362   AliCDBEntry*  entry = AliCDBManager::Instance()->Get("GRP/GRP/Data",run);
363   AliGRPObject* grpData = 0x0;
364    
365   if(entry){
366     TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
367     if(m){
368        m->Print();
369        grpData = new AliGRPObject();
370        grpData->ReadValuesFromMap(m);
371     }
372     else{
373        grpData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
374        entry->SetOwner(0);
375     }
376     AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
377   }
378
379   if(!grpData) AliError("No GRP entry found in OCDB!");
380
381 // Retrieval of collision mode
382
383   TString beamType = grpData->GetBeamType();
384   if(beamType==AliGRPObject::GetInvalidString()){
385      AliError("GRP/GRP/Data entry:  missing value for the beam type !");
386      AliError("\t VZERO cannot retrieve beam type\n");
387      return;
388   }
389
390    if( (beamType.CompareTo("P-P") ==0)  || (beamType.CompareTo("p-p") ==0) ){
391        fCollisionMode=0;
392   }
393    else if( (beamType.CompareTo("Pb-Pb") ==0)  || (beamType.CompareTo("A-A") ==0) ){
394        fCollisionMode=1;
395    }
396     
397   fBeamEnergy = grpData->GetBeamEnergy();
398   if(fBeamEnergy==AliGRPObject::GetInvalidFloat()) {
399      AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
400      fBeamEnergy = 0.;
401   }
402   
403 //     printf("\n ++++++ Beam type and collision mode retrieved as %s %d @ %1.3f GeV ++++++\n\n",beamType.Data(), fCollisionMode, fBeamEnergy);
404
405 }
406
407 //____________________________________________________________________________
408 AliVZEROCalibData* AliVZERODigitizer::GetCalibData() const
409
410 {
411   AliCDBManager *man = AliCDBManager::Instance();
412
413   AliCDBEntry *entry=0;
414
415   entry = man->Get("VZERO/Calib/Data");
416
417 //   if(!entry){
418 //     AliWarning("Load of calibration data from default storage failed!");
419 //     AliWarning("Calibration data will be loaded from local storage ($ALICE_ROOT)");
420 //     Int_t runNumber = man->GetRun();
421 //     entry = man->GetStorage("local://$ALICE_ROOT/OCDB")
422 //       ->Get("VZERO/Calib/Data",runNumber);
423 //      
424 //   }
425
426   // Retrieval of data in directory VZERO/Calib/Data:
427
428
429   AliVZEROCalibData *calibdata = 0;
430
431   if (entry) calibdata = (AliVZEROCalibData*) entry->GetObject();
432   if (!calibdata)  AliFatal("No calibration data from calibration database !");
433
434   return calibdata;
435
436 }
437
438 //____________________________________________________________________________
439 Int_t AliVZERODigitizer::GetPMNumber(Int_t cell) const
440
441 {
442    
443   Int_t pmNumber[80] = { 0,  1,  2,  3,  4,  5,  6,  7,
444                           8,  9, 10, 11, 12, 13, 14, 15, 
445                          16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 
446                          24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31,
447                          32, 33, 34, 35, 36, 37, 38, 39,
448                          40, 41, 42, 43, 44, 45, 46, 47, 
449                          48, 49, 50, 51, 52, 53, 54, 55,
450                          56, 57, 58, 59, 60, 61, 62, 63 };
451                               
452   return pmNumber[cell];        
453 }
454
455 double AliVZERODigitizer::SignalShape(double *x, double *par)
456 {
457   Double_t xx = x[0];
458   if (xx <= par[0]) return 0;
459   Double_t a = 1./TMath::Power((xx-par[0])/par[1],1./par[2]);
460   if (xx <= par[3]) return a;
461   Double_t b = 1./TMath::Power((xx-par[3])/par[4],1./par[5]);
462   //  Double_t f = TMath::Min(a,b);
463   Double_t f = a*b/(a+b);
464   AliDebug(100,Form("x=%f func=%f",xx,f));
465   return f;
466 }