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