Merging and reconstruction code review
[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 //int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
47
48 //____________________________________________________________________________
49 AliZDCMerger::AliZDCMerger()
50 {
51 // Default constructor    
52     //fMerge       = kDigitize  -> Only digitization
53     //fMerge       = kMerge     -> Digitization + Merging
54     fMerge       = kMerge;      
55     fFnBgr       = 0;
56     fBgrFile     = 0;
57     fNEvBgr      = 0;
58     fTrHBgr      = 0;
59     fHitsBgr     = 0;
60     fTrSDBgr     = 0;
61     fTrDBgr      = 0;
62     fImpPar      = 0;
63     fSpecn       = 0;
64     fSpecp       = 0;
65     fFreeSpn     = 0;
66     fFreeSpp     = 0;
67     fFnSpecn     = 0;
68     fSpecnFile   = 0;
69     fFnSpecp     = 0;
70     fSpecpFile   = 0;
71     fNMhits      = 0;
72      
73 }
74
75 //____________________________________________________________________________
76 AliZDCMerger::~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 //____________________________________________________________________________
88 void AliZDCMerger::InitMerging()
89 {
90     // Hits tree, impact parameter, num. of spectators n & p 
91     //          in background (full Hijing) event
92     Float_t b;
93     Background(b, fSpecn, fSpecp);
94
95     // Production of nuclear fragments -> num. of FREE spectators n & p
96     Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp);
97     
98     // Extract from spectators distribution the signal events:
99     // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p 
100     Mixing();
101 }
102
103 //____________________________________________________________________________
104 void 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
114 //    // Get AliRun object from file or create it if not on file
115 //    gAlice = (AliRun*)fBgrFile->Get("gAlice");
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 //    }
124     
125     gAlice->GetEvent(fNEvBgr);
126     AliHeader *header = gAlice->GetHeader();
127     AliGenEventHeader* mcHeader = header->GenEventHeader();
128     fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
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);
136 }
137
138 //____________________________________________________________________________
139 TFile* AliZDCMerger::OpenBgrFile()
140 {
141     // Initialise background event
142     TFile *file = new TFile(fFnBgr,"UPDATE");
143     printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
144     fHitsBgr = new TClonesArray("AliZDCHit",1000);
145     fMHits   = new TClonesArray("AliZDCMergedHit",1000);
146     return file;
147 }
148
149 //____________________________________________________________________________
150 void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
151                                  Int_t &fFreeSpn, Int_t &fFreeSpp)
152 {
153     //printf("\n        Fragmentation -> fSpecn = %d, fSpecp = %d\n",fSpecn,fSpecp);
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);
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);
173     printf("\n  Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
174 }
175
176 //____________________________________________________________________________
177 void AliZDCMerger::Mixing()
178 {
179     
180     //printf("\n        AliZDCMerger->Mixing\n");
181        
182     // ### Background event Hits ###########################################
183     fBgrFile->cd();
184 //    fBgrFile->ls();
185     
186    AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
187 //    if(ZDC) printf("\n        Ho trovato lo ZDC!\n");
188
189     // Hits tree
190     if(fTrHBgr) delete fTrHBgr;
191     fTrHBgr = 0;        
192     // SDigits tree
193     if(fTrSDBgr) delete fTrSDBgr;
194     fTrSDBgr = 0;    
195
196 //    fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
197     // Hits tree
198     char treeBgrName[20];
199     sprintf(treeBgrName,"TreeH%d",fNEvBgr);
200     fTrHBgr = (TTree*)gDirectory->Get(treeBgrName); 
201     if(!fTrHBgr){
202       printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
203     }    
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     
216     Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
217     //printf("\n        --- ntracks = %d\n\n", ntracks);
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;
225
226     // --- Tracks loop
227     for(itrack=0; itrack<ntracks; itrack++){
228        fTrHBgr->GetEvent(itrack);
229        
230        nhits = fHitsBgr->GetEntries();
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();
240             MHits[4] = zdcHit->GetLightPMQ();
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           
251     } // Tracks loop
252     //printf("  fNMhits (after bckg) = %d, \n",fNMhits);     
253         
254     // ### Signal event Hits ###########################################
255     // --- Neutrons
256     ExtractSignal(1);
257         
258     // --- Protons
259     ExtractSignal(2);
260     //printf("  fNMhits (after signal) = %d \n",fNMhits); 
261         
262 }
263
264 //____________________________________________________________________________
265 void AliZDCMerger::ExtractSignal(Int_t SpecType)
266 {
267
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();
275    printf("\n   --- ExtractSignal x n: file %s opened\n", fFnSpecn);
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();
282    printf("\n   --- ExtractSignal x p: file %s opened\n", fFnSpecp);
283    NumEvents = fFreeSpp;
284  }
285  //printf("\n           # of free spectators = %d\n", NumEvents);
286  //printf("\n         fNMhits (before adding signal) = %d\n",fNMhits);
287
288   TNtuple *ZDCSignal = (TNtuple*) gDirectory->Get("ZDCSignal");
289   Int_t nentries = (Int_t) ZDCSignal->GetEntries();
290   //printf("\n   # entries = %d\n", nentries);
291   
292   AliZDCMergedHit *MHit; 
293   Float_t *entry, HitsSpec[7];
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{
309      for(i=0; i<nentries; i++){  
310         ZDCSignal->GetEvent(i);
311         entry = ZDCSignal->GetArgs();
312         if(entry[0] == rnd[iev]){
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           }
317           //printf("\n   i = %d, iev = %d, entry[0] = %f, rnd[%d] = %d ",i,iev,entry[0],iev,rnd[iev]);
318           MHit = new AliZDCMergedHit(Volume, HitsSpec);
319           new((*fMHits)[fNMhits++]) AliZDCMergedHit(*MHit);
320           delete MHit;
321         }
322         else if(entry[0] > rnd[iev]){
323           iev++;
324           continue;
325         }
326      }
327   }while(iev<NumEvents);
328   
329   if(SpecType ==1){
330     //printf("\n         fNMhits (after n signal) = %d\n",fNMhits);
331     fSpecnFile->Close();
332   }
333   else if(SpecType == 2){
334     //printf("\n         fNMhits (after p signal) = %d\n",fNMhits);
335     fSpecpFile->Close();
336   }
337       
338 }
339
340 //____________________________________________________________________________
341 void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
342 {
343
344   printf("\n    AliZDCMerger->Digitize()");
345
346   AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
347 //  if(ZDC) printf("\n  Ho trovato lo ZDC!\n");
348
349   Int_t lightQ, lightC, sector[2], digit;
350   Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM1 = 0, PMZEM2 = 0;
351   Int_t i;
352   for(i=0; i<4; i++){
353      PMQZN[i] = 0;
354      PMQZP[i] = 0;
355   }
356   
357   AliZDCMergedHit *MHit;
358   Int_t imhit;
359   printf("         fNMHits = %d\n", fNMhits);
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);
366    // tmp -> c'erano i quadranti cannati!
367   if((sector[1]!=1) && (sector[1]!=2) && (sector[1]!=3) && (sector[1]!=4)){
368      printf("\n  *** ERROR!!! sector[0] = %d, sector[1] = %d\n",
369               sector[0], sector[1]);
370        sector[1] = 0;
371    }
372      lightQ    = Int_t(MHit->GetLightPMQ());
373      lightC    = Int_t(MHit->GetLightPMC());
374 //     printf("    imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d\n", 
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
390       
391   // ### Digits creation ###############################################
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);
397   printf("\n\n  ZN ###  PMCZN = %d      ADCZN = %d",PMCZN, digit);
398   PedValue = AddPedestal();
399   digit += PedValue;
400 //  printf("    PedValue = %d",PedValue);
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]);
406     printf("\n          PMQZN[%d] = %d  phe     ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit);
407     PedValue = AddPedestal();
408     digit += PedValue;
409 //    printf("  PedValue = %d",PedValue);
410     ZDC->AddDigit(sector, digit);
411   }
412   printf("\n");
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);
418   printf("\n    ZP --- PMCZP = %d       phe     ADCZP = %d ADCch",PMCZP,digit);
419   PedValue = AddPedestal();
420   digit += PedValue;
421 //  printf("    PedValue = %d",PedValue);
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]);
426     printf("\n         PMQZP[%d] = %d   phe     ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit);
427     PedValue = AddPedestal();
428     digit += PedValue;
429 //    printf("  PedValue = %d",PedValue);
430     ZDC->AddDigit(sector, digit);
431   }
432   printf("\n");
433   
434   // Create digits for ZEM
435   sector[0] = 3; 
436   sector[1] = 1; // Detector = ZEM1
437   digit  = Phe2ADCch(3, 1, PMZEM1);
438   printf("\n    ZEM *** PMZEM1 = %d      phe     ADCZEM1 = %d ADCch",PMZEM1,digit);
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);
445   printf("\n    ZEM *** PMZEM2 = %d      phe     ADCZEM2 = %d ADCch\n",PMZEM2,digit);
446   PedValue = AddPedestal();
447   digit += PedValue;
448 //  printf("  PedValue = %d\n",PedValue);
449   ZDC->AddDigit(sector, digit); 
450 }
451
452 //_____________________________________________________________________________
453 Int_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;
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 //_____________________________________________________________________________
477 Int_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.;
484   Float_t PedWidth = 5.;
485   
486   PedValue    = (Int_t) gRandom->Gaus(PedMean,PedWidth);
487   
488   return PedValue;
489 }