Some function moved to AliZDC
[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 <TH1.h>
29
30 // --- AliRoot header files
31 #include "AliZDCMerger.h"
32 #include "AliZDC.h"
33 #include "AliZDCHit.h"
34 #include "AliZDCDigit.h"
35 #include "AliZDCFragment.h"
36 #include "AliRun.h"
37 #include "AliDetector.h"
38 #include "AliHeader.h"
39 #include "AliGenEventHeader.h"
40 #include "/home/oppedisa/AliRoot/pro/EVGEN/AliGenHijingEventHeader.h"
41
42 ClassImp(AliZDCMerger)
43
44 //____________________________________________________________________________
45 AliZDCMerger::AliZDCMerger()
46 {
47 // Default constructor    
48     fMerge       = kDigitize;
49     fNEvBgr      = 0;
50     fFnBgr       = 0;
51 //    fFnBgr       = "/home/oppedisa/alirun/hijing.root";
52     fBgrFile     = 0;
53     fTrHBgr      = 0;
54     fTrSDBgr     = 0;
55     fImpPar      = 0;
56     fSpecn       = 0;
57     fSpecp       = 0;
58     fFreeSpn     = 0;
59     fFreeSpp     = 0;
60     fFnSpecn     = 0;
61 //    fFnSpecn     = "/home/oppedisa/alirun/ZNsignal.root";
62     fSpecnFile   = 0;
63     fFnSpecp     = 0;
64 //    fFnSpecp     = "/home/oppedisa/alirun/ZPsignal.root";
65     fSpecpFile   = 0;
66     
67 }
68
69 //____________________________________________________________________________
70 AliZDCMerger::~AliZDCMerger()
71 {
72 // Destructor
73     if (fBgrFile)    delete fBgrFile;
74     if (fTrHBgr)     delete fTrHBgr;
75     if (fTrSDBgr)    delete fTrSDBgr;
76     if (fHitsBgr)    delete fHitsBgr;
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     Int_t nspecn, nspecp;
88     Background(b, nspecn, nspecp);
89
90     // Production of nuclear fragments -> num. of FREE spectators n & p
91     Fragmentation(b, nspecn, nspecp, fFreeSpn, fFreeSpp);
92     
93     // Extract from spectators distribution the signal events:
94     // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p 
95     Mixing(fFreeSpn, fFreeSpp);
96 }
97
98 //____________________________________________________________________________
99 void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
100 {
101     
102     // --- Open the background file
103     if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
104     
105     // --- Read from the TreeE impact parameter (b),
106     //     # of spectators n and p (fSpecn, fSpecp)
107     fBgrFile->cd();
108
109     // Get AliRun object from file or create it if not on file
110     gAlice = (AliRun*)fBgrFile->Get("gAlice");
111     if (!gAlice) {
112       gAlice = (AliRun*)fBgrFile->Get("gAlice");
113       if (gAlice) printf("AliRun object found on file\n");
114       if (!gAlice) {
115           printf("\n create new gAlice object");
116           gAlice = new AliRun("gAlice","Alice test program");
117       }
118     }
119
120 //    gAlice = (AliRun*)fBgrFile->Get("gAlice");
121     AliHeader *header = gAlice->GetHeader();
122     AliGenEventHeader* mcHeader = header->GenEventHeader();
123     fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
124     fSpecn  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
125     fSpecp  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
126 //      printf("\n      HIJING simulation - b = %f fm, Nspecn = %d, Nspecp = %d\n",fImpPar,fSpecn,fSpecp);
127 }
128
129 //____________________________________________________________________________
130 TFile* AliZDCMerger::OpenBgrFile()
131 {
132     // Initialise background event
133     TFile *file = new TFile(fFnBgr,"UPDATE");
134 //      printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
135     fHitsBgr = new TClonesArray("AliZDCHit",1000);
136     return file;
137 }
138
139 //____________________________________________________________________________
140 void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
141                                  Int_t &fFreeSpn, Int_t &fFreeSpp)
142 {
143     Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot;
144     AliZDCFragment *frag = new AliZDCFragment(fImpPar);
145     for(j=0; j<=99; j++){
146        zz[j] =0;
147        nn[j] =0;
148     }
149
150     // Fragments generation
151     frag->GenerateIMF(zz, nAlpha);
152
153     // Attach neutrons
154     Ztot=0;
155     Ntot=0;
156     frag->AttachNeutrons(zz, nn, Ztot, Ntot);
157     fFreeSpn = fSpecn-Ztot-2*nAlpha;
158     fFreeSpp = fSpecp-Ntot-2*nAlpha;
159 //      printf("\n      Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
160 }
161
162 //____________________________________________________________________________
163 void AliZDCMerger::Mixing(Int_t fFreeSpn, Int_t fFreeSpp)
164 {
165
166 //      printf("\n      AliZDCMerger->Mixing\n");
167         
168     // ### Background event Hits ###########################################
169     fBgrFile->cd();
170 //    fBgrFile->ls();
171     
172     AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
173 //    if(ZDC) printf("\n        Ho trovato lo ZDC!\n");
174
175     // Hits tree
176     if(fTrHBgr) delete fTrHBgr;
177     fTrHBgr = 0;    
178     // SDigits tree
179     if(fTrSDBgr) delete fTrSDBgr;
180     fTrSDBgr = 0;    
181
182 //    // Read from Event Tree the # of events of the Hijing file 
183 //    TTree *thead = (TTree *) fBgrFile->Get("TE");
184 //    if(!te){
185 //      printf("\n      ERROR! -> File %s does not contain TreeE\n");
186 //    }
187 //    fNEvBgr = (Int_t) te->GetEntries();
188     fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
189
190     // Hits Tree
191     char treeBgrName[20];
192     sprintf(treeBgrName,"TreeH%d",fNEvBgr);
193     fTrHBgr = (TTree*)gDirectory->Get(treeBgrName); // TreeH
194     if(!fTrHBgr){
195       printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
196     }    
197     Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
198 //    printf("\n        --- ntracks = %d\n",ntracks);
199
200     // SDigits Tree
201     char treeSBgrName[20];
202     sprintf(treeSBgrName,"TreeS%d",fNEvBgr);
203     fTrSDBgr = (TTree*)gDirectory->Get(treeSBgrName); // TreeH
204     if(!fTrSDBgr){
205       printf("\n ERROR -> Can't find TreeS%d in background file\n",fNEvBgr);
206     }    
207
208     Int_t itrack, i, volume[2], detector, quadrant;
209     TClonesArray &sdigits = *fHitsBgr;  // SDigits TCArray
210     Float_t hits[10];
211     // --- Tracks loop
212     for(itrack=0; itrack<ntracks; itrack++){
213 //       printf("               itrack = %d \n", itrack);
214        fTrHBgr->GetEvent(itrack);
215        
216        fTrack = itrack;
217        Int_t NMhits  = 0;
218        for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1); 
219            zdcHit;
220            zdcHit = (AliZDCHit*)ZDC->NextHit()){ 
221
222           for(i=0; i<2; i++) volume[i] = zdcHit->GetVolume(i);
223 //        printf("\n            volume[0] = %d volume[1]= %d \n",volume[0], volume[1]);
224           if(volume[0] != 3){
225             // Background hits  
226             hits[7] = zdcHit->GetLightPMQ();
227             hits[8] = zdcHit->GetLightPMC();
228 //          printf("\n  Prima ### Background -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
229
230             //Signal hits
231             detector = volume[0];
232             quadrant = volume[1];
233             ExtractSignal(detector, quadrant, fQuadLight, fComLight);      
234             hits[7] += fQuadLight;
235             hits[8] += fComLight;
236 //          printf("\n  Bckg + Signal -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
237
238             zdcHit->SetLightPMQ(hits[7]);
239             zdcHit->SetLightPMC(hits[8]);
240             hits[7] = zdcHit->GetLightPMQ();
241             hits[8] = zdcHit->GetLightPMC();
242 //          printf("\n  Dopo ### Background -> PMQ = %f, PMC = %f\n",hits[7],hits[8]);
243                     
244             new (sdigits[NMhits++]) AliZDCHit(zdcHit);
245 //            printf("\n        NMhits = %d\n",NMhits);
246           }
247           fBgrFile->cd();
248           fTrSDBgr->Fill();
249           fTrSDBgr->Write(0,TObject::kOverwrite);
250        } //Hits loop
251        Digitize();
252     } // Tracks loop
253         
254 //    fBgrFile->Close();
255         
256 }
257
258 //____________________________________________________________________________
259 void AliZDCMerger::ExtractSignal(Int_t detector, Int_t quadrant,
260                                     Float_t &fQuadLight, Float_t &fComLight)
261 {
262
263 //    printf("\n        Entering ExtractSignal method -> detector = %d quadrant = %d\n",
264 //           detector, quadrant);
265
266     // Connect spectator n histo's file
267     fSpecnFile = TFile::Open("/home/oppedisa/alirun/ZNsignal.root");
268     fSpecnFile->cd();
269 //    fSpecnFile->ls();
270     TH1F *hPMQ1zn = (TH1F*) gDirectory->Get("hPMQ1zn;1");
271     TH1F *hPMQ2zn = (TH1F*) gDirectory->Get("hPMQ2zn;1");
272     TH1F *hPMQ3zn = (TH1F*) gDirectory->Get("hPMQ3zn;1");
273     TH1F *hPMQ4zn = (TH1F*) gDirectory->Get("hPMQ4zn;1");
274     TH1F *hPMC1zn = (TH1F*) gDirectory->Get("hPMC1zn;1");
275     TH1F *hPMC2zn = (TH1F*) gDirectory->Get("hPMC2zn;1");
276     TH1F *hPMC3zn = (TH1F*) gDirectory->Get("hPMC3zn;1");
277     TH1F *hPMC4zn = (TH1F*) gDirectory->Get("hPMC4zn;1");
278 //    Axis_t x = hPMQ1zn -> GetRandom();
279 //    printf("  hPMQ1zn -> GetRandom() = %f\n",x);
280     
281     // Connect spectator p histo's file
282     fSpecpFile = TFile::Open("/home/oppedisa/alirun/ZPsignal.root");
283     fSpecpFile->cd(); 
284 //    fSpecpFile->ls();
285     TH1F *hPMQ1zp = (TH1F*) gDirectory->Get("hPMQ1zp;1");
286     TH1F *hPMQ2zp = (TH1F*) gDirectory->Get("hPMQ2zp;1");
287     TH1F *hPMQ3zp = (TH1F*) gDirectory->Get("hPMQ3zp;1");
288     TH1F *hPMQ4zp = (TH1F*) gDirectory->Get("hPMQ4zp;1");
289     TH1F *hPMC1zp = (TH1F*) gDirectory->Get("hPMC1zp;1");
290     TH1F *hPMC2zp = (TH1F*) gDirectory->Get("hPMC2zp;1");
291     TH1F *hPMC3zp = (TH1F*) gDirectory->Get("hPMC3zp;1");
292     TH1F *hPMC4zp = (TH1F*) gDirectory->Get("hPMC4zp;1");
293
294     if(detector == 1){ // --- ZN
295       if(quadrant == 1){
296         fQuadLight = (Float_t) hPMQ1zn -> GetRandom();
297         fComLight = (Float_t) hPMC1zn -> GetRandom();
298       }
299       else if(quadrant == 2){
300         fQuadLight = (Float_t) hPMQ2zn -> GetRandom();
301         fComLight = (Float_t) hPMC2zn -> GetRandom();
302       }
303       else if(quadrant == 3){
304         fQuadLight = (Float_t) hPMQ3zn -> GetRandom();
305         fComLight = (Float_t) hPMC3zn -> GetRandom();
306       }
307       else if(quadrant == 4){
308         fQuadLight = (Float_t) hPMQ4zn -> GetRandom();
309         fComLight = (Float_t) hPMC4zn -> GetRandom();
310       }
311     }
312
313     else if(detector == 2){ // --- ZP
314       fSpecpFile->cd(); // Connect spectator p histo's file
315 //      fSpecpFile->ls();
316       if(quadrant == 1){
317         fQuadLight = (Float_t) hPMQ1zp -> GetRandom();
318         fComLight = (Float_t) hPMC1zp -> GetRandom();
319       }
320       else if(quadrant == 2){
321         fQuadLight = (Float_t) hPMQ2zp -> GetRandom();
322         fComLight = (Float_t) hPMC2zp -> GetRandom();
323       }
324       else if(quadrant == 3){
325         fQuadLight = (Float_t) hPMQ3zp -> GetRandom();
326         fComLight = (Float_t) hPMC3zp -> GetRandom();
327       }
328       else if(quadrant == 4){
329         fQuadLight = (Float_t) hPMQ4zp -> GetRandom();
330         fComLight = (Float_t) hPMC4zp -> GetRandom();
331       }
332     }
333 //    printf("  ---     Exiting ExtractSignal -> fQuadLight = %f, fComLight = %f",
334 //           fQuadLight,fComLight);
335
336 }
337
338 //____________________________________________________________________________
339 void AliZDCMerger::Digitize()
340 {
341
342 //  printf("\n  AliZDCMerger->Digitize()");
343
344   AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
345 //  if(ZDC) printf("\n  Ho trovato lo ZDC!\n");
346
347   Int_t itrack, lightQ, lightC, sector[2], digit;
348   Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM = 0;
349   Int_t i;
350   for(i=0; i<4; i++){
351      PMQZN[i] = 0;
352      PMQZP[i] = 0;
353   }
354
355   // ### Digitization "on the flight" (no merging) #######################
356   if(fMerge == kDigitize){
357 //    printf("\n                fMerge == kDigitize\n");
358     fTrHBgr = gAlice->TreeH(); // TTree
359     TTree *treeH = gAlice->TreeH();
360     if(!fTrHBgr){
361       printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
362     }    
363     Int_t ntracks =  (Int_t) treeH->GetEntries();
364 //    printf("\n        --- ntracks = %d\n",ntracks);
365         
366     // Loop over tracks
367     for(itrack=0; itrack<ntracks; itrack++){
368        gAlice->ResetHits();
369        gAlice->TreeH()->GetEvent(itrack);
370        
371        // Loop over hits
372        for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1); 
373            zdcHit;
374            zdcHit = (AliZDCHit*)ZDC->NextHit()){ 
375           sector[0] = zdcHit->GetVolume(0);
376           sector[1] = zdcHit->GetVolume(1);
377           lightQ    = Int_t(zdcHit->GetLightPMQ());
378           lightC    = Int_t(zdcHit->GetLightPMC());
379 //        printf("\n    Digitise -> DET. = %d, quad = %d", sector[0], sector[1]);
380 //        printf("                  PMQ = %d, PMC = %d",  lightQ, lightC);
381           
382           if(sector[0] == 1){   //ZN 
383             PMCZN = PMCZN + lightC;
384             PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
385           }
386           else if(sector[0] == 2){      //ZP 
387             PMCZP = PMCZP + lightC;
388             PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
389           }
390           else if(sector[0] == 3){      //ZEM 
391             PMZEM = PMZEM + lightC;
392           }
393        } // hits loop
394     } // tracks loop
395   } // if(fMerge)
396
397   // ### Merging and digitization #######################################
398   else if(fMerge == kMerge){
399 //    printf("\n                fMerge == kMerge\n");
400
401 //    Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
402 //    printf("\n        --- ntracks = %d\n",ntracks);
403 //        
404 //    // Loop over tracks
405 //    for(itrack=0; itrack<ntracks; itrack++){
406        fTrHBgr->GetEvent(fTrack);
407        
408 //       printf("\n\n   Track # %d --- Digitise -> ", fTrack);
409        // Loop over hits
410        for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1); 
411            zdcHit;
412            zdcHit = (AliZDCHit*)ZDC->NextHit()){ 
413           sector[0] = zdcHit->GetVolume(0);
414           sector[1] = zdcHit->GetVolume(1);
415           lightQ    = Int_t(zdcHit->GetLightPMQ());
416           lightC    = Int_t(zdcHit->GetLightPMC());
417 //        printf("\n    DET. = %d, quad = %d,PMQ = %d, PMC = %d", 
418 //               sector[0], sector[1],lightQ, lightC);
419           
420           if(sector[0] == 1){   //ZN 
421             PMCZN = PMCZN + lightC;
422             PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
423           }
424           else if(sector[0] == 2){      //ZP 
425             PMCZP = PMCZP + lightC;
426             PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
427           }
428           else if(sector[0] == 3){      //ZEM 
429             PMZEM = PMZEM + lightC;
430           }
431        } // hits loop
432 //    } // tracks loop
433   } // if(fMerge)
434
435       
436   // Create digits for ZN
437   Int_t PedValue;
438   sector[0] = 1; // Detector = ZN
439   sector[1] = 0; // Common PM ADC
440   digit = Phe2ADCch(1, 0, PMCZN);
441 //  printf("\n\n        ZN ###  PMCZN = %d      ADCZN = %d",PMCZN, digit);
442   ZDC->AddDigit(sector, digit);
443   Int_t j;
444   for(j=0; j<4; j++){
445     sector[1] = j+1; // Towers PM ADCs
446     digit = Phe2ADCch(1, j+1, PMQZN[j]);
447 //    printf("\n                PMQZN[%d] = %d  ADCZN[%d] = %d",j,PMQZN[j],j,digit);
448     PedValue = AddPedestal();
449     digit += PedValue;
450 //    printf("  PedValue = %d",PedValue);
451     ZDC->AddDigit(sector, digit);
452   }
453 //    printf("\n");
454   
455   // Create digits for ZP
456   sector[0] = 2; // Detector = ZP
457   sector[1] = 0; // Common PM ADC
458   digit = Phe2ADCch(2, 0, PMCZP);
459 //  printf("\n  ZP --- PMCZP = %d       ADCZP = %d",PMCZP,digit);
460   ZDC->AddDigit(sector, digit);
461   for(j=0; j<4; j++){
462     sector[1] = j+1; // Towers PM ADCs
463     digit = Phe2ADCch(2, j+1, PMQZP[j]);
464 //    printf("\n               PMQZP[%d] = %d   ADCZP[%d] = %d",j,PMQZP[j],j,digit);
465     PedValue = AddPedestal();
466     digit += PedValue;
467 //    printf("  PedValue = %d",PedValue);
468     ZDC->AddDigit(sector, digit);
469   }
470 //    printf("\n");
471   
472   // Create digits for ZEM
473   sector[0] = 3; // Detector = ZZEM
474   sector[1] = 0; // Single PM ADC
475   digit  = Phe2ADCch(3, 0, PMZEM);
476 //  printf("\n  ZEM *** PMZEM = %d      ADCZEM = %d",PMZEM,digit);
477     PedValue = AddPedestal();
478     digit += PedValue;
479 //    printf("  PedValue = %d\n",PedValue);
480   ZDC->AddDigit(sector, digit); 
481
482 }
483
484 //_____________________________________________________________________________
485 Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light)
486 {
487   // Evaluation of the ADC channel corresponding to the light yield Light
488
489   if(gAlice->GetDebug() > 0){
490 //    printf("\n  Phe2ADCch -> Detector = %d, Quadrant = %d, Light = %d\n", Det, Quad, Light);
491   }
492   
493   Int_t ADCch = 0;
494   // Parameters for conversion of light yield in ADC channels
495   Float_t fPMGain[3][5];      // PM gain
496   Float_t fADCRes;            // ADC conversion factor
497   
498   Int_t j,i;
499   for(i=0; i<3; i++){
500      for(j=0; j<5; j++){
501         fPMGain[i][j]   = 100000.;
502      }
503   }
504   fADCRes   = 0.00000064; // ADC Resolution: 250 fC/ADCch
505   
506   ADCch = (Int_t) (Light*fPMGain[Det-1][Quad]*fADCRes);
507      
508   return ADCch;
509 }
510
511 //_____________________________________________________________________________
512 Int_t AliZDCMerger::AddPedestal()
513 {
514   // --- Pedestal value -> extracted from a gaussian distribution
515   // obtained from the beam test on the ZEM prototype (Aug. 2000)
516   
517   Int_t PedValue;
518   Float_t PedMean  = 50.;
519   Float_t PedWidth = 10.;
520   
521   PedValue    = (Int_t) gRandom->Gaus(PedMean,PedWidth);
522   
523   return PedValue;
524 }