e23e3600d0e473590667d1b8767f24bb240bf4cb
[u/mrichter/AliRoot.git] / ZDC / AliZDCMerger.cxx
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     fMHits       = 0;
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 }