Macros used to produce the OCDB entries for MC. One for default (p-p) and one for...
[u/mrichter/AliRoot.git] / VZERO / AliVZERODigitizer.cxx
CommitLineData
9e04c3b6 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
090026bf 16/* $Id$ */
17
b0d2c2d3 18///_________________________________________________________________________
19///
20/// This class constructs Digits out of Hits
21///
22///
9e04c3b6 23
24// --- Standard library ---
9e04c3b6 25
26// --- ROOT system ---
090026bf 27#include <TMath.h>
9e04c3b6 28#include <TTree.h>
fe0adf2a 29#include <TMap.h>
8bd3e27a 30#include <TGeoManager.h>
31#include <TGeoPhysicalNode.h>
32#include <AliGeomManager.h>
e939a978 33#include <TRandom.h>
9672d66e 34#include <TF1.h>
9e04c3b6 35
36// --- AliRoot header files ---
37#include "AliVZEROConst.h"
38#include "AliRun.h"
39#include "AliVZERO.h"
40#include "AliVZEROhit.h"
9e04c3b6 41#include "AliRunLoader.h"
42#include "AliLoader.h"
fe0adf2a 43#include "AliGRPObject.h"
9e04c3b6 44#include "AliRunDigitizer.h"
ce7090f5 45#include "AliCDBManager.h"
46#include "AliCDBStorage.h"
47#include "AliCDBEntry.h"
48#include "AliVZEROCalibData.h"
49
9e04c3b6 50#include "AliVZEROdigit.h"
b0d2c2d3 51#include "AliVZERODigitizer.h"
9e04c3b6 52
53ClassImp(AliVZERODigitizer)
54
55 AliVZERODigitizer::AliVZERODigitizer()
fe0adf2a 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),
9672d66e 64 fBeamEnergy(0.),
65 fSignalShape(NULL)
0b2bea8b 66
9e04c3b6 67{
b0d2c2d3 68 // default constructor
69
0b2bea8b 70// fNdigits = 0;
71// fDigits = 0;
72//
73// fPhotoCathodeEfficiency = 0.18;
74// fPMVoltage = 768.0;
75// fPMGain = TMath::Power((fPMVoltage / 112.5) ,7.04277);
ce7090f5 76
0b2bea8b 77// fCalibData = GetCalibData();
9672d66e 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);
9e04c3b6 81}
82
83//____________________________________________________________________________
84 AliVZERODigitizer::AliVZERODigitizer(AliRunDigitizer* manager)
0b2bea8b 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),
fe0adf2a 91 fDigits(0),
92 fCollisionMode(0),
9672d66e 93 fBeamEnergy(0.),
94 fSignalShape(NULL)
0b2bea8b 95
9e04c3b6 96{
97 // constructor
98
0b2bea8b 99// fNdigits = 0;
100// fDigits = 0;
101//
102// fPhotoCathodeEfficiency = 0.18;
103// fPMVoltage = 768.0;
104// fPMGain = TMath::Power( (fPMVoltage / 112.5) ,7.04277 );
ce7090f5 105
0b2bea8b 106// fCalibData = GetCalibData();
9672d66e 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);
ce7090f5 110
9e04c3b6 111}
112
113//____________________________________________________________________________
114 AliVZERODigitizer::~AliVZERODigitizer()
115{
116 // destructor
117
118 if (fDigits) {
119 fDigits->Delete();
120 delete fDigits;
b0d2c2d3 121 fDigits=0;
122 }
9672d66e 123
124 if (fSignalShape) {
125 delete fSignalShape;
126 fSignalShape = NULL;
127 }
9e04c3b6 128}
129
b0d2c2d3 130//_____________________________________________________________________________
131Bool_t AliVZERODigitizer::Init()
9e04c3b6 132{
b0d2c2d3 133 // Initialises the digitizer
9e04c3b6 134
b0d2c2d3 135 // Initialises the Digit array
9e04c3b6 136 fDigits = new TClonesArray ("AliVZEROdigit", 1000);
8bd3e27a 137
c61a7285 138 // TGeoHMatrix *im = AliGeomManager::GetMatrix("VZERO/V0C");
139 // im->Print();
fe0adf2a 140
141 GetCollisionMode();
b0d2c2d3 142 return kTRUE;
9e04c3b6 143}
144
145//____________________________________________________________________________
b0d2c2d3 146void AliVZERODigitizer::Exec(Option_t* /*option*/)
ce7090f5 147{
b0d2c2d3 148 // Creates digits from hits
ce7090f5 149
db0db003 150 Float_t map[80]; // 48 values on V0C + 32 on V0A
f25f4990 151// Int_t pmNumber[80];
db0db003 152 Float_t adc[64]; // 32 PMs on V0C + 32 PMs on V0A
8bd3e27a 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];
ce7090f5 156 fNdigits = 0;
8adc9b44 157 Float_t pmGain_smeared[64];
8bd3e27a 158 Float_t cPM[80];
159
160 // Smearing of the PM tubes intrinsic characteristics
161
162 for(Int_t ii=0; ii<64; ii++) {
8adc9b44 163 pmGain_smeared[ii] = gRandom->Gaus(fPMGain, fPMGain/5); }
8bd3e27a 164
165 // Retrieval of ADC gain values and pedestal information from local CDB
ce7090f5 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
b0d2c2d3 168
ce7090f5 169 // Reminder : We have 16 scintillating cells mounted on 8 PMs
20277079 170 // on Ring 3 and Ring 4 in V0C - added to produce ADC outputs
171 // on these rings...
ce7090f5 172
8bd3e27a 173 for(Int_t i=0; i<16; i++) {
fad64858 174 adc_gain[i] = fCalibData->GetGain(i);
175 cPM[i] = fPhotoCathodeEfficiency * pmGain_smeared[i];
176 }
ce7090f5 177
178 for(Int_t j=16; j<48; j=j+2) {
fad64858 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 }
8bd3e27a 185
186 for(Int_t i=48; i<80; i++){
db0db003 187 adc_gain[i] = fCalibData->GetGain(i-16);
fad64858 188 cPM[i] = fPhotoCathodeEfficiency * pmGain_smeared[i-16];
189 };
ce7090f5 190
fad64858 191 for(Int_t i=0; i<64; i++){
192 adc_pedestal[i] = fCalibData->GetPedestal(i);
193 adc_sigma[i] = fCalibData->GetSigma(i);
194 };
8bd3e27a 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] );}
ce7090f5 198
fad64858 199 AliRunLoader* outRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
b0d2c2d3 200 if (!outRunLoader) {
201 Error("Exec", "Can not get output Run Loader");
fad64858 202 return;
203 }
ce7090f5 204
b0d2c2d3 205 AliLoader* outLoader = outRunLoader->GetLoader("VZEROLoader");
fad64858 206
b0d2c2d3 207 if (!outLoader) {
208 Error("Exec", "Can not get output VZERO Loader");
fad64858 209 return;
210 }
b0d2c2d3 211
e539620f 212 const char* mode = "update";
213 if(outRunLoader->GetEventNumber() == 0) mode = "recreate";
214 outLoader->LoadDigits(mode);
fad64858 215
b0d2c2d3 216 if (!outLoader->TreeD()) outLoader->MakeTree("D");
217 outLoader->MakeDigitsContainer();
ce7090f5 218 TTree* treeD = outLoader->TreeD();
b0d2c2d3 219 Int_t bufsize = 16000;
220 treeD->Branch("VZERODigit", &fDigits, bufsize);
c299dbe4 221
b0d2c2d3 222 for (Int_t iInput = 0; iInput < fManager->GetNinputs(); iInput++) {
fad64858 223 AliRunLoader* runLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
8bd3e27a 224 AliLoader* loader = runLoader->GetLoader("VZEROLoader");
225 if (!loader) {
226 Error("Exec", "Can not get VZERO Loader for input %d", iInput);
fad64858 227 continue;
228 }
ce7090f5 229
8bd3e27a 230 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
b0d2c2d3 231
8bd3e27a 232 AliVZERO* vzero = (AliVZERO*) runLoader->GetAliRun()->GetDetector("VZERO");
233 if (!vzero) {
234 Error("Exec", "No VZERO detector for input %d", iInput);
fad64858 235 continue;
236 }
9e04c3b6 237
8bd3e27a 238 loader->LoadHits();
239 TTree* treeH = loader->TreeH();
240 if (!treeH) {
241 Error("Exec", "Cannot get TreeH for input %d", iInput);
fad64858 242 continue;
243 }
c299dbe4 244
8bd3e27a 245 for(Int_t i=0; i<80; i++) {map[i] = 0; time[i] = 0.0;}
fe0adf2a 246
8bd3e27a 247 TClonesArray* hits = vzero->Hits();
9e04c3b6 248
b0d2c2d3 249// Now makes Digits from hits
fe0adf2a 250
8bd3e27a 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++) {
fad64858 258 AliVZEROhit* hit = (AliVZEROhit *)hits->UncheckedAt(iHit);
259 Int_t nPhot = hit->Nphot();
260 Int_t cell = hit->Cell();
db0db003 261 map[cell] += Float_t(nPhot);
fad64858 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 }
8bd3e27a 268 } // hit loop
269 } // track loop
b0d2c2d3 270
8bd3e27a 271 loader->UnloadHits();
b0d2c2d3 272
273 } // input loop
20277079 274
275// Now builds the scintillator cell response (80 cells i.e. 80 responses)
841137ce 276
20277079 277 for (Int_t i=0; i<80; i++) {
8adc9b44 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))
8bd3e27a 281 + noise*1e-3;
fe0adf2a 282 if(fCollisionMode >0) adc_gain[i] = adc_gain[i]/70.0; // reduce dynamics in Ion Collision Mode
db0db003 283 map[i] = pmResponse * adc_gain[i];
fe0adf2a 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 );
db0db003 287 if(map[i] > (MIP/2.) )
288 {map[i] = gRandom->Gaus(map[i], (MIP/6.) );}
8adc9b44 289 }
c299dbe4 290
8bd3e27a 291// Now transforms 80 cell responses into 64 photomultiplier responses
292// Also adds the ADC pedestals taken out of the calibration data base
20277079 293
8bd3e27a 294 for (Int_t j=0; j<16; j++){
db0db003 295 adc[j] = map [j] + gRandom->Gaus(adc_pedestal[j], adc_sigma[j]);
20277079 296 time2[j]= time[j];}
297
298 for (Int_t j=48; j<80; j++){
db0db003 299 adc[j-16] = map [j] + gRandom->Gaus(adc_pedestal[j-16],adc_sigma[j-16]);
8bd3e27a 300 time2[j-16]= time[j]; }
20277079 301
302 for (Int_t j=0; j<16; j++){
db0db003 303 adc[16+j] = map [16+2*j]+ map [16+2*j+1] + gRandom->Gaus(adc_pedestal[16+j], adc_sigma[16+j]);
c299dbe4 304 Float_t min_time = TMath::Min(time [16+2*j],time [16+2*j+1]);
305 time2[16+j] = min_time;
8bd3e27a 306 if(min_time==0.0){time2[16+j]=TMath::Max(time[16+2*j],time[16+2*j+1]);}
c299dbe4 307 }
308
309
20277079 310// Now add digits to the digit Tree
311
312 for (Int_t i=0; i<64; i++) {
313 if(adc[i] > 0) {
fad64858 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));
8bd3e27a 316// multiply by 10 to have 100 ps per channel :
fad64858 317
4682ffee 318 AddDigit(i, adc[i], time2[i] ) ;
fad64858 319 }
20277079 320 }
b0d2c2d3 321 treeD->Fill();
322 outLoader->WriteDigits("OVERWRITE");
323 outLoader->UnloadDigits();
324 ResetDigit();
9e04c3b6 325}
326
327//____________________________________________________________________________
db0db003 328void AliVZERODigitizer::AddDigit(Int_t PMnumber, Float_t adc, Float_t time)
9e04c3b6 329 {
330
331// Adds Digit
332
333 TClonesArray &ldigits = *fDigits;
fad64858 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
9e04c3b6 340}
341//____________________________________________________________________________
342void AliVZERODigitizer::ResetDigit()
343{
fe0adf2a 344
9e04c3b6 345// Clears Digits
fe0adf2a 346
9e04c3b6 347 fNdigits = 0;
b0d2c2d3 348 if (fDigits) fDigits->Delete();
9e04c3b6 349}
ce7090f5 350
351//____________________________________________________________________________
fe0adf2a 352void 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//____________________________________________________________________________
ce7090f5 408AliVZEROCalibData* AliVZERODigitizer::GetCalibData() const
409
410{
c0b82b5a 411 AliCDBManager *man = AliCDBManager::Instance();
ce7090f5 412
c0b82b5a 413 AliCDBEntry *entry=0;
ce7090f5 414
c0b82b5a 415 entry = man->Get("VZERO/Calib/Data");
416
df791951 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();
162637e4 421// entry = man->GetStorage("local://$ALICE_ROOT/OCDB")
df791951 422// ->Get("VZERO/Calib/Data",runNumber);
423//
424// }
c0b82b5a 425
426 // Retrieval of data in directory VZERO/Calib/Data:
ce7090f5 427
ce7090f5 428
c0b82b5a 429 AliVZEROCalibData *calibdata = 0;
ce7090f5 430
c0b82b5a 431 if (entry) calibdata = (AliVZEROCalibData*) entry->GetObject();
df791951 432 if (!calibdata) AliFatal("No calibration data from calibration database !");
ce7090f5 433
c0b82b5a 434 return calibdata;
ce7090f5 435
436}
437
8adc9b44 438//____________________________________________________________________________
439Int_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
9672d66e 455double 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}