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