]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ZDC/AliZDCMerger.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
File paths corrected
[u/mrichter/AliRoot.git] / ZDC / AliZDCMerger.cxx
... / ...
CommitLineData
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>
28#include <TH1.h>
29#include <TSystem.h>
30
31// --- AliRoot header files
32#include "AliZDCMerger.h"
33#include "AliZDC.h"
34#include "AliZDCHit.h"
35#include "AliZDCDigit.h"
36#include "AliZDCFragment.h"
37#include "AliRun.h"
38#include "AliDetector.h"
39#include "AliHeader.h"
40#include "AliGenEventHeader.h"
41#include "AliGenHijingEventHeader.h"
42
43ClassImp(AliZDCMerger)
44
45//____________________________________________________________________________
46AliZDCMerger::AliZDCMerger()
47{
48// Default constructor
49 fMerge = kDigitize;
50 fNEvBgr = 0;
51 fFnBgr = 0;
52 fBgrFile = 0;
53 fTrHBgr = 0;
54 fTrSDBgr = 0;
55 fImpPar = 0;
56 fSpecn = 0;
57 fSpecp = 0;
58 fFreeSpn = 0;
59 fFreeSpp = 0;
60 fFnSpecn = 0;
61 fSpecnFile = 0;
62 fFnSpecp = 0;
63 fSpecpFile = 0;
64
65}
66
67//____________________________________________________________________________
68AliZDCMerger::~AliZDCMerger()
69{
70// Destructor
71 if (fBgrFile) delete fBgrFile;
72 if (fTrHBgr) delete fTrHBgr;
73 if (fTrSDBgr) delete fTrSDBgr;
74 if (fHitsBgr) delete fHitsBgr;
75 if (fSpecnFile) delete fSpecnFile;
76 if (fSpecpFile) delete fSpecpFile;
77}
78
79//____________________________________________________________________________
80void AliZDCMerger::InitMerging()
81{
82 // Hits tree, impact parameter, num. of spectators n & p
83 // in background (full Hijing) event
84 Float_t b;
85 Int_t nspecn, nspecp;
86 Background(b, nspecn, nspecp);
87
88 // Production of nuclear fragments -> num. of FREE spectators n & p
89 Fragmentation(b, nspecn, nspecp, fFreeSpn, fFreeSpp);
90
91 // Extract from spectators distribution the signal events:
92 // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p
93 Mixing(fFreeSpn, fFreeSpp);
94}
95
96//____________________________________________________________________________
97void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
98{
99
100 // --- Open the background file
101 if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
102
103 // --- Read from the TreeE impact parameter (b),
104 // # of spectators n and p (fSpecn, fSpecp)
105 fBgrFile->cd();
106
107 // Get AliRun object from file or create it if not on file
108 gAlice = (AliRun*)fBgrFile->Get("gAlice");
109 if (!gAlice) {
110 gAlice = (AliRun*)fBgrFile->Get("gAlice");
111 if (gAlice) printf("AliRun object found on file\n");
112 if (!gAlice) {
113 printf("\n create new gAlice object");
114 gAlice = new AliRun("gAlice","Alice test program");
115 }
116 }
117
118// gAlice = (AliRun*)fBgrFile->Get("gAlice");
119 AliHeader *header = gAlice->GetHeader();
120 AliGenEventHeader* mcHeader = header->GenEventHeader();
121 fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
122 fSpecn = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
123 fSpecp = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
124// printf("\n HIJING simulation - b = %f fm, Nspecn = %d, Nspecp = %d\n",fImpPar,fSpecn,fSpecp);
125}
126
127//____________________________________________________________________________
128TFile* AliZDCMerger::OpenBgrFile()
129{
130 // Initialise background event
131 TFile *file = new TFile(fFnBgr,"UPDATE");
132// printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
133 fHitsBgr = new TClonesArray("AliZDCHit",1000);
134 return file;
135}
136
137//____________________________________________________________________________
138void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
139 Int_t &fFreeSpn, Int_t &fFreeSpp)
140{
141 Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot;
142 AliZDCFragment *frag = new AliZDCFragment(fImpPar);
143 for(j=0; j<=99; j++){
144 zz[j] =0;
145 nn[j] =0;
146 }
147
148 // Fragments generation
149 frag->GenerateIMF(zz, nAlpha);
150
151 // Attach neutrons
152 Ztot=0;
153 Ntot=0;
154 frag->AttachNeutrons(zz, nn, Ztot, Ntot);
155 fFreeSpn = fSpecn-Ztot-2*nAlpha;
156 fFreeSpp = fSpecp-Ntot-2*nAlpha;
157// printf("\n Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
158}
159
160//____________________________________________________________________________
161void AliZDCMerger::Mixing(Int_t fFreeSpn, Int_t fFreeSpp)
162{
163
164// printf("\n AliZDCMerger->Mixing\n");
165
166 // ### Background event Hits ###########################################
167 fBgrFile->cd();
168// fBgrFile->ls();
169
170 AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
171// if(ZDC) printf("\n Ho trovato lo ZDC!\n");
172
173 // Hits tree
174 if(fTrHBgr) delete fTrHBgr;
175 fTrHBgr = 0;
176 // SDigits tree
177 if(fTrSDBgr) delete fTrSDBgr;
178 fTrSDBgr = 0;
179
180// // Read from Event Tree the # of events of the Hijing file
181// TTree *thead = (TTree *) fBgrFile->Get("TE");
182// if(!te){
183// printf("\n ERROR! -> File %s does not contain TreeE\n");
184// }
185// fNEvBgr = (Int_t) te->GetEntries();
186 fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
187
188 // Hits Tree
189 char treeBgrName[20];
190 sprintf(treeBgrName,"TreeH%d",fNEvBgr);
191 fTrHBgr = (TTree*)gDirectory->Get(treeBgrName); // TreeH
192 if(!fTrHBgr){
193 printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
194 }
195 Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
196// printf("\n --- ntracks = %d\n",ntracks);
197
198 // SDigits Tree
199 char treeSBgrName[20];
200 sprintf(treeSBgrName,"TreeS%d",fNEvBgr);
201 fTrSDBgr = (TTree*)gDirectory->Get(treeSBgrName); // TreeH
202 if(!fTrSDBgr){
203 printf("\n ERROR -> Can't find TreeS%d in background file\n",fNEvBgr);
204 }
205
206 Int_t itrack, i, volume[2], detector, quadrant;
207 TClonesArray &sdigits = *fHitsBgr; // SDigits TCArray
208 Float_t hits[10];
209 // --- Tracks loop
210 for(itrack=0; itrack<ntracks; itrack++){
211// printf(" itrack = %d \n", itrack);
212 fTrHBgr->GetEvent(itrack);
213
214 fTrack = itrack;
215 Int_t NMhits = 0;
216 for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);
217 zdcHit;
218 zdcHit = (AliZDCHit*)ZDC->NextHit()){
219
220 for(i=0; i<2; i++) volume[i] = zdcHit->GetVolume(i);
221// printf("\n volume[0] = %d volume[1]= %d \n",volume[0], volume[1]);
222 if(volume[0] != 3){
223 // Background hits
224 hits[7] = zdcHit->GetLightPMQ();
225 hits[8] = zdcHit->GetLightPMC();
226// printf("\n Prima ### Background -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
227
228 //Signal hits
229 detector = volume[0];
230 quadrant = volume[1];
231 ExtractSignal(detector, quadrant, fQuadLight, fComLight);
232 hits[7] += fQuadLight;
233 hits[8] += fComLight;
234// printf("\n Bckg + Signal -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
235
236 zdcHit->SetLightPMQ(hits[7]);
237 zdcHit->SetLightPMC(hits[8]);
238 hits[7] = zdcHit->GetLightPMQ();
239 hits[8] = zdcHit->GetLightPMC();
240// printf("\n Dopo ### Background -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
241
242 new (sdigits[NMhits++]) AliZDCHit(zdcHit);
243// printf("\n NMhits = %d\n",NMhits);
244 }
245 fBgrFile->cd();
246 fTrSDBgr->Fill();
247 fTrSDBgr->Write(0,TObject::kOverwrite);
248 } //Hits loop
249 Digitize();
250 } // Tracks loop
251
252// fBgrFile->Close();
253
254}
255
256//____________________________________________________________________________
257void AliZDCMerger::ExtractSignal(Int_t detector, Int_t quadrant,
258 Float_t &fQuadLight, Float_t &fComLight)
259{
260
261// printf("\n Entering ExtractSignal method -> detector = %d quadrant = %d\n",
262// detector, quadrant);
263
264 // Connect spectator n histo's file
265 fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignal.root");
266 fSpecnFile = TFile::Open(fFnSpecn,"R");
267 fSpecnFile->cd();
268// fSpecnFile->ls();
269 TH1F *hPMQ1zn = (TH1F*) gDirectory->Get("hPMQ1zn;1");
270 TH1F *hPMQ2zn = (TH1F*) gDirectory->Get("hPMQ2zn;1");
271 TH1F *hPMQ3zn = (TH1F*) gDirectory->Get("hPMQ3zn;1");
272 TH1F *hPMQ4zn = (TH1F*) gDirectory->Get("hPMQ4zn;1");
273 TH1F *hPMC1zn = (TH1F*) gDirectory->Get("hPMC1zn;1");
274 TH1F *hPMC2zn = (TH1F*) gDirectory->Get("hPMC2zn;1");
275 TH1F *hPMC3zn = (TH1F*) gDirectory->Get("hPMC3zn;1");
276 TH1F *hPMC4zn = (TH1F*) gDirectory->Get("hPMC4zn;1");
277// Axis_t x = hPMQ1zn -> GetRandom();
278// printf(" hPMQ1zn -> GetRandom() = %f\n",x);
279
280 // Connect spectator p histo's file
281 fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignal.root");
282 fSpecpFile = TFile::Open(fFnSpecp,"R");
283 fSpecpFile->cd();
284// fSpecpFile->ls();
285 TH1F *hPMQ1zp = (TH1F*) gDirectory->Get("hPMQ1zp;1");
286 TH1F *hPMQ2zp = (TH1F*) gDirectory->Get("hPMQ2zp;1");
287 TH1F *hPMQ3zp = (TH1F*) gDirectory->Get("hPMQ3zp;1");
288 TH1F *hPMQ4zp = (TH1F*) gDirectory->Get("hPMQ4zp;1");
289 TH1F *hPMC1zp = (TH1F*) gDirectory->Get("hPMC1zp;1");
290 TH1F *hPMC2zp = (TH1F*) gDirectory->Get("hPMC2zp;1");
291 TH1F *hPMC3zp = (TH1F*) gDirectory->Get("hPMC3zp;1");
292 TH1F *hPMC4zp = (TH1F*) gDirectory->Get("hPMC4zp;1");
293
294 if(detector == 1){ // --- ZN
295 if(quadrant == 1){
296 fQuadLight = (Float_t) hPMQ1zn -> GetRandom();
297 fComLight = (Float_t) hPMC1zn -> GetRandom();
298 }
299 else if(quadrant == 2){
300 fQuadLight = (Float_t) hPMQ2zn -> GetRandom();
301 fComLight = (Float_t) hPMC2zn -> GetRandom();
302 }
303 else if(quadrant == 3){
304 fQuadLight = (Float_t) hPMQ3zn -> GetRandom();
305 fComLight = (Float_t) hPMC3zn -> GetRandom();
306 }
307 else if(quadrant == 4){
308 fQuadLight = (Float_t) hPMQ4zn -> GetRandom();
309 fComLight = (Float_t) hPMC4zn -> GetRandom();
310 }
311 }
312
313 else if(detector == 2){ // --- ZP
314 fSpecpFile->cd(); // Connect spectator p histo's file
315// fSpecpFile->ls();
316 if(quadrant == 1){
317 fQuadLight = (Float_t) hPMQ1zp -> GetRandom();
318 fComLight = (Float_t) hPMC1zp -> GetRandom();
319 }
320 else if(quadrant == 2){
321 fQuadLight = (Float_t) hPMQ2zp -> GetRandom();
322 fComLight = (Float_t) hPMC2zp -> GetRandom();
323 }
324 else if(quadrant == 3){
325 fQuadLight = (Float_t) hPMQ3zp -> GetRandom();
326 fComLight = (Float_t) hPMC3zp -> GetRandom();
327 }
328 else if(quadrant == 4){
329 fQuadLight = (Float_t) hPMQ4zp -> GetRandom();
330 fComLight = (Float_t) hPMC4zp -> GetRandom();
331 }
332 }
333// printf(" --- Exiting ExtractSignal -> fQuadLight = %f, fComLight = %f",
334// fQuadLight,fComLight);
335
336}
337
338//____________________________________________________________________________
339void AliZDCMerger::Digitize()
340{
341
342// printf("\n AliZDCMerger->Digitize()");
343
344 AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
345// if(ZDC) printf("\n Ho trovato lo ZDC!\n");
346
347 Int_t itrack, lightQ, lightC, sector[2], digit;
348 Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM = 0;
349 Int_t i;
350 for(i=0; i<4; i++){
351 PMQZN[i] = 0;
352 PMQZP[i] = 0;
353 }
354
355 // ### Digitization "on the flight" (no merging) #######################
356 if(fMerge == kDigitize){
357// printf("\n fMerge == kDigitize\n");
358 fTrHBgr = gAlice->TreeH(); // TTree
359 TTree *treeH = gAlice->TreeH();
360 if(!fTrHBgr){
361 printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
362 }
363 Int_t ntracks = (Int_t) treeH->GetEntries();
364// printf("\n --- ntracks = %d\n",ntracks);
365
366 // Loop over tracks
367 for(itrack=0; itrack<ntracks; itrack++){
368 gAlice->ResetHits();
369 gAlice->TreeH()->GetEvent(itrack);
370
371 // Loop over hits
372 for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);
373 zdcHit;
374 zdcHit = (AliZDCHit*)ZDC->NextHit()){
375 sector[0] = zdcHit->GetVolume(0);
376 sector[1] = zdcHit->GetVolume(1);
377 lightQ = Int_t(zdcHit->GetLightPMQ());
378 lightC = Int_t(zdcHit->GetLightPMC());
379// printf("\n Digitise -> DET. = %d, quad = %d", sector[0], sector[1]);
380// printf(" PMQ = %d, PMC = %d", lightQ, lightC);
381
382 if(sector[0] == 1){ //ZN
383 PMCZN = PMCZN + lightC;
384 PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
385 }
386 else if(sector[0] == 2){ //ZP
387 PMCZP = PMCZP + lightC;
388 PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
389 }
390 else if(sector[0] == 3){ //ZEM
391 PMZEM = PMZEM + lightC;
392 }
393 } // hits loop
394 } // tracks loop
395 } // if(fMerge)
396
397 // ### Merging and digitization #######################################
398 else if(fMerge == kMerge){
399// printf("\n fMerge == kMerge\n");
400
401// Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
402// printf("\n --- ntracks = %d\n",ntracks);
403//
404// // Loop over tracks
405// for(itrack=0; itrack<ntracks; itrack++){
406 fTrHBgr->GetEvent(fTrack);
407
408// printf("\n\n Track # %d --- Digitise -> ", fTrack);
409 // Loop over hits
410 for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);
411 zdcHit;
412 zdcHit = (AliZDCHit*)ZDC->NextHit()){
413 sector[0] = zdcHit->GetVolume(0);
414 sector[1] = zdcHit->GetVolume(1);
415 lightQ = Int_t(zdcHit->GetLightPMQ());
416 lightC = Int_t(zdcHit->GetLightPMC());
417// printf("\n DET. = %d, quad = %d,PMQ = %d, PMC = %d",
418// sector[0], sector[1],lightQ, lightC);
419
420 if(sector[0] == 1){ //ZN
421 PMCZN = PMCZN + lightC;
422 PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
423 }
424 else if(sector[0] == 2){ //ZP
425 PMCZP = PMCZP + lightC;
426 PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
427 }
428 else if(sector[0] == 3){ //ZEM
429 PMZEM = PMZEM + lightC;
430 }
431 } // hits loop
432// } // tracks loop
433 } // if(fMerge)
434
435
436 // Create digits for ZN
437 Int_t PedValue;
438 sector[0] = 1; // Detector = ZN
439 sector[1] = 0; // Common PM ADC
440 digit = Phe2ADCch(1, 0, PMCZN);
441// printf("\n\n ZN ### PMCZN = %d ADCZN = %d",PMCZN, digit);
442 ZDC->AddDigit(sector, digit);
443 Int_t j;
444 for(j=0; j<4; j++){
445 sector[1] = j+1; // Towers PM ADCs
446 digit = Phe2ADCch(1, j+1, PMQZN[j]);
447// printf("\n PMQZN[%d] = %d ADCZN[%d] = %d",j,PMQZN[j],j,digit);
448 PedValue = AddPedestal();
449 digit += PedValue;
450// printf(" PedValue = %d",PedValue);
451 ZDC->AddDigit(sector, digit);
452 }
453// printf("\n");
454
455 // Create digits for ZP
456 sector[0] = 2; // Detector = ZP
457 sector[1] = 0; // Common PM ADC
458 digit = Phe2ADCch(2, 0, PMCZP);
459// printf("\n ZP --- PMCZP = %d ADCZP = %d",PMCZP,digit);
460 ZDC->AddDigit(sector, digit);
461 for(j=0; j<4; j++){
462 sector[1] = j+1; // Towers PM ADCs
463 digit = Phe2ADCch(2, j+1, PMQZP[j]);
464// printf("\n PMQZP[%d] = %d ADCZP[%d] = %d",j,PMQZP[j],j,digit);
465 PedValue = AddPedestal();
466 digit += PedValue;
467// printf(" PedValue = %d",PedValue);
468 ZDC->AddDigit(sector, digit);
469 }
470// printf("\n");
471
472 // Create digits for ZEM
473 sector[0] = 3; // Detector = ZZEM
474 sector[1] = 0; // Single PM ADC
475 digit = Phe2ADCch(3, 0, PMZEM);
476// printf("\n ZEM *** PMZEM = %d ADCZEM = %d",PMZEM,digit);
477 PedValue = AddPedestal();
478 digit += PedValue;
479// printf(" PedValue = %d\n",PedValue);
480 ZDC->AddDigit(sector, digit);
481
482}
483
484//_____________________________________________________________________________
485Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light)
486{
487 // Evaluation of the ADC channel corresponding to the light yield Light
488
489 if(gAlice->GetDebug() > 0){
490// printf("\n Phe2ADCch -> Detector = %d, Quadrant = %d, Light = %d\n", Det, Quad, Light);
491 }
492
493 Int_t ADCch = 0;
494 // Parameters for conversion of light yield in ADC channels
495 Float_t fPMGain[3][5]; // PM gain
496 Float_t fADCRes; // ADC conversion factor
497
498 Int_t j,i;
499 for(i=0; i<3; i++){
500 for(j=0; j<5; j++){
501 fPMGain[i][j] = 100000.;
502 }
503 }
504 fADCRes = 0.00000064; // ADC Resolution: 250 fC/ADCch
505
506 ADCch = (Int_t) (Light*fPMGain[Det-1][Quad]*fADCRes);
507
508 return ADCch;
509}
510
511//_____________________________________________________________________________
512Int_t AliZDCMerger::AddPedestal()
513{
514 // --- Pedestal value -> extracted from a gaussian distribution
515 // obtained from the beam test on the ZEM prototype (Aug. 2000)
516
517 Int_t PedValue;
518 Float_t PedMean = 50.;
519 Float_t PedWidth = 10.;
520
521 PedValue = (Int_t) gRandom->Gaus(PedMean,PedWidth);
522
523 return PedValue;
524}