Important bugfix. Missing reset of the equipment header data. It was causing a wrong...
[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     fSpecn  = ((AliGenHijingEventHeader*) mcHeader)->ProjSpectatorsn();
124     fSpecp  = ((AliGenHijingEventHeader*) mcHeader)->ProjSpectatorsp();
125     // 
126     printf("\n  HIJING ev. #%d - b = %f fm, Nspecn = %d, Nspecp = %d\n",
127             fNEvBgr,fImpPar,fSpecn,fSpecp);
128 }
129
130 //____________________________________________________________________________
131 TFile* AliZDCMerger::OpenBgrFile()
132 {
133     // Initialise background event
134   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(fFnBgr);
135   if(!file)cerr<<"AliZDCMerger: background file "<<fFnBgr<<" not found\n";
136   //    TFile *file = new TFile(fFnBgr,"UPDATE");
137     printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
138     fHitsBgr = new TClonesArray("AliZDCHit",1000);
139     fMHits   = new TClonesArray("AliZDCMergedHit",1000);
140     return file;
141 }
142
143 //____________________________________________________________________________
144 void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
145                                  Int_t &fFreeSpn, Int_t &fFreeSpp)
146 {
147     //printf("\n        Fragmentation -> fSpecn = %d, fSpecp = %d\n",fSpecn,fSpecp);
148     Int_t j, zz[100], nn[100], nAlpha, ztot, ntot;
149     AliZDCFragment *frag = new AliZDCFragment(fImpPar);
150     for(j=0; j<=99; j++){
151        zz[j] =0;
152        nn[j] =0;
153     }
154
155     // Fragments generation
156     frag->GenerateIMF(zz, nAlpha);
157
158     // Attach neutrons
159     ztot=0;
160     ntot=0;
161     frag->AttachNeutrons(zz, nn, ztot, ntot);
162     fFreeSpn = fSpecn-ntot-2*nAlpha;
163     fFreeSpp = fSpecp-ztot-2*nAlpha;
164     if(fFreeSpn<0) fFreeSpn=0;
165     if(fFreeSpp<0) fFreeSpp=0;
166     //printf("\n                        2*nAlpha = %d, ztot = %d, ntot = %d\n",2*nAlpha, ztot,    ntot);
167     printf("\n  Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
168 }
169
170 //____________________________________________________________________________
171 void AliZDCMerger::Mixing()
172 {
173     
174     //printf("\n        AliZDCMerger->Mixing\n");
175        
176     // ### Background event Hits ###########################################
177     fBgrFile->cd();
178 //    fBgrFile->ls();
179     
180    AliZDC *zdc = (AliZDC *)gAlice->GetModule("ZDC");
181 //    if(zdc) printf("\n        Ho trovato lo ZDC!\n");
182
183 //    fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
184     // Hits tree
185     char treeBgrName[20];
186     sprintf(treeBgrName,"TreeH%d",fNEvBgr);
187     fTrHBgr = (TTree*)gDirectory->Get(treeBgrName); 
188     if(!fTrHBgr){
189       printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
190     }    
191 //    fTrHBgr->Print();
192     
193     // Branch address
194     TBranch *branch;
195     char branchname[20];
196     sprintf(branchname,"%s",zdc->GetName());
197     if(fTrHBgr && fHitsBgr){
198 //      printf("\n      fTrHBgr!=0 && fHitsBgr!=0\n");
199       branch = fTrHBgr->GetBranch(branchname);
200       if(branch) branch->SetAddress(&fHitsBgr);
201     }
202     
203     Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
204     //printf("\n        --- ntracks = %d\n\n", ntracks);
205     
206     Int_t itrack, nhits, ihit, j, sector[2];
207     AliZDCHit* zdcHit;
208     AliZDCMergedHit *mHit;
209     Float_t mHits[7];
210     TClonesArray &sdigits = *fMHits;    // SDigits TCArray
211     fNMhits = 0;
212
213     // --- Tracks loop
214     for(itrack=0; itrack<ntracks; itrack++){
215        fTrHBgr->GetEvent(itrack);
216        
217        nhits = fHitsBgr->GetEntries();
218 //       printf("             nhits = %d \n", nhits);
219        for(ihit=0; ihit<nhits; ihit++){    
220           zdcHit = (AliZDCHit*) fHitsBgr->UncheckedAt(ihit);
221
222             for(j=0; j<2; j++) sector[j] = zdcHit->GetVolume(j);
223             mHits[0] = zdcHit->GetPrimKinEn();
224             mHits[1] = zdcHit->GetXImpact();
225             mHits[2] = zdcHit->GetYImpact();
226             mHits[3] = zdcHit->GetSFlag();
227             mHits[4] = zdcHit->GetLightPMQ();
228             mHits[5] = zdcHit->GetLightPMC();
229             mHits[6] = zdcHit->GetEnergy();
230             mHit = new AliZDCMergedHit(sector, mHits);
231 //          mHit->Print("");
232             new((*fMHits)[fNMhits]) AliZDCMergedHit(*mHit);
233             new (sdigits[fNMhits])  AliZDCMergedHit(*mHit);
234             delete mHit;
235             fNMhits++;
236           }//Hits loop
237           
238     } // Tracks loop
239     //printf("  fNMhits (after bckg) = %d, \n",fNMhits);     
240         
241     // ### Signal event Hits ###########################################
242     // --- Neutrons
243     ExtractSignal(1);
244         
245     // --- Protons
246     ExtractSignal(2);
247     //printf("  fNMhits (after signal) = %d \n",fNMhits); 
248         
249 }
250
251 //____________________________________________________________________________
252 void AliZDCMerger::ExtractSignal(Int_t SpecType)
253 {
254
255 // printf("\n   Entering in Extract Signal\n");
256  
257  Int_t numEvents = 0;
258  if(SpecType == 1){             // --- Signal for spectator neutrons
259    fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
260    fSpecnFile = TFile::Open(fFnSpecn,"R");
261    fSpecnFile->cd();
262    printf("\n   --- ExtractSignal x n: file %s opened\n", fFnSpecn);
263    numEvents = fFreeSpn;
264  }
265  else if(SpecType == 2){        // --- Signal for spectator protons
266    fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
267    fSpecpFile = TFile::Open(fFnSpecp,"R");
268    fSpecpFile->cd();
269    printf("\n   --- ExtractSignal x p: file %s opened\n", fFnSpecp);
270    numEvents = fFreeSpp;
271  }
272  //printf("\n           # of free spectators = %d\n", numEvents);
273  //printf("\n         fNMhits (before adding signal) = %d\n",fNMhits);
274
275   TNtuple *zdcSignal = (TNtuple*) gDirectory->Get("ZDCSignal");
276   Int_t nentries = (Int_t) zdcSignal->GetEntries();
277   //printf("\n   # entries = %d\n", nentries);
278   
279   AliZDCMergedHit *mHit; 
280   Float_t *entry, hitsSpec[7];
281   Int_t pl, i, j, k, iev=0, rnd[125], volume[2];
282   for(pl=0;pl<125;pl++){
283      rnd[pl] = 0;
284   }
285   for(pl=0;pl<numEvents;pl++){
286      rnd[pl] = (Int_t) (9999*gRandom->Rndm());
287      if(rnd[pl] >= 9999) rnd[pl] = 9998;
288      //printf(" rnd[%d] = %d\n",pl,rnd[pl]);     
289   }
290   // Sorting vector in ascending order with C function QSORT 
291   qsort((void*)rnd,numEvents,sizeof(Int_t),comp);
292   //for(pl=0;pl<numEvents;pl++){
293   ////printf("  rnd[%d] = %d\n",pl,rnd[pl]);     
294   //}
295   do{
296      for(i=0; i<nentries; i++){  
297         zdcSignal->GetEvent(i);
298         entry = zdcSignal->GetArgs();
299         if(entry[0] == rnd[iev]){
300           for(k=0; k<2; k++) volume[k] = (Int_t) entry[k+1];
301           for(j=0; j<7; j++){
302              hitsSpec[j] = entry[j+3];
303           }
304           //printf("\n   i = %d, iev = %d, entry[0] = %f, rnd[%d] = %d ",i,iev,entry[0],iev,rnd[iev]);
305           mHit = new AliZDCMergedHit(volume, hitsSpec);
306           new((*fMHits)[fNMhits++]) AliZDCMergedHit(*mHit);
307           delete mHit;
308         }
309         else if(entry[0] > rnd[iev]){
310           iev++;
311           continue;
312         }
313      }
314   }while(iev<numEvents);
315   
316   if(SpecType ==1){
317     //printf("\n         fNMhits (after n signal) = %d\n",fNMhits);
318     fSpecnFile->Close();
319   }
320   else if(SpecType == 2){
321     //printf("\n         fNMhits (after p signal) = %d\n",fNMhits);
322     fSpecpFile->Close();
323   }
324       
325 }