1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///_________________________________________________________________________
20 /// This class constructs Digits out of Hits
24 // --- Standard library ---
26 // --- ROOT system ---
30 #include <TGeoManager.h>
31 #include <TGeoPhysicalNode.h>
32 #include <AliGeomManager.h>
36 // --- AliRoot header files ---
37 #include "AliVZEROConst.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"
50 #include "AliVZEROdigit.h"
51 #include "AliVZERODigitizer.h"
53 ClassImp(AliVZERODigitizer)
55 AliVZERODigitizer::AliVZERODigitizer()
57 fCalibData(GetCalibData()),
58 fPhotoCathodeEfficiency(0.18),
60 fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)),
68 // default constructor
73 // fPhotoCathodeEfficiency = 0.18;
74 // fPMVoltage = 768.0;
75 // fPMGain = TMath::Power((fPMVoltage / 112.5) ,7.04277);
77 // fCalibData = GetCalibData();
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);
83 //____________________________________________________________________________
84 AliVZERODigitizer::AliVZERODigitizer(AliRunDigitizer* manager)
85 :AliDigitizer(manager),
86 fCalibData(GetCalibData()),
87 fPhotoCathodeEfficiency(0.18),
89 fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)),
102 // fPhotoCathodeEfficiency = 0.18;
103 // fPMVoltage = 768.0;
104 // fPMGain = TMath::Power( (fPMVoltage / 112.5) ,7.04277 );
106 // fCalibData = GetCalibData();
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);
113 //____________________________________________________________________________
114 AliVZERODigitizer::~AliVZERODigitizer()
130 //_____________________________________________________________________________
131 Bool_t AliVZERODigitizer::Init()
133 // Initialises the digitizer
135 // Initialises the Digit array
136 fDigits = new TClonesArray ("AliVZEROdigit", 1000);
138 // TGeoHMatrix *im = AliGeomManager::GetMatrix("VZERO/V0C");
145 //____________________________________________________________________________
146 void AliVZERODigitizer::Exec(Option_t* /*option*/)
148 // Creates digits from hits
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];
157 Float_t pmGain_smeared[64];
160 // Smearing of the PM tubes intrinsic characteristics
162 for(Int_t ii=0; ii<64; ii++) {
163 pmGain_smeared[ii] = gRandom->Gaus(fPMGain, fPMGain/5); }
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
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
173 for(Int_t i=0; i<16; i++) {
174 adc_gain[i] = fCalibData->GetGain(i);
175 cPM[i] = fPhotoCathodeEfficiency * pmGain_smeared[i];
178 for(Int_t j=16; j<48; j=j+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];
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];
191 for(Int_t i=0; i<64; i++){
192 adc_pedestal[i] = fCalibData->GetPedestal(i);
193 adc_sigma[i] = fCalibData->GetSigma(i);
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] );}
199 AliRunLoader* outRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
201 Error("Exec", "Can not get output Run Loader");
205 AliLoader* outLoader = outRunLoader->GetLoader("VZEROLoader");
208 Error("Exec", "Can not get output VZERO Loader");
212 const char* mode = "update";
213 if(outRunLoader->GetEventNumber() == 0) mode = "recreate";
214 outLoader->LoadDigits(mode);
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);
222 for (Int_t iInput = 0; iInput < fManager->GetNinputs(); iInput++) {
223 AliRunLoader* runLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
224 AliLoader* loader = runLoader->GetLoader("VZEROLoader");
226 Error("Exec", "Can not get VZERO Loader for input %d", iInput);
230 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
232 AliVZERO* vzero = (AliVZERO*) runLoader->GetAliRun()->GetDetector("VZERO");
234 Error("Exec", "No VZERO detector for input %d", iInput);
239 TTree* treeH = loader->TreeH();
241 Error("Exec", "Cannot get TreeH for input %d", iInput);
245 for(Int_t i=0; i<80; i++) {map[i] = 0; time[i] = 0.0;}
247 TClonesArray* hits = vzero->Hits();
249 // Now makes Digits from hits
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;}
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();
265 if(t < time_ref[cell]) time_ref[cell] = t;
266 time[cell] = TMath::Min(t,time_ref[cell]);
271 loader->UnloadHits();
275 // Now builds the scintillator cell response (80 cells i.e. 80 responses)
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))
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.) );}
291 // Now transforms 80 cell responses into 64 photomultiplier responses
292 // Also adds the ADC pedestals taken out of the calibration data base
294 for (Int_t j=0; j<16; j++){
295 adc[j] = map [j] + gRandom->Gaus(adc_pedestal[j], adc_sigma[j]);
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]; }
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]);}
310 // Now add digits to the digit Tree
312 for (Int_t i=0; i<64; i++) {
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 :
318 AddDigit(i, adc[i], time2[i] ) ;
322 outLoader->WriteDigits("OVERWRITE");
323 outLoader->UnloadDigits();
327 //____________________________________________________________________________
328 void AliVZERODigitizer::AddDigit(Int_t PMnumber, Float_t adc, Float_t time)
333 TClonesArray &ldigits = *fDigits;
335 if (((Int_t) gRandom->Uniform(2))<1) integrator = kFALSE;
336 else integrator = kTRUE;
338 new(ldigits[fNdigits++]) AliVZEROdigit(PMnumber,adc,time,0,kFALSE,kFALSE,integrator);
341 //____________________________________________________________________________
342 void AliVZERODigitizer::ResetDigit()
348 if (fDigits) fDigits->Delete();
351 //____________________________________________________________________________
352 void AliVZERODigitizer::GetCollisionMode()
354 // Retrieves the collision mode from GRP data
356 // Initialization of the GRP entry
358 Int_t run = AliCDBManager::Instance()->GetRun();
360 // printf("\n ++++++ Run Number retrieved as %d \n",run);
362 AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data",run);
363 AliGRPObject* grpData = 0x0;
366 TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry
369 grpData = new AliGRPObject();
370 grpData->ReadValuesFromMap(m);
373 grpData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry
376 AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
379 if(!grpData) AliError("No GRP entry found in OCDB!");
381 // Retrieval of collision mode
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");
390 if( (beamType.CompareTo("P-P") ==0) || (beamType.CompareTo("p-p") ==0) ){
393 else if( (beamType.CompareTo("Pb-Pb") ==0) || (beamType.CompareTo("A-A") ==0) ){
397 fBeamEnergy = grpData->GetBeamEnergy();
398 if(fBeamEnergy==AliGRPObject::GetInvalidFloat()) {
399 AliError("GRP/GRP/Data entry: missing value for the beam energy ! Using 0");
403 // printf("\n ++++++ Beam type and collision mode retrieved as %s %d @ %1.3f GeV ++++++\n\n",beamType.Data(), fCollisionMode, fBeamEnergy);
407 //____________________________________________________________________________
408 AliVZEROCalibData* AliVZERODigitizer::GetCalibData() const
411 AliCDBManager *man = AliCDBManager::Instance();
413 AliCDBEntry *entry=0;
415 entry = man->Get("VZERO/Calib/Data");
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);
426 // Retrieval of data in directory VZERO/Calib/Data:
429 AliVZEROCalibData *calibdata = 0;
431 if (entry) calibdata = (AliVZEROCalibData*) entry->GetObject();
432 if (!calibdata) AliFatal("No calibration data from calibration database !");
438 //____________________________________________________________________________
439 Int_t AliVZERODigitizer::GetPMNumber(Int_t cell) const
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 };
452 return pmNumber[cell];
455 double AliVZERODigitizer::SignalShape(double *x, double *par)
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));