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