Pointer fMHits set to 0 in the default constructor
[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
46//____________________________________________________________________________
47AliZDCMerger::AliZDCMerger()
48{
49// Default constructor
50 fMerge = kDigitize;
1acbaaa8 51// fMerge = kMerge;
359cdddc 52 fFnBgr = 0;
359cdddc 53 fBgrFile = 0;
1acbaaa8 54 fNEvBgr = 0;
359cdddc 55 fTrHBgr = 0;
2685bf00 56 fHitsBgr = 0;
359cdddc 57 fTrSDBgr = 0;
2685bf00 58 fTrDBgr = 0;
359cdddc 59 fImpPar = 0;
60 fSpecn = 0;
61 fSpecp = 0;
62 fFreeSpn = 0;
63 fFreeSpp = 0;
64 fFnSpecn = 0;
359cdddc 65 fSpecnFile = 0;
66 fFnSpecp = 0;
359cdddc 67 fSpecpFile = 0;
1acbaaa8 68 fNMhits = 0;
bfeb9a74 69 fMHits = 0;
359cdddc 70}
71
72//____________________________________________________________________________
73AliZDCMerger::~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//____________________________________________________________________________
85void 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
1acbaaa8 98 Mixing();
359cdddc 99}
100
101//____________________________________________________________________________
102void 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
1acbaaa8 112// // Get AliRun object from file or create it if not on file
359cdddc 113// gAlice = (AliRun*)fBgrFile->Get("gAlice");
1acbaaa8 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
359cdddc 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();
1acbaaa8 128 printf("\n HIJING simulation - b = %f fm, Nspecn = %d, Nspecp = %d\n",fImpPar,fSpecn,fSpecp);
359cdddc 129}
130
131//____________________________________________________________________________
132TFile* AliZDCMerger::OpenBgrFile()
133{
134 // Initialise background event
135 TFile *file = new TFile(fFnBgr,"UPDATE");
1acbaaa8 136 printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
359cdddc 137 fHitsBgr = new TClonesArray("AliZDCHit",1000);
1acbaaa8 138 fMHits = new TClonesArray("AliZDCMergedHit",1000);
359cdddc 139 return file;
140}
141
142//____________________________________________________________________________
143void 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;
1acbaaa8 162 printf("\n Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
359cdddc 163}
164
165//____________________________________________________________________________
1acbaaa8 166void AliZDCMerger::Mixing()
359cdddc 167{
1acbaaa8 168
169// printf("\n AliZDCMerger->Mixing\n");
170
359cdddc 171 // ### Background event Hits ###########################################
172 fBgrFile->cd();
173// fBgrFile->ls();
174
1acbaaa8 175 AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
359cdddc 176// if(ZDC) printf("\n Ho trovato lo ZDC!\n");
177
178 // Hits tree
179 if(fTrHBgr) delete fTrHBgr;
1acbaaa8 180 fTrHBgr = 0;
359cdddc 181 // SDigits tree
182 if(fTrSDBgr) delete fTrSDBgr;
183 fTrSDBgr = 0;
184
1acbaaa8 185// fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
359cdddc 186
1acbaaa8 187 // Hits tree
359cdddc 188 char treeBgrName[20];
189 sprintf(treeBgrName,"TreeH%d",fNEvBgr);
1acbaaa8 190 fTrHBgr = (TTree*)gDirectory->Get(treeBgrName);
359cdddc 191 if(!fTrHBgr){
192 printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
193 }
1acbaaa8 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
359cdddc 206 Int_t ntracks = (Int_t) fTrHBgr->GetEntries();
1acbaaa8 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;
359cdddc 215
359cdddc 216 // --- Tracks loop
217 for(itrack=0; itrack<ntracks; itrack++){
1acbaaa8 218// printf(" itrack = %d", itrack);
219// gAlice->ResetHits();
359cdddc 220 fTrHBgr->GetEvent(itrack);
221
1acbaaa8 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
359cdddc 248 } // Tracks loop
1acbaaa8 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);
359cdddc 258
1acbaaa8 259 // --- Protons
260 ExtractSignal(2);
261// printf("\n fNMhits (after signal) = %d \n",fNMhits);
359cdddc 262
263}
264
265//____________________________________________________________________________
1acbaaa8 266void AliZDCMerger::ExtractSignal(Int_t SpecType)
359cdddc 267{
268
1acbaaa8 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
359cdddc 326}
327
328//____________________________________________________________________________
1acbaaa8 329void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
359cdddc 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
1acbaaa8 337 Int_t lightQ, lightC, sector[2], digit;
338 Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM1 = 0, PMZEM2 = 0;
359cdddc 339 Int_t i;
340 for(i=0; i<4; i++){
341 PMQZN[i] = 0;
342 PMQZP[i] = 0;
343 }
1acbaaa8 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
359cdddc 377
1acbaaa8 378 // ### Digits creation ###############################################
359cdddc 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);
1acbaaa8 385 PedValue = AddPedestal();
386 digit += PedValue;
387// printf(" PedValue = %d",PedValue);
359cdddc 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]);
1acbaaa8 393// printf("\n PMQZN[%d] = %d phe ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit);
359cdddc 394 PedValue = AddPedestal();
395 digit += PedValue;
396// printf(" PedValue = %d",PedValue);
397 ZDC->AddDigit(sector, digit);
398 }
1acbaaa8 399// printf("\n");
359cdddc 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);
1acbaaa8 405// printf("\n ZP --- PMCZP = %d phe ADCZP = %d ADCch",PMCZP,digit);
406 PedValue = AddPedestal();
407 digit += PedValue;
408 printf(" PedValue = %d",PedValue);
359cdddc 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]);
1acbaaa8 413// printf("\n PMQZP[%d] = %d phe ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit);
359cdddc 414 PedValue = AddPedestal();
415 digit += PedValue;
416// printf(" PedValue = %d",PedValue);
417 ZDC->AddDigit(sector, digit);
418 }
1acbaaa8 419// printf("\n");
359cdddc 420
421 // Create digits for ZEM
1acbaaa8 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);
359cdddc 436 ZDC->AddDigit(sector, digit);
359cdddc 437}
438
439//_____________________________________________________________________________
440Int_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//_____________________________________________________________________________
467Int_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}