]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ZDC/AliZDCDigitizer.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
MC calculation of Q vector for AMPT
[u/mrichter/AliRoot.git] / ZDC / AliZDCDigitizer.cxx
CommitLineData
8309c1ab 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// ZDC digitizer class //
21// //
22///////////////////////////////////////////////////////////////////////////////
23
3d7e0d0c 24// --- Standard libraries
8309c1ab 25#include <stdlib.h>
3d7e0d0c 26#include <stdio.h>
8309c1ab 27
28// --- ROOT system
29#include <TTree.h>
30#include <TFile.h>
31#include <TNtuple.h>
32#include <TRandom.h>
33
34// --- AliRoot header files
8309c1ab 35#include "AliRun.h"
36#include "AliHeader.h"
37#include "AliGenHijingEventHeader.h"
07e38cef 38#include "AliGenCocktailEventHeader.h"
f21fc003 39#include "AliDigitizationInput.h"
8309c1ab 40#include "AliRunLoader.h"
c93255fe 41#include "AliLoader.h"
12434381 42#include "AliGRPObject.h"
78d18275 43#include "AliCDBManager.h"
48642b09 44#include "AliCDBEntry.h"
8309c1ab 45#include "AliZDCSDigit.h"
46#include "AliZDCDigit.h"
47#include "AliZDCFragment.h"
a18a0543 48#include "AliZDCv3.h"
8309c1ab 49#include "AliZDCDigitizer.h"
8a2624cc 50
51class AliCDBStorage;
6024ec85 52class AliZDCPedestals;
8309c1ab 53
54ClassImp(AliZDCDigitizer)
55
56
57//____________________________________________________________________________
a718c993 58AliZDCDigitizer::AliZDCDigitizer() :
59 fIsCalibration(0),
fd9afd60 60 fIsSignalInADCGate(kFALSE),
61 fFracLostSignal(0.),
a718c993 62 fPedData(0),
2d9c70ab 63 fSpectators2Track(kFALSE),
f2e2aa97 64 fBeamEnergy(0.),
5d24ec95 65 fBeamType(""),
66 fIspASystem(kFALSE),
94ff33bf 67 fIsRELDISgen(kFALSE),
68 fSpectatorData(0)
8309c1ab 69{
698e394d 70 // Default constructor
83534754 71 for(Int_t i=0; i<2; i++) fADCRes[i]=0.;
8309c1ab 72}
73
74//____________________________________________________________________________
f21fc003 75AliZDCDigitizer::AliZDCDigitizer(AliDigitizationInput* digInput):
76 AliDigitizer(digInput),
a718c993 77 fIsCalibration(0), //By default the simulation doesn't create calib. data
fd9afd60 78 fIsSignalInADCGate(kFALSE),
79 fFracLostSignal(0.),
a718c993 80 fPedData(GetPedData()),
2d9c70ab 81 fSpectators2Track(kFALSE),
f2e2aa97 82 fBeamEnergy(0.),
5d24ec95 83 fBeamType(""),
84 fIspASystem(kFALSE),
94ff33bf 85 fIsRELDISgen(kFALSE),
86 fSpectatorData(NULL)
8309c1ab 87{
78d18275 88 // Get calibration data
d27e20dd 89 if(fIsCalibration!=0) printf("\n\t AliZDCDigitizer -> Creating calibration data (pedestals)\n");
edb4e893 90 for(Int_t i=0; i<5; i++){
698e394d 91 for(Int_t j=0; j<5; j++)
92 fPMGain[i][j] = 0.;
93 }
8309c1ab 94}
95
96//____________________________________________________________________________
97AliZDCDigitizer::~AliZDCDigitizer()
98{
99// Destructor
94ff33bf 100 if(fSpectatorData!=NULL){
101 fSpectatorData->Close();
102 delete fSpectatorData;
103 }
8309c1ab 104}
105
106
cc2abffd 107//____________________________________________________________________________
108AliZDCDigitizer::AliZDCDigitizer(const AliZDCDigitizer &digitizer):
a718c993 109 AliDigitizer(),
110 fIsCalibration(digitizer.fIsCalibration),
fd9afd60 111 fIsSignalInADCGate(digitizer.fIsSignalInADCGate),
112 fFracLostSignal(digitizer.fFracLostSignal),
a718c993 113 fPedData(digitizer.fPedData),
2d9c70ab 114 fSpectators2Track(digitizer.fSpectators2Track),
f2e2aa97 115 fBeamEnergy(digitizer.fBeamEnergy),
5d24ec95 116 fBeamType(digitizer.fBeamType),
117 fIspASystem(digitizer.fIspASystem),
94ff33bf 118 fIsRELDISgen(digitizer.fIsRELDISgen),
119 fSpectatorData(digitizer.fSpectatorData)
cc2abffd 120{
121 // Copy constructor
122
edb4e893 123 for(Int_t i=0; i<5; i++){
cc2abffd 124 for(Int_t j=0; j<5; j++){
125 fPMGain[i][j] = digitizer.fPMGain[i][j];
126 }
127 }
128 for(Int_t i=0; i<2; i++) fADCRes[i] = digitizer.fADCRes[i];
cc2abffd 129
12434381 130
cc2abffd 131}
132
8309c1ab 133//____________________________________________________________________________
134Bool_t AliZDCDigitizer::Init()
135{
48642b09 136 // Initialize the digitizer
12434381 137
5d24ec95 138
139 //printf(" **** Initializing AliZDCDigitizer fBeamEnergy = %1.0f GeV\n\n", fBeamEnergy);
12434381 140 AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
83534754 141 if(!entry) AliFatal("No calibration data loaded!");
12434381 142 AliGRPObject* grpData = 0x0;
143 if(entry){
144 TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry
145 if(m){
146 //m->Print();
147 grpData = new AliGRPObject();
148 grpData->ReadValuesFromMap(m);
149 }
150 else{
151 grpData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry
152 }
153 entry->SetOwner(0);
154 AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
155 }
83534754 156 if(!grpData){
157 AliError("No GRP entry found in OCDB! \n ");
158 return kFALSE;
159 }
12434381 160
5d24ec95 161 fBeamType = grpData->GetBeamType();
162 if(fBeamType==AliGRPObject::GetInvalidString()){
12434381 163 AliError("\t UNKNOWN beam type from GRP obj -> PMT gains not set in ZDC digitizer!!!\n");
164 }
165
5d24ec95 166 // If beam energy is not set from Config.C (RELDIS / spectator generators)
167 if(fBeamEnergy<0.01){
168 fBeamEnergy = grpData->GetBeamEnergy();
169 if(fBeamEnergy==AliGRPObject::GetInvalidFloat()){
170 AliWarning("GRP/GRP/Data entry: missing value for the beam energy ! Using 0.");
171 AliError("\t UNKNOWN beam type from GRP obj -> PMT gains not set in ZDC digitizer!!!\n");
172 fBeamEnergy = 0.;
a736fd93 173 }
12434381 174 }
f2e2aa97 175
f6c61b40 176 /*if(fIspASystem){
5d24ec95 177 fBeamType = "p-A";
178 AliInfo(" AliZDCDigitizer -> Manually setting beam type to p-A\n");
e40b9288 179 }*/
5d24ec95 180
181 // Setting beam type for spectator generator and RELDIS generator
182 if(((fBeamType.CompareTo("UNKNOWN")) == 0) || fIsRELDISgen){
183 fBeamType = "A-A";
184 AliInfo(" AliZDCDigitizer -> Manually setting beam type to A-A\n");
185 }
8f190cf4 186 printf("\n\t AliZDCDigitizer -> beam type %s - beam energy = %f GeV\n", fBeamType.Data(), fBeamEnergy);
94ff33bf 187 if(fSpectators2Track) printf("\n\t AliZDCDigitizer -> spectator signal added at digit level\n\n");
5d24ec95 188
db247fc3 189 if(fBeamEnergy>0.1){
190 ReadPMTGains();
f6c61b40 191 //CalculatePMTGains();
db247fc3 192 }
193 else{
194 AliWarning("\n Beam energy is 0 -> ZDC PMT gains can't be set -> NO ZDC DIGITS!!!\n");
195 }
196
698e394d 197 // ADC Caen V965
7f73eb6b 198 fADCRes[0] = 0.0000008; // ADC Resolution high gain: 200 fC/adcCh
199 fADCRes[1] = 0.0000064; // ADC Resolution low gain: 25 fC/adcCh
8309c1ab 200
201 return kTRUE;
202}
203
204//____________________________________________________________________________
f21fc003 205void AliZDCDigitizer::Digitize(Option_t* /*option*/)
8309c1ab 206{
48642b09 207 // Execute digitization
8309c1ab 208
d79f8d50 209 // ------------------------------------------------------------
f91d1e35 210 // !!! 2nd ZDC set added
211 // *** 1st 3 arrays are digits from REAL (simulated) hits
212 // *** last 2 are copied from simulated digits
fd9afd60 213 // --- pm[0][...] = light in ZN side C [C, Q1, Q2, Q3, Q4]
214 // --- pm[1][...] = light in ZP side C [C, Q1, Q2, Q3, Q4]
48642b09 215 // --- pm[2][...] = light in ZEM [x, 1, 2, x, x]
fd9afd60 216 // --- pm[3][...] = light in ZN side A [C, Q1, Q2, Q3, Q4] ->NEW!
217 // --- pm[4][...] = light in ZP side A [C, Q1, Q2, Q3, Q4] ->NEW!
d79f8d50 218 // ------------------------------------------------------------
219 Float_t pm[5][5];
8574b308 220 for(Int_t iSector1=0; iSector1<5; iSector1++)
221 for(Int_t iSector2=0; iSector2<5; iSector2++){
8309c1ab 222 pm[iSector1][iSector2] = 0;
223 }
d79f8d50 224
225 // ------------------------------------------------------------
226 // ### Out of time ADC added (22 channels)
227 // --- same codification as for signal PTMs (see above)
228 // ------------------------------------------------------------
3d7ab74d 229 // Float_t pmoot[5][5];
230 // for(Int_t iSector1=0; iSector1<5; iSector1++)
231 // for(Int_t iSector2=0; iSector2<5; iSector2++){
232 // pmoot[iSector1][iSector2] = 0;
233 // }
8309c1ab 234
235 // impact parameter and number of spectators
131b5d31 236 Float_t impPar = 0;
fd9afd60 237 Int_t specNTarg = 0, specPTarg = 0;
238 Int_t specNProj = 0, specPProj = 0;
a18a0543 239 Float_t signalTime0 = 0.;
8309c1ab 240
241 // loop over input streams
f21fc003 242 for(Int_t iInput = 0; iInput<fDigInput->GetNinputs(); iInput++){
8309c1ab 243
244 // get run loader and ZDC loader
245 AliRunLoader* runLoader =
f21fc003 246 AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(iInput));
8309c1ab 247 AliLoader* loader = runLoader->GetLoader("ZDCLoader");
8574b308 248 if(!loader) continue;
8309c1ab 249
250 // load sdigits
251 loader->LoadSDigits();
252 TTree* treeS = loader->TreeS();
8574b308 253 if(!treeS) continue;
8309c1ab 254 AliZDCSDigit sdigit;
255 AliZDCSDigit* psdigit = &sdigit;
256 treeS->SetBranchAddress("ZDC", &psdigit);
257
258 // loop over sdigits
8574b308 259 for(Int_t iSDigit=0; iSDigit<treeS->GetEntries(); iSDigit++){
8309c1ab 260 treeS->GetEntry(iSDigit);
7f73eb6b 261 //
8574b308 262 if(!psdigit) continue;
263 if((sdigit.GetSector(1) < 0) || (sdigit.GetSector(1) > 4)){
8309c1ab 264 AliError(Form("\nsector[0] = %d, sector[1] = %d\n",
265 sdigit.GetSector(0), sdigit.GetSector(1)));
266 continue;
267 }
fd9afd60 268 // Checking if signal is inside ADC gate
269 if(iSDigit==0) signalTime0 = sdigit.GetTrackTime();
270 else{
271 // Assuming a signal lenght of 20 ns, signal is in gate if
272 // signal ENDS in signalTime0+50. -> sdigit.GetTrackTime()+20<=signalTime0+50.
273 if(sdigit.GetTrackTime()<=signalTime0+30.) fIsSignalInADCGate = kTRUE;
274 if(sdigit.GetTrackTime()>signalTime0+30.){
275 fIsSignalInADCGate = kFALSE;
5d24ec95 276 // Vedi quaderno per spiegazione approx. usata nel calcolo della fraz. di segnale perso
fd9afd60 277 fFracLostSignal = (sdigit.GetTrackTime()-30)*(sdigit.GetTrackTime()-30)/280.;
278 }
279 }
698e394d 280
fd9afd60 281 Float_t sdSignal = sdigit.GetLightPM();
282 if(fIsSignalInADCGate == kFALSE){
7e858b1d 283 AliDebug(2,Form("\t Signal time %f -> fraction %f of ZDC signal for det.(%d, %d) out of ADC gate\n",
0fba8107 284 sdigit.GetTrackTime(),fFracLostSignal,sdigit.GetSector(0),sdigit.GetSector(1)));
fd9afd60 285 sdSignal = (1-fFracLostSignal)*sdSignal;
286 }
287
48642b09 288 pm[(sdigit.GetSector(0))-1][sdigit.GetSector(1)] += sdigit.GetLightPM();
698e394d 289 //Ch. debug
290 /*printf("\t Detector %d, Tower %d -> pm[%d][%d] = %.0f \n",
48642b09 291 sdigit.GetSector(0), sdigit.GetSector(1),sdigit.GetSector(0)-1,
698e394d 292 sdigit.GetSector(1), pm[sdigit.GetSector(0)-1][sdigit.GetSector(1)]);
878fa0b3 293 */
fd9afd60 294
8309c1ab 295 }
296
8309c1ab 297 loader->UnloadSDigits();
298
299 // get the impact parameter and the number of spectators in case of hijing
8574b308 300 if(!runLoader->GetAliRun()) runLoader->LoadgAlice();
d3b3a3b2 301 AliHeader* header = runLoader->GetHeader();
8574b308 302 if(!header) continue;
8309c1ab 303 AliGenEventHeader* genHeader = header->GenEventHeader();
8574b308 304 if(!genHeader) continue;
07e38cef 305 AliGenHijingEventHeader *hijingHeader = 0;
306 if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())) hijingHeader = dynamic_cast <AliGenHijingEventHeader*> (genHeader);
307 else if(genHeader->InheritsFrom(AliGenCocktailEventHeader::Class())){
308 TList* listOfHeaders = ((AliGenCocktailEventHeader*) genHeader)->GetHeaders();
c47d6938 309 if(listOfHeaders){
4e97feaf 310 for(Int_t iH = 0; iH < listOfHeaders->GetEntries(); ++iH) {
311 AliGenEventHeader *currHeader = dynamic_cast <AliGenEventHeader *> (listOfHeaders->At(iH));
312 if (currHeader && currHeader->InheritsFrom(AliGenHijingEventHeader::Class())) {
313 hijingHeader = dynamic_cast <AliGenHijingEventHeader*> (currHeader);
314 break;
315 }
316 }
c47d6938 317 }
01e01248 318 else{
8f190cf4 319 printf(" No list of headers from generator \n");
01e01248 320 }
07e38cef 321 }
c47d6938 322 if(!hijingHeader){
c47d6938 323 printf(" No HIJING header found in list of headers from generator\n");
c47d6938 324 }
07e38cef 325
a20bc87c 326 if(hijingHeader && fSpectators2Track==kTRUE){
8fced614 327 impPar = hijingHeader->ImpactParameter();
07e38cef 328 specNProj = hijingHeader->ProjSpectatorsn();
329 specPProj = hijingHeader->ProjSpectatorsp();
330 specNTarg = hijingHeader->TargSpectatorsn();
331 specPTarg = hijingHeader->TargSpectatorsp();
94ff33bf 332 /*printf("\t AliZDCDigitizer: b = %1.2f fm\n"
8fced614 333 " \t PROJECTILE: #spectator n %d, #spectator p %d\n"
334 " \t TARGET: #spectator n %d, #spectator p %d\n",
94ff33bf 335 impPar, specNProj, specPProj, specNTarg, specPTarg);*/
a18a0543 336
5d24ec95 337 // Applying fragmentation algorithm and adding spectator signal
338 Int_t freeSpecNProj=0, freeSpecPProj=0;
339 if(specNProj!=0 || specPProj!=0) Fragmentation(impPar, specNProj, specPProj, freeSpecNProj, freeSpecPProj);
340 Int_t freeSpecNTarg=0, freeSpecPTarg=0;
341 if(specNTarg!=0 || specPTarg!=0) Fragmentation(impPar, specNTarg, specPTarg, freeSpecNTarg, freeSpecPTarg);
342 if(freeSpecNProj!=0) SpectatorSignal(1, freeSpecNProj, pm);
343 if(freeSpecPProj!=0) SpectatorSignal(2, freeSpecPProj, pm);
94ff33bf 344 //printf("\t AliZDCDigitizer -> Adding spectator signal for PROJECTILE: %d free n and %d free p\n",freeSpecNProj,freeSpecPProj);
5d24ec95 345 if(freeSpecNTarg!=0) SpectatorSignal(3, freeSpecNTarg, pm);
346 if(freeSpecPTarg!=0) SpectatorSignal(4, freeSpecPTarg, pm);
94ff33bf 347 //printf("\t AliZDCDigitizer -> Adding spectator signal for TARGET: %d free n and %d free p\n",freeSpecNTarg,freeSpecPTarg);
5d24ec95 348 }
8309c1ab 349 }
350
8a2624cc 351
8309c1ab 352 // get the output run loader and loader
353 AliRunLoader* runLoader =
f21fc003 354 AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
8309c1ab 355 AliLoader* loader = runLoader->GetLoader("ZDCLoader");
8574b308 356 if(!loader) {
8309c1ab 357 AliError("no ZDC loader found");
358 return;
359 }
360
361 // create the output tree
362 const char* mode = "update";
8574b308 363 if(runLoader->GetEventNumber() == 0) mode = "recreate";
8309c1ab 364 loader->LoadDigits(mode);
365 loader->MakeTree("D");
366 TTree* treeD = loader->TreeD();
367 AliZDCDigit digit;
368 AliZDCDigit* pdigit = &digit;
369 const Int_t kBufferSize = 4000;
370 treeD->Branch("ZDC", "AliZDCDigit", &pdigit, kBufferSize);
371
372 // Create digits
80e87581 373 Int_t sector[2];
374 Int_t digi[2], digioot[2];
375 for(sector[0]=1; sector[0]<6; sector[0]++){
abf60186 376 for(sector[1]=0; sector[1]<5; sector[1]++){
377 if((sector[0]==3) && ((sector[1]<1) || (sector[1]>2))) continue;
8574b308 378 for(Int_t res=0; res<2; res++){
abf60186 379 digi[res] = Phe2ADCch(sector[0], sector[1], pm[sector[0]-1][sector[1]], res)
48642b09 380 + Pedestal(sector[0], sector[1], res);
7f73eb6b 381 }
abf60186 382 new(pdigit) AliZDCDigit(sector, digi);
8309c1ab 383 treeD->Fill();
698e394d 384
385 //Ch. debug
386 //printf("\t DIGIT added -> det %d quad %d - digi[0,1] = [%d, %d]\n",
e40b9288 387 // sector[0], sector[1], digi[0], digi[1]); // Chiara debugging!
698e394d 388
8309c1ab 389 }
83347831 390 } // Loop over detector
391 // Adding in-time digits for 2 reference PTM signals (after signal ch.)
392 // (for the moment the ref. signal is completely invented assuming a PMgain of 5*10^4!)
393 Int_t sectorRef[2];
394 sectorRef[1] = 5;
9d38ced8 395 Int_t sigRef[2];
80e87581 396 // Reference signal are set to 100 (high gain chain) and 800 (low gain chain)
9d38ced8 397 if(fIsCalibration==0) {sigRef[0]=100; sigRef[1]=800;}
80e87581 398 else {sigRef[0]=0; sigRef[1]=0;} // calibration -> simulation of pedestal values
9d38ced8 399 //
83347831 400 for(Int_t iref=0; iref<2; iref++){
401 sectorRef[0] = 3*iref+1;
402 for(Int_t res=0; res<2; res++){
403 sigRef[res] += Pedestal(sectorRef[0], sectorRef[1], res);
404 }
83347831 405 new(pdigit) AliZDCDigit(sectorRef, sigRef);
406 treeD->Fill();
698e394d 407
408 //Ch. debug
409 //printf("\t RefDigit added -> det = %d, quad = %d - digi[0,1] = [%d, %d]\n",
410 // sectorRef[0], sectorRef[1], sigRef[0], sigRef[1]); // Chiara debugging!
abf60186 411 }
d79f8d50 412 //
413 // --- Adding digits for out-of-time channels after signal digits
80e87581 414 for(sector[0]=1; sector[0]<6; sector[0]++){
d79f8d50 415 for(sector[1]=0; sector[1]<5; sector[1]++){
416 if((sector[0]==3) && ((sector[1]<1) || (sector[1]>2))) continue;
8574b308 417 for(Int_t res=0; res<2; res++){
d79f8d50 418 digioot[res] = Pedestal(sector[0], sector[1], res); // out-of-time ADCs
419 }
d79f8d50 420 new(pdigit) AliZDCDigit(sector, digioot);
421 treeD->Fill();
698e394d 422
423 //Ch. debug
424 //printf("\t DIGIToot added -> det = %d, quad = %d - digi[0,1] = [%d, %d]\n",
425 // sector[0], sector[1], digioot[0], digioot[1]); // Chiara debugging!
d79f8d50 426 }
427 }
83347831 428 // Adding out-of-time digits for 2 reference PTM signals (after out-of-time ch.)
429 Int_t sigRefoot[2];
430 for(Int_t iref=0; iref<2; iref++){
431 sectorRef[0] = 3*iref+1;
432 for(Int_t res=0; res<2; res++){
433 sigRefoot[res] = Pedestal(sectorRef[0], sectorRef[1], res);
434 }
83347831 435 new(pdigit) AliZDCDigit(sectorRef, sigRefoot);
436 treeD->Fill();
698e394d 437 //Ch. debug
438 //printf("\t RefDigitoot added -> det = %d, quad = %d - digi[0,1] = [%d, %d]\n",
439 // sectorRef[0], sectorRef[1], sigRefoot[0], sigRefoot[1]); // Chiara debugging!
440
83347831 441 }
442 //printf("\t AliZDCDigitizer -> TreeD has %d entries\n",(Int_t) treeD->GetEntries());
443
8309c1ab 444 // write the output tree
445 loader->WriteDigits("OVERWRITE");
446 loader->UnloadDigits();
447}
448
449
25175f38 450//_____________________________________________________________________________
451void AliZDCDigitizer::ReadPMTGains()
452{
453// Read PMT gain from an external file
454
455 char *fname = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/PMTGainsdata.txt");
456 FILE *fdata = fopen(fname,"r");
457 if(fdata==NULL){
458 AliWarning(" Can't open file $ALICE_ROOT/ZDC/PMTGainsdata.txt to read ZDC PMT Gains\n");
459 AliWarning(" -> ZDC signal will be pedestal!!!!!!!!!!!!\n\n");
460 return;
461 }
3d7e0d0c 462 int read=1;
25175f38 463 Float_t data[5];
464 Int_t beam[12], det[12];
465 Float_t gain[12], aEne[12], bEne[12];
466 for(int ir=0; ir<12; ir++){
467 for(int ic=0; ic<5; ic++){
3d7e0d0c 468 read = fscanf(fdata,"%f ",&data[ic]);
469 if(read==0) AliDebug(3, " Error in reading PMT gains from external file ");
25175f38 470 }
144a0d1b 471 beam[ir] = (int) (data[0]);
472 det[ir] = (int) (data[1]);
25175f38 473 gain[ir] = data[2];
474 aEne[ir] = data[3];
475 bEne[ir] = data[4];
476 }
3d7ab74d 477 fclose(fdata);
25175f38 478
03f508de 479 if(((fBeamType.CompareTo("P-P")) == 0) || ((fBeamType.CompareTo("p-p")) == 0)){
25175f38 480 for(int i=0; i<12; i++){
481 if(beam[i]==0 && fBeamEnergy!=0.){
482 if(det[i]!=31 && det[i]!=32){
483 for(Int_t j=0; j<5; j++) fPMGain[det[i]-1][j] = gain[i]*(aEne[i]/fBeamEnergy+bEne[i]);
484 }
485 else if(det[i] == 31) fPMGain[2][1] = gain[i]*(aEne[i]-fBeamEnergy*bEne[i]);
486 else if(det[i] == 32) fPMGain[2][2] = gain[i]*(aEne[i]-fBeamEnergy*bEne[i]);
487 }
488 }
489 //
0e9277fd 490 printf("\n AliZDCDigitizer::ReadPMTGains -> ZDC PMT gains for p-p @ %1.0f+%1.0f GeV: ZNC(%1.0f), ZPC(%1.0f), ZEM(%1.0f), ZNA(%1.0f) ZPA(%1.0f)\n",
491 fBeamEnergy, fBeamEnergy, fPMGain[0][0], fPMGain[1][0], fPMGain[2][1], fPMGain[3][0], fPMGain[4][0]);
25175f38 492 }
493 else if(((fBeamType.CompareTo("A-A")) == 0)){
494 for(int i=0; i<12; i++){
495 if(beam[i]==1){
496 Float_t scalGainFactor = fBeamEnergy/2760.;
497 if(det[i]!=31 && det[i]!=32){
498 for(Int_t j=0; j<5; j++) fPMGain[det[i]-1][j] = gain[i]/(aEne[i]*scalGainFactor);
499 }
500 else{
501 for(int iq=1; iq<3; iq++) fPMGain[2][iq] = gain[i]/(aEne[i]*scalGainFactor);
502 }
503 }
504 }
505 //
0e9277fd 506 printf("\n AliZDCDigitizer::ReadPMTGains -> ZDC PMT gains for Pb-Pb @ %1.0f+%1.0f A GeV: ZN(%1.0f), ZP(%1.0f), ZEM(%1.0f)\n",
507 fBeamEnergy, fBeamEnergy, fPMGain[0][0], fPMGain[1][0], fPMGain[2][1]);
25175f38 508 }
0e9277fd 509 else if(((fBeamType.CompareTo("p-A")) == 0) || ((fBeamType.CompareTo("P-A")) == 0)
510 || ((fBeamType.CompareTo("A-p")) == 0) || ((fBeamType.CompareTo("A-P")) == 0)){
25175f38 511 for(int i=0; i<12; i++){
512 if(beam[i]==0 && fBeamEnergy!=0.){
513 if(det[i]==1 || det[i]==2){
514 for(Int_t j=0; j<5; j++) fPMGain[det[i]-1][j] = gain[i]*(aEne[i]/fBeamEnergy+bEne[i]);
515 }
516 }
517 if(beam[i]==1){
518 Float_t scalGainFactor = fBeamEnergy/2760.;
519 Float_t npartScalingFactor = 208./15.;
520 if(det[i]==4 || det[i]==5){
521 for(Int_t j=0; j<5; j++) fPMGain[det[i]-1][j] = npartScalingFactor*gain[i]/(aEne[i]*scalGainFactor);
522 }
523 else if(det[i]==31 || det[i]==32){
524 for(int iq=1; iq<3; iq++) fPMGain[2][iq] = npartScalingFactor*gain[i]/(aEne[i]*scalGainFactor);
525 }
526 }
527 }
0e9277fd 528 printf("\n AliZDCDigitizer::ReadPMTGains -> ZDC PMT gains for p-Pb: ZNC(%1.0f), ZPC(%1.0f), ZEM(%1.0f), ZNA(%1.0f) ZPA(%1.0f)\n",
529 fPMGain[0][0], fPMGain[1][0], fPMGain[2][1], fPMGain[3][0], fPMGain[4][0]);
25175f38 530 }
531}
532
5d24ec95 533//_____________________________________________________________________________
534void AliZDCDigitizer::CalculatePMTGains()
535{
536// Calculate PMT gain according to beam type and beam energy
58b2fa8a 537 if( ((fBeamType.CompareTo("P-P")) == 0) || ((fBeamType.CompareTo("p-p"))) ){
5d24ec95 538 // PTM gains rescaled to beam energy for p-p
539 // New correction coefficients for PMT gains needed
540 // to reproduce experimental spectra (from Grazia Jul 2010)
541 if(fBeamEnergy != 0){
542 for(Int_t j = 0; j < 5; j++){
543 fPMGain[0][j] = 1.515831*(661.444/fBeamEnergy+0.000740671)*10000000;
544 fPMGain[1][j] = 0.674234*(864.350/fBeamEnergy+0.00234375)*10000000;
545 fPMGain[3][j] = 1.350938*(661.444/fBeamEnergy+0.000740671)*10000000;
546 fPMGain[4][j] = 0.678597*(864.350/fBeamEnergy+0.00234375)*10000000;
547 }
548 fPMGain[2][1] = 0.869654*(1.32312-0.000101515*fBeamEnergy)*10000000;
549 fPMGain[2][2] = 1.030883*(1.32312-0.000101515*fBeamEnergy)*10000000;
550 //
0e9277fd 551 printf("\n AliZDCDigitizer::CalculatePMTGains -> ZDC PMT gains for p-p @ %1.0f+%1.0f GeV: ZNC(%1.0f), ZPC(%1.0f), ZEM(%1.0f), ZNA(%1.0f) ZPA(%1.0f)\n",
552 fBeamEnergy, fBeamEnergy, fPMGain[0][0], fPMGain[1][0], fPMGain[2][1], fPMGain[3][0], fPMGain[4][0]);
5d24ec95 553
554 }
555 }
556 else if(((fBeamType.CompareTo("A-A")) == 0)){
557 // PTM gains for Pb-Pb @ 2.7+2.7 A TeV ***************
558 // rescaled for Pb-Pb @ 1.38+1.38 A TeV ***************
559 // Values corrected after 2010 Pb-Pb data taking (7/2/2011 - Ch.)
560 // Experimental data compared to EMD simulation for single nucleon peaks:
561 // ZN gains must be divided by 4, ZP gains by 10!
562 Float_t scalGainFactor = fBeamEnergy/2760.;
563 for(Int_t j = 0; j < 5; j++){
564 fPMGain[0][j] = 50000./(4*scalGainFactor); // ZNC
565 fPMGain[1][j] = 100000./(5*scalGainFactor); // ZPC
566 fPMGain[2][j] = 100000./scalGainFactor; // ZEM
567 fPMGain[3][j] = 50000./(4*scalGainFactor); // ZNA
568 fPMGain[4][j] = 100000./(5*scalGainFactor); // ZPA
569 }
0e9277fd 570 printf("\n AliZDCDigitizer::CalculatePMTGains -> ZDC PMT gains for Pb-Pb @ %1.0f+%1.0f A GeV: ZN(%1.0f), ZP(%1.0f), ZEM(%1.0f)\n",
571 fBeamEnergy, fBeamEnergy, fPMGain[0][0], fPMGain[1][0], fPMGain[2][1]);
5d24ec95 572 }
f6c61b40 573 else if(((fBeamType.CompareTo("p-A")) == 0) || ((fBeamType.CompareTo("P-A"))) ){
5d24ec95 574 // PTM gains for Pb-Pb @ 1.38+1.38 A TeV on side A
575 // PTM gains rescaled to beam energy for p-p on side C
576 // WARNING! Energies are set by hand for 2011 pA RUN!!!
312b8ed7 577 Float_t scalGainFactor = fBeamEnergy/2760.;
f6c61b40 578 Float_t npartScalingFactor = 208./15.;
5d24ec95 579
580 for(Int_t j = 0; j < 5; j++){
581 fPMGain[0][j] = 1.515831*(661.444/fBeamEnergy+0.000740671)*10000000; //ZNC (p)
582 fPMGain[1][j] = 0.674234*(864.350/fBeamEnergy+0.00234375)*10000000; //ZPC (p)
8f190cf4 583 if(j<2) fPMGain[2][j] = npartScalingFactor*100000./scalGainFactor; // ZEM (Pb)
312b8ed7 584 // Npart max scales from 400 in Pb-Pb to ~8 in pPb -> *40.
f6c61b40 585 fPMGain[3][j] = npartScalingFactor*50000/(4*scalGainFactor); // ZNA (Pb)
586 fPMGain[4][j] = npartScalingFactor*100000/(5*scalGainFactor); // ZPA (Pb)
5d24ec95 587 }
0e9277fd 588 printf("\n AliZDCDigitizer::CalculatePMTGains -> ZDC PMT gains for p-Pb: ZNC(%1.0f), ZPC(%1.0f), ZEM(%1.0f), ZNA(%1.0f) ZPA(%1.0f)\n",
589 fPMGain[0][0], fPMGain[1][0], fPMGain[2][1], fPMGain[3][0], fPMGain[4][0]);
5d24ec95 590 }
8f190cf4 591 else if(((fBeamType.CompareTo("A-p")) == 0) || ((fBeamType.CompareTo("A-P"))) ){
592 // PTM gains for Pb-Pb @ 1.38+1.38 A TeV on side
593 // PTM gains rescaled to beam energy for p-p on side C
594 // WARNING! Energies are set by hand for 2011 pA RUN!!!
595 Float_t scalGainFactor = fBeamEnergy/2760.;
596 Float_t npartScalingFactor = 208./15.;
597
598 for(Int_t j = 0; j < 5; j++){
599 fPMGain[3][j] = 1.350938*(661.444/fBeamEnergy+0.000740671)*10000000; //ZNA (p)
600 fPMGain[4][j] = 0.678597*(864.350/fBeamEnergy+0.00234375)*10000000; //ZPA (p)
601 // Npart max scales from 400 in Pb-Pb to ~8 in pPb -> *40.
602 fPMGain[1][j] = npartScalingFactor*50000/(4*scalGainFactor); // ZNC (Pb)
603 fPMGain[2][j] = npartScalingFactor*100000/(5*scalGainFactor); // ZPC (Pb)
604 }
605 fPMGain[2][1] = 0.869654*(1.32312-0.000101515*fBeamEnergy)*10000000; // ZEM (pp)
606 fPMGain[2][2] = 1.030883*(1.32312-0.000101515*fBeamEnergy)*10000000; // ZEM (pp)
607 printf("\n AliZDCDigitizer::CalculatePMTGains -> ZDC PMT gains for p-Pb: ZNC(%1.0f), ZPC(%1.0f), ZEM(%1.0f), ZNA(%1.0f) ZPA(%1.0f)\n",
608 fPMGain[0][0], fPMGain[1][0], fPMGain[2][1], fPMGain[3][0], fPMGain[4][0]);
609 }
5d24ec95 610}
611
8309c1ab 612//_____________________________________________________________________________
613void AliZDCDigitizer::Fragmentation(Float_t impPar, Int_t specN, Int_t specP,
614 Int_t &freeSpecN, Int_t &freeSpecP) const
615{
616// simulate fragmentation of spectators
617
8309c1ab 618 AliZDCFragment frag(impPar);
8309c1ab 619
620 // Fragments generation
92339a90 621 frag.GenerateIMF();
622 Int_t nAlpha = frag.GetNalpha();
8309c1ab 623
624 // Attach neutrons
98711e01 625 frag.AttachNeutrons();
92339a90 626 Int_t ztot = frag.GetZtot();
627 Int_t ntot = frag.GetNtot();
98711e01 628
629 // Removing fragments and alpha pcs
8309c1ab 630 freeSpecN = specN-ntot-2*nAlpha;
631 freeSpecP = specP-ztot-2*nAlpha;
98711e01 632
3848c666 633 // Removing deuterons
a0607c1f 634 Int_t ndeu = (Int_t) (freeSpecN*frag.DeuteronNumber());
3848c666 635 freeSpecN -= ndeu;
98711e01 636 freeSpecP -= ndeu;
3848c666 637 //
8309c1ab 638 if(freeSpecN<0) freeSpecN=0;
639 if(freeSpecP<0) freeSpecP=0;
640 AliDebug(2, Form("FreeSpn = %d, FreeSpp = %d", freeSpecN, freeSpecP));
641}
642
643//_____________________________________________________________________________
94ff33bf 644void AliZDCDigitizer::SpectatorSignal(Int_t SpecType, Int_t numEvents, Float_t pm[5][5])
8309c1ab 645{
646// add signal of the spectators
fd9afd60 647
94ff33bf 648 if(!fSpectatorData) fSpectatorData = TFile::Open("$ALICE_ROOT/ZDC/SpectatorSignal.root");
649 if(!fSpectatorData || !fSpectatorData->IsOpen()) {
2d9c70ab 650 AliError((" Opening file $ALICE_ROOT/ZDC/SpectatorSignal.root failed\n"));
651 return;
8309c1ab 652 }
2d9c70ab 653
e7eabc21 654 TNtuple* zdcSignal=0x0;
2d9c70ab 655
656 Float_t sqrtS = 2*fBeamEnergy;
657 //
658 if(TMath::Abs(sqrtS-5500) < 100.){
5d24ec95 659 AliInfo(" Extracting signal from SpectatorSignal/energy5500 ");
94ff33bf 660 fSpectatorData->cd("energy5500");
2d9c70ab 661 //
662 if(SpecType == 1) { // --- Signal for projectile spectator neutrons
94ff33bf 663 fSpectatorData->GetObject("energy5500/ZNCSignal;1",zdcSignal);
2d9c70ab 664 if(!zdcSignal) AliError(" PROBLEM!!! Can't retrieve ZNCSignal from SpectatorSignal.root file");
665 }
666 else if(SpecType == 2) { // --- Signal for projectile spectator protons
94ff33bf 667 fSpectatorData->GetObject("energy5500/ZPCSignal;1",zdcSignal);
2d9c70ab 668 if(!zdcSignal) AliError(" PROBLEM!!! Can't retrieve ZPCSignal from SpectatorSignal.root file");
669 }
670 else if(SpecType == 3) { // --- Signal for target spectator neutrons
94ff33bf 671 fSpectatorData->GetObject("energy5500/ZNASignal;1",zdcSignal);
2d9c70ab 672 if(!zdcSignal) AliError(" PROBLEM!!! Can't retrieve ZNASignal from SpectatorSignal.root file");
673 }
674 else if(SpecType == 4) { // --- Signal for target spectator protons
94ff33bf 675 fSpectatorData->GetObject("energy5500/ZPASignal;1",zdcSignal);
2d9c70ab 676 if(!zdcSignal) AliError(" PROBLEM!!! Can't retrieve ZPASignal from SpectatorSignal.root file");
677 }
fd9afd60 678 }
2d9c70ab 679 else if(TMath::Abs(sqrtS-2760) < 100.){
5d24ec95 680 AliInfo(" Extracting signal from SpectatorSignal/energy2760 ");
94ff33bf 681 fSpectatorData->cd("energy2760");
2d9c70ab 682 //
683 if(SpecType == 1) { // --- Signal for projectile spectator neutrons
94ff33bf 684 fSpectatorData->GetObject("energy2760/ZNCSignal;1",zdcSignal);
2d9c70ab 685 if(!zdcSignal) AliError(" PROBLEM!!! Can't retrieve ZNCSignal from SpectatorSignal.root file");
686 }
687 else if(SpecType == 2) { // --- Signal for projectile spectator protons
94ff33bf 688 fSpectatorData->GetObject("energy2760/ZPCSignal;1",zdcSignal);
2d9c70ab 689 if(!zdcSignal) AliError(" PROBLEM!!! Can't retrieve ZPCSignal from SpectatorSignal.root file");
690 }
691 else if(SpecType == 3) { // --- Signal for target spectator neutrons
94ff33bf 692 fSpectatorData->GetObject("energy2760/ZNASignal;1",zdcSignal);
2d9c70ab 693 if(!zdcSignal) AliError(" PROBLEM!!! Can't retrieve ZNASignal from SpectatorSignal.root file");
694 }
695 else if(SpecType == 4) { // --- Signal for target spectator protons
94ff33bf 696 fSpectatorData->GetObject("energy2760/ZPASignal;1",zdcSignal);
2d9c70ab 697 if(!zdcSignal) AliError(" PROBLEM!!! Can't retrieve ZPASignal from SpectatorSignal.root file");
698 }
fd9afd60 699 }
700
83534754 701 if(!zdcSignal){
702 printf("\n No spectator signal available for ZDC digitization\n");
703 return;
704 }
705
8309c1ab 706 Int_t nentries = (Int_t) zdcSignal->GetEntries();
707
fd9afd60 708 Float_t *entry;
709 Int_t pl, i, k, iev=0, rnd[125], volume[2];
48642b09 710 for(pl=0;pl<125;pl++) rnd[pl] = 0;
8574b308 711 if(numEvents > 125) {
fea9b334 712 AliDebug(2,Form("numEvents (%d) is larger than 125", numEvents));
8309c1ab 713 numEvents = 125;
714 }
715 for(pl=0;pl<numEvents;pl++){
716 rnd[pl] = (Int_t) (9999*gRandom->Rndm());
85a96223 717 if(rnd[pl] >= 9999) rnd[pl] = 9998;
8309c1ab 718 //printf(" rnd[%d] = %d\n",pl,rnd[pl]);
719 }
720 // Sorting vector in ascending order with C function QSORT
721 qsort((void*)rnd,numEvents,sizeof(Int_t),comp);
722 do{
723 for(i=0; i<nentries; i++){
724 zdcSignal->GetEvent(i);
725 entry = zdcSignal->GetArgs();
726 if(entry[0] == rnd[iev]){
727 for(k=0; k<2; k++) volume[k] = (Int_t) entry[k+1];
48642b09 728 //
fd9afd60 729 Float_t lightQ = entry[7];
730 Float_t lightC = entry[8];
731 //
acf9550c 732 if(volume[0] != 3) { // ZN or ZP
8309c1ab 733 pm[volume[0]-1][0] += lightC;
734 pm[volume[0]-1][volume[1]] += lightQ;
48642b09 735 //printf("\n pm[%d][0] = %.0f, pm[%d][%d] = %.0f\n",(volume[0]-1),pm[volume[0]-1][0],
736 // (volume[0]-1),volume[1],pm[volume[0]-1][volume[1]]);
737 }
738 else {
8574b308 739 if(volume[1] == 1) pm[2][1] += lightC; // ZEM 1
acf9550c 740 else pm[2][2] += lightQ; // ZEM 2
48642b09 741 //printf("\n pm[2][1] = %.0f, pm[2][2] = %.0f\n",pm[2][1],pm[2][2]);
742 }
8309c1ab 743 }
744 else if(entry[0] > rnd[iev]){
745 iev++;
746 continue;
747 }
748 }
749 }while(iev<numEvents);
750
8309c1ab 751}
752
753
754//_____________________________________________________________________________
755Int_t AliZDCDigitizer::Phe2ADCch(Int_t Det, Int_t Quad, Float_t Light,
756 Int_t Res) const
757{
48642b09 758 // Evaluation of the ADC channel corresponding to the light yield Light
8a2624cc 759 Int_t vADCch = (Int_t) (Light * fPMGain[Det-1][Quad] * fADCRes[Res]);
698e394d 760 // Ch. debug
761 //printf("\t Phe2ADCch -> det %d quad %d - PMGain[%d][%d] %1.0f phe %1.2f ADC %d\n",
762 // Det,Quad,Det-1,Quad,fPMGain[Det-1][Quad],Light,vADCch);
7f73eb6b 763
8a2624cc 764 return vADCch;
48642b09 765}
8309c1ab 766
48642b09 767//_____________________________________________________________________________
768Int_t AliZDCDigitizer::Pedestal(Int_t Det, Int_t Quad, Int_t Res) const
769{
8a2624cc 770 // Returns a pedestal for detector det, PM quad, channel with res.
771 //
65b16e87 772 Float_t pedValue;
abf60186 773 // Normal run
774 if(fIsCalibration == 0){
c35ed519 775 Int_t index=0, kNch=24;
83347831 776 if(Quad!=5){
c35ed519 777 if(Det==1) index = Quad+kNch*Res; // ZNC
778 else if(Det==2) index = (Quad+5)+kNch*Res; // ZPC
779 else if(Det==3) index = (Quad+9)+kNch*Res; // ZEM
780 else if(Det==4) index = (Quad+12)+kNch*Res; // ZNA
781 else if(Det==5) index = (Quad+17)+kNch*Res; // ZPA
83347831 782 }
c35ed519 783 else index = (Det-1)/3+22+kNch*Res; // Reference PMs
83347831 784 //
c35ed519 785 Float_t meanPed = fPedData->GetMeanPed(index);
786 Float_t pedWidth = fPedData->GetMeanPedWidth(index);
65b16e87 787 pedValue = gRandom->Gaus(meanPed,pedWidth);
abf60186 788 //
58e12fee 789 /*printf("\t AliZDCDigitizer::Pedestal -> det %d quad %d res %d - Ped[%d] = %d\n",
65b16e87 790 Det, Quad, Res, index,(Int_t) pedValue); // Chiara debugging!
abf60186 791 */
792 }
7f73eb6b 793 // To create calibration object
83347831 794 else{
65b16e87 795 if(Res == 0) pedValue = gRandom->Gaus((35.+10.*gRandom->Rndm()),(0.5+0.2*gRandom->Rndm())); //High gain
796 else pedValue = gRandom->Gaus((250.+100.*gRandom->Rndm()),(3.5+2.*gRandom->Rndm())); //Low gain
83347831 797 }
dbc8d7fc 798
65b16e87 799 return (Int_t) pedValue;
8309c1ab 800}
801
802//_____________________________________________________________________________
78d18275 803AliCDBStorage* AliZDCDigitizer::SetStorage(const char *uri)
8309c1ab 804{
8309c1ab 805
78d18275 806 Bool_t deleteManager = kFALSE;
48642b09 807
78d18275 808 AliCDBManager *manager = AliCDBManager::Instance();
809 AliCDBStorage *defstorage = manager->GetDefaultStorage();
48642b09 810
78d18275 811 if(!defstorage || !(defstorage->Contains("ZDC"))){
812 AliWarning("No default storage set or default storage doesn't contain ZDC!");
813 manager->SetDefaultStorage(uri);
814 deleteManager = kTRUE;
815 }
816
817 AliCDBStorage *storage = manager->GetDefaultStorage();
818
819 if(deleteManager){
820 AliCDBManager::Instance()->UnsetDefaultStorage();
821 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
822 }
823
824 return storage;
825}
48642b09 826
78d18275 827//_____________________________________________________________________________
6024ec85 828AliZDCPedestals* AliZDCDigitizer::GetPedData() const
829{
830
831 // Getting pedestal calibration object for ZDC set
832
833 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
834 if(!entry) AliFatal("No calibration data loaded!");
be7fc6f8 835 else{
836 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
837 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
6024ec85 838
be7fc6f8 839 return calibdata;
840 }
6024ec85 841}
f91d1e35 842