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