AliZDCDigits inherits from TObject
[u/mrichter/AliRoot.git] / ZDC / AliZDCMerger.cxx
CommitLineData
359cdddc 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
17///////////////////////////////////////////////////////////////////////////////
18// //
19// ZDC event merging class //
20// //
21///////////////////////////////////////////////////////////////////////////////
22
23// --- ROOT system
24#include <iostream.h>
25#include <TTree.h>
26#include <TFile.h>
27#include <TDirectory.h>
1acbaaa8 28#include <TNtuple.h>
53c30603 29#include <TSystem.h>
359cdddc 30
31// --- AliRoot header files
32#include "AliZDCMerger.h"
33#include "AliZDC.h"
34#include "AliZDCHit.h"
1acbaaa8 35#include "AliZDCMergedHit.h"
359cdddc 36#include "AliZDCDigit.h"
37#include "AliZDCFragment.h"
38#include "AliRun.h"
39#include "AliDetector.h"
40#include "AliHeader.h"
41#include "AliGenEventHeader.h"
53c30603 42#include "AliGenHijingEventHeader.h"
359cdddc 43
44ClassImp(AliZDCMerger)
45
1450a7cd 46//int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
47
359cdddc 48//____________________________________________________________________________
49AliZDCMerger::AliZDCMerger()
50{
51// Default constructor
1450a7cd 52 //fMerge = kDigitize -> Only digitization
53 //fMerge = kMerge -> Digitization + Merging
54 fMerge = kMerge;
359cdddc 55 fFnBgr = 0;
359cdddc 56 fBgrFile = 0;
1acbaaa8 57 fNEvBgr = 0;
359cdddc 58 fTrHBgr = 0;
2685bf00 59 fHitsBgr = 0;
359cdddc 60 fTrSDBgr = 0;
2685bf00 61 fTrDBgr = 0;
359cdddc 62 fImpPar = 0;
63 fSpecn = 0;
64 fSpecp = 0;
65 fFreeSpn = 0;
66 fFreeSpp = 0;
67 fFnSpecn = 0;
359cdddc 68 fSpecnFile = 0;
69 fFnSpecp = 0;
359cdddc 70 fSpecpFile = 0;
1acbaaa8 71 fNMhits = 0;
1450a7cd 72
359cdddc 73}
74
75//____________________________________________________________________________
76AliZDCMerger::~AliZDCMerger()
77{
78// Destructor
79 if (fBgrFile) delete fBgrFile;
80 if (fTrHBgr) delete fTrHBgr;
81 if (fTrSDBgr) delete fTrSDBgr;
82 if (fHitsBgr) delete fHitsBgr;
83 if (fSpecnFile) delete fSpecnFile;
84 if (fSpecpFile) delete fSpecpFile;
85}
86
87//____________________________________________________________________________
88void AliZDCMerger::InitMerging()
89{
90 // Hits tree, impact parameter, num. of spectators n & p
91 // in background (full Hijing) event
92 Float_t b;
1450a7cd 93 Background(b, fSpecn, fSpecp);
359cdddc 94
95 // Production of nuclear fragments -> num. of FREE spectators n & p
1450a7cd 96 Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp);
359cdddc 97
98 // Extract from spectators distribution the signal events:
99 // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p
1acbaaa8 100 Mixing();
359cdddc 101}
102
103//____________________________________________________________________________
104void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
105{
106
107 // --- Open the background file
108 if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
109
110 // --- Read from the TreeE impact parameter (b),
111 // # of spectators n and p (fSpecn, fSpecp)
112 fBgrFile->cd();
113
1acbaaa8 114// // Get AliRun object from file or create it if not on file
359cdddc 115// gAlice = (AliRun*)fBgrFile->Get("gAlice");
1acbaaa8 116// if (!gAlice) {
117// gAlice = (AliRun*)fBgrFile->Get("gAlice");
118// if (gAlice) printf("AliRun object found on file\n");
119// if (!gAlice) {
120// printf("\n create new gAlice object");
121// gAlice = new AliRun("gAlice","Alice test program");
122// }
123// }
1450a7cd 124
125 gAlice->GetEvent(fNEvBgr);
359cdddc 126 AliHeader *header = gAlice->GetHeader();
127 AliGenEventHeader* mcHeader = header->GenEventHeader();
128 fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
1450a7cd 129 Int_t dSpecn = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
130 Int_t dSpecp = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
131 // Until there is only 1 ZDC set the # of spectators must be divided by 2!!!
132 fSpecn = dSpecn/2;
133 fSpecp = dSpecp/2;
134 printf("\n HIJING ev. #%d - b = %f fm, Nspecn = %d, Nspecp = %d\n",
135 fNEvBgr,fImpPar,fSpecn,fSpecp);
359cdddc 136}
137
138//____________________________________________________________________________
139TFile* AliZDCMerger::OpenBgrFile()
140{
141 // Initialise background event
142 TFile *file = new TFile(fFnBgr,"UPDATE");
1acbaaa8 143 printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
359cdddc 144 fHitsBgr = new TClonesArray("AliZDCHit",1000);
1acbaaa8 145 fMHits = new TClonesArray("AliZDCMergedHit",1000);
359cdddc 146 return file;
147}
148
149//____________________________________________________________________________
150void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
151 Int_t &fFreeSpn, Int_t &fFreeSpp)
152{
1450a7cd 153 //printf("\n Fragmentation -> fSpecn = %d, fSpecp = %d\n",fSpecn,fSpecp);
359cdddc 154 Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot;
155 AliZDCFragment *frag = new AliZDCFragment(fImpPar);
156 for(j=0; j<=99; j++){
157 zz[j] =0;
158 nn[j] =0;
159 }
160
161 // Fragments generation
162 frag->GenerateIMF(zz, nAlpha);
163
164 // Attach neutrons
165 Ztot=0;
166 Ntot=0;
167 frag->AttachNeutrons(zz, nn, Ztot, Ntot);
1450a7cd 168 fFreeSpn = fSpecn-Ntot-2*nAlpha;
169 fFreeSpp = fSpecp-Ztot-2*nAlpha;
170 if(fFreeSpn<0) fFreeSpn=0;
171 if(fFreeSpp<0) fFreeSpp=0;
172 //printf("\n 2*nAlpha = %d, Ztot = %d, Ntot = %d\n",2*nAlpha, Ztot, Ntot);
1acbaaa8 173 printf("\n Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
359cdddc 174}
175
176//____________________________________________________________________________
1acbaaa8 177void AliZDCMerger::Mixing()
359cdddc 178{
1acbaaa8 179
1450a7cd 180 //printf("\n AliZDCMerger->Mixing\n");
1acbaaa8 181
359cdddc 182 // ### Background event Hits ###########################################
183 fBgrFile->cd();
184// fBgrFile->ls();
185
1acbaaa8 186 AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
359cdddc 187// if(ZDC) printf("\n Ho trovato lo ZDC!\n");
188
189 // Hits tree
190 if(fTrHBgr) delete fTrHBgr;
1acbaaa8 191 fTrHBgr = 0;
359cdddc 192 // SDigits tree
193 if(fTrSDBgr) delete fTrSDBgr;
194 fTrSDBgr = 0;
195
1acbaaa8 196// fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
1acbaaa8 197 // Hits tree
359cdddc 198 char treeBgrName[20];
199 sprintf(treeBgrName,"TreeH%d",fNEvBgr);
1acbaaa8 200 fTrHBgr = (TTree*)gDirectory->Get(treeBgrName);
359cdddc 201 if(!fTrHBgr){
202 printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
203 }
1acbaaa8 204// fTrHBgr->Print();
205
206 // Branch address
207 TBranch *branch;
208 char branchname[20];
209 sprintf(branchname,"%s",ZDC->GetName());
210 if(fTrHBgr && fHitsBgr){
211// printf("\n fTrHBgr!=0 && fHitsBgr!=0\n");
212 branch = fTrHBgr->GetBranch(branchname);
213 if(branch) branch->SetAddress(&fHitsBgr);
214 }
215
359cdddc 216 Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
1450a7cd 217 //printf("\n --- ntracks = %d\n\n", ntracks);
1acbaaa8 218
219 Int_t itrack, nhits, ihit, j, sector[2];
220 AliZDCHit* zdcHit;
221 AliZDCMergedHit *MHit;
222 Float_t MHits[7];
223 TClonesArray &sdigits = *fMHits; // SDigits TCArray
224 fNMhits = 0;
359cdddc 225
359cdddc 226 // --- Tracks loop
227 for(itrack=0; itrack<ntracks; itrack++){
359cdddc 228 fTrHBgr->GetEvent(itrack);
229
1acbaaa8 230 nhits = fHitsBgr->GetEntries();
1acbaaa8 231// printf(" nhits = %d \n", nhits);
232 for(ihit=0; ihit<nhits; ihit++){
233 zdcHit = (AliZDCHit*) fHitsBgr->UncheckedAt(ihit);
234
235 for(j=0; j<2; j++) sector[j] = zdcHit->GetVolume(j);
236 MHits[0] = zdcHit->GetPrimKinEn();
237 MHits[1] = zdcHit->GetXImpact();
238 MHits[2] = zdcHit->GetYImpact();
239 MHits[3] = zdcHit->GetSFlag();
1450a7cd 240 MHits[4] = zdcHit->GetLightPMQ();
1acbaaa8 241 MHits[5] = zdcHit->GetLightPMC();
242 MHits[6] = zdcHit->GetEnergy();
243 MHit = new AliZDCMergedHit(sector, MHits);
244// MHit->Print("");
245 new((*fMHits)[fNMhits]) AliZDCMergedHit(*MHit);
246 new (sdigits[fNMhits]) AliZDCMergedHit(*MHit);
247 delete MHit;
248 fNMhits++;
249 }//Hits loop
250
359cdddc 251 } // Tracks loop
1450a7cd 252 //printf(" fNMhits (after bckg) = %d, \n",fNMhits);
1acbaaa8 253
254 // ### Signal event Hits ###########################################
255 // --- Neutrons
256 ExtractSignal(1);
359cdddc 257
1acbaaa8 258 // --- Protons
259 ExtractSignal(2);
1450a7cd 260 //printf(" fNMhits (after signal) = %d \n",fNMhits);
359cdddc 261
262}
263
264//____________________________________________________________________________
1acbaaa8 265void AliZDCMerger::ExtractSignal(Int_t SpecType)
359cdddc 266{
267
1acbaaa8 268// printf("\n Entering in Extract Signal\n");
269
270 Int_t NumEvents = 0;
271 if(SpecType == 1){ // --- Signal for spectator neutrons
272 fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
273 fSpecnFile = TFile::Open(fFnSpecn,"R");
274 fSpecnFile->cd();
1450a7cd 275 printf("\n --- ExtractSignal x n: file %s opened\n", fFnSpecn);
1acbaaa8 276 NumEvents = fFreeSpn;
277 }
278 else if(SpecType == 2){ // --- Signal for spectator protons
279 fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
280 fSpecpFile = TFile::Open(fFnSpecp,"R");
281 fSpecpFile->cd();
1450a7cd 282 printf("\n --- ExtractSignal x p: file %s opened\n", fFnSpecp);
1acbaaa8 283 NumEvents = fFreeSpp;
284 }
1450a7cd 285 //printf("\n # of free spectators = %d\n", NumEvents);
286 //printf("\n fNMhits (before adding signal) = %d\n",fNMhits);
1acbaaa8 287
288 TNtuple *ZDCSignal = (TNtuple*) gDirectory->Get("ZDCSignal");
289 Int_t nentries = (Int_t) ZDCSignal->GetEntries();
1450a7cd 290 //printf("\n # entries = %d\n", nentries);
1acbaaa8 291
292 AliZDCMergedHit *MHit;
293 Float_t *entry, HitsSpec[7];
1450a7cd 294 Int_t pl, i, j, k, iev=0, rnd[125], Volume[2];
295 for(pl=0;pl<125;pl++){
296 rnd[pl] = 0;
297 }
298 for(pl=0;pl<NumEvents;pl++){
299 rnd[pl] = (Int_t) (9999*gRandom->Rndm());
300 if(rnd[pl] >= 9998) rnd[pl] = 9997;
301 //printf(" rnd[%d] = %d\n",pl,rnd[pl]);
302 }
303 // Sorting vector in ascending order with C function QSORT
304 qsort((void*)rnd,NumEvents,sizeof(Int_t),comp);
305 //for(pl=0;pl<NumEvents;pl++){
306 ////printf(" rnd[%d] = %d\n",pl,rnd[pl]);
307 //}
308 do{
1acbaaa8 309 for(i=0; i<nentries; i++){
310 ZDCSignal->GetEvent(i);
311 entry = ZDCSignal->GetArgs();
1450a7cd 312 if(entry[0] == rnd[iev]){
1acbaaa8 313 for(k=0; k<2; k++) Volume[k] = (Int_t) entry[k+1];
314 for(j=0; j<7; j++){
315 HitsSpec[j] = entry[j+3];
316 }
1450a7cd 317 //printf("\n i = %d, iev = %d, entry[0] = %f, rnd[%d] = %d ",i,iev,entry[0],iev,rnd[iev]);
1acbaaa8 318 MHit = new AliZDCMergedHit(Volume, HitsSpec);
319 new((*fMHits)[fNMhits++]) AliZDCMergedHit(*MHit);
320 delete MHit;
321 }
1450a7cd 322 else if(entry[0] > rnd[iev]){
323 iev++;
324 continue;
325 }
1acbaaa8 326 }
1450a7cd 327 }while(iev<NumEvents);
1acbaaa8 328
329 if(SpecType ==1){
1450a7cd 330 //printf("\n fNMhits (after n signal) = %d\n",fNMhits);
1acbaaa8 331 fSpecnFile->Close();
332 }
333 else if(SpecType == 2){
1450a7cd 334 //printf("\n fNMhits (after p signal) = %d\n",fNMhits);
1acbaaa8 335 fSpecpFile->Close();
336 }
337
359cdddc 338}
339
340//____________________________________________________________________________
1acbaaa8 341void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
359cdddc 342{
343
1450a7cd 344 printf("\n AliZDCMerger->Digitize()");
359cdddc 345
346 AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
347// if(ZDC) printf("\n Ho trovato lo ZDC!\n");
348
1acbaaa8 349 Int_t lightQ, lightC, sector[2], digit;
350 Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM1 = 0, PMZEM2 = 0;
359cdddc 351 Int_t i;
352 for(i=0; i<4; i++){
353 PMQZN[i] = 0;
354 PMQZP[i] = 0;
355 }
1acbaaa8 356
357 AliZDCMergedHit *MHit;
358 Int_t imhit;
1450a7cd 359 printf(" fNMHits = %d\n", fNMhits);
1acbaaa8 360 // Loop over SDigits
361 for(imhit=0; imhit<fNMhits; imhit++){
362
363 MHit = (AliZDCMergedHit*) fMHits->UncheckedAt(imhit);
364 sector[0] = MHit->GetSector(0);
365 sector[1] = MHit->GetSector(1);
1450a7cd 366 // tmp -> c'erano i quadranti cannati!
1acbaaa8 367 if((sector[1]!=1) && (sector[1]!=2) && (sector[1]!=3) && (sector[1]!=4)){
1450a7cd 368 printf("\n *** ERROR!!! sector[0] = %d, sector[1] = %d\n",
369 sector[0], sector[1]);
1acbaaa8 370 sector[1] = 0;
371 }
372 lightQ = Int_t(MHit->GetLightPMQ());
373 lightC = Int_t(MHit->GetLightPMC());
1450a7cd 374// printf(" imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d\n",
1acbaaa8 375// imhit,sector[0], sector[1],lightQ, lightC);
376
377 if(sector[0] == 1){ //ZN
378 PMCZN = PMCZN + lightC;
379 PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
380 }
381 else if(sector[0] == 2){ //ZP
382 PMCZP = PMCZP + lightC;
383 PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
384 }
385 else if(sector[0] == 3){ //ZEM
386 if(sector[1] ==1) PMZEM1 = PMZEM1 + lightC;
387 else PMZEM2 = PMZEM2 + lightQ;
388 }
389 } // SDigits loop
359cdddc 390
1acbaaa8 391 // ### Digits creation ###############################################
359cdddc 392 // Create digits for ZN
393 Int_t PedValue;
394 sector[0] = 1; // Detector = ZN
395 sector[1] = 0; // Common PM ADC
396 digit = Phe2ADCch(1, 0, PMCZN);
1450a7cd 397 printf("\n\n ZN ### PMCZN = %d ADCZN = %d",PMCZN, digit);
1acbaaa8 398 PedValue = AddPedestal();
399 digit += PedValue;
400// printf(" PedValue = %d",PedValue);
359cdddc 401 ZDC->AddDigit(sector, digit);
402 Int_t j;
403 for(j=0; j<4; j++){
404 sector[1] = j+1; // Towers PM ADCs
405 digit = Phe2ADCch(1, j+1, PMQZN[j]);
1450a7cd 406 printf("\n PMQZN[%d] = %d phe ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit);
359cdddc 407 PedValue = AddPedestal();
408 digit += PedValue;
409// printf(" PedValue = %d",PedValue);
410 ZDC->AddDigit(sector, digit);
411 }
1450a7cd 412 printf("\n");
359cdddc 413
414 // Create digits for ZP
415 sector[0] = 2; // Detector = ZP
416 sector[1] = 0; // Common PM ADC
417 digit = Phe2ADCch(2, 0, PMCZP);
1450a7cd 418 printf("\n ZP --- PMCZP = %d phe ADCZP = %d ADCch",PMCZP,digit);
1acbaaa8 419 PedValue = AddPedestal();
420 digit += PedValue;
1450a7cd 421// printf(" PedValue = %d",PedValue);
359cdddc 422 ZDC->AddDigit(sector, digit);
423 for(j=0; j<4; j++){
424 sector[1] = j+1; // Towers PM ADCs
425 digit = Phe2ADCch(2, j+1, PMQZP[j]);
1450a7cd 426 printf("\n PMQZP[%d] = %d phe ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit);
359cdddc 427 PedValue = AddPedestal();
428 digit += PedValue;
429// printf(" PedValue = %d",PedValue);
430 ZDC->AddDigit(sector, digit);
431 }
1450a7cd 432 printf("\n");
359cdddc 433
434 // Create digits for ZEM
1acbaaa8 435 sector[0] = 3;
436 sector[1] = 1; // Detector = ZEM1
437 digit = Phe2ADCch(3, 1, PMZEM1);
1450a7cd 438 printf("\n ZEM *** PMZEM1 = %d phe ADCZEM1 = %d ADCch",PMZEM1,digit);
1acbaaa8 439 PedValue = AddPedestal();
440 digit += PedValue;
441// printf(" PedValue = %d\n",PedValue);
442 ZDC->AddDigit(sector, digit);
443 sector[1] = 2; // Detector = ZEM2
444 digit = Phe2ADCch(3, 2, PMZEM2);
1450a7cd 445 printf("\n ZEM *** PMZEM2 = %d phe ADCZEM2 = %d ADCch\n",PMZEM2,digit);
1acbaaa8 446 PedValue = AddPedestal();
447 digit += PedValue;
448// printf(" PedValue = %d\n",PedValue);
359cdddc 449 ZDC->AddDigit(sector, digit);
359cdddc 450}
451
452//_____________________________________________________________________________
453Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light)
454{
455 // Evaluation of the ADC channel corresponding to the light yield Light
456
457 if(gAlice->GetDebug() > 0){
458// printf("\n Phe2ADCch -> Detector = %d, Quadrant = %d, Light = %d\n", Det, Quad, Light);
459 }
460
461 Int_t ADCch = 0;
359cdddc 462
463 Int_t j,i;
464 for(i=0; i<3; i++){
465 for(j=0; j<5; j++){
466 fPMGain[i][j] = 100000.;
467 }
468 }
469 fADCRes = 0.00000064; // ADC Resolution: 250 fC/ADCch
470
471 ADCch = (Int_t) (Light*fPMGain[Det-1][Quad]*fADCRes);
472
473 return ADCch;
474}
475
476//_____________________________________________________________________________
477Int_t AliZDCMerger::AddPedestal()
478{
479 // --- Pedestal value -> extracted from a gaussian distribution
480 // obtained from the beam test on the ZEM prototype (Aug. 2000)
481
482 Int_t PedValue;
483 Float_t PedMean = 50.;
1450a7cd 484 Float_t PedWidth = 5.;
359cdddc 485
486 PedValue = (Int_t) gRandom->Gaus(PedMean,PedWidth);
487
488 return PedValue;
489}