Avoid warning.
[u/mrichter/AliRoot.git] / RICH / AliRICHMerger.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2000, 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 /* $Id$ */
17
18 #include <Riostream.h> 
19
20 #include <TTree.h> 
21 #include <TObjArray.h>
22 #include <TFile.h>
23 #include <TDirectory.h>
24 #include <TParticle.h>
25
26 #include "AliRICHMerger.h"
27 #include "AliRICHChamber.h"
28 #include "AliHitMap.h"
29 #include "AliRICHHitMapA1.h"
30 #include "AliRICH.h"
31 #include "AliRICHHit.h"
32 #include "AliRICHSDigit.h"
33 #include "AliRICHDigit.h"
34 #include "AliRICHTransientDigit.h"
35 #include "AliRun.h"
36 #include "AliPDG.h"
37
38 ClassImp(AliRICHMerger)
39
40 //___________________________________________
41   AliRICHMerger::AliRICHMerger()
42 {
43 // Default constructor    
44   fEvNrSig = 0;
45   fEvNrBgr = 0;
46   fMerge   = kDigitize;
47   fFnBgr   = 0;
48   fHitMap  = 0;
49   fList    = 0;
50   fTrH1    = 0;
51   fHitsBgr = 0;
52   fSDigitsBgr = 0;
53   fHitMap     = 0;
54   fList       = 0;
55   fBgrFile    = 0;
56 }
57
58 //------------------------------------------------------------------------
59 AliRICHMerger::~AliRICHMerger()
60 {
61 // Destructor
62   if (fTrH1)       delete fTrH1;
63   if (fHitsBgr)    delete fHitsBgr;
64   if (fSDigitsBgr) delete fSDigitsBgr;
65   if (fHitMap)     delete fHitMap;
66   if (fList)       delete fList;
67   if (fBgrFile)    delete fBgrFile;
68 }
69
70 //------------------------------------------------------------------------
71 Bool_t AliRICHMerger::Exists(const AliRICHSDigit *padhit)
72 {
73   return (fHitMap[fNch]->TestHit(padhit->PadX(),padhit->PadY()));
74 }
75
76 //------------------------------------------------------------------------
77 void AliRICHMerger::Update(AliRICHSDigit *padhit)
78 {
79   AliRICHTransientDigit *pdigit = 
80     static_cast<AliRICHTransientDigit*>(
81       fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
82
83   // update charge
84   //
85   Int_t iqpad    = Int_t(padhit->QPad()); // charge per pad
86   pdigit->AddSignal(iqpad);
87   pdigit->AddPhysicsSignal(iqpad);              
88
89   // update list of tracks
90   //
91   Int_t track, charge;    
92   if (fSignal) {
93     track = fTrack;
94     charge = iqpad;
95   } else {
96     track = kBgTag;
97     charge = kBgTag;
98   }
99   pdigit->UpdateTrackList(track,charge);
100 }
101
102 //------------------------------------------------------------------------
103 void AliRICHMerger::CreateNew(AliRICHSDigit *padhit)
104 {
105   fList->AddAtAndExpand(
106     new AliRICHTransientDigit(fNch,fDigits),fCounter);
107   fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
108
109   AliRICHTransientDigit* pdigit = 
110     static_cast<AliRICHTransientDigit*>(fList->Last());
111   // list of tracks
112   Int_t track, charge;
113   if (fSignal) {
114     track = fTrack;
115     charge = padhit->QPad();
116   } else {
117     track = kBgTag;
118     charge = kBgTag;
119   }
120   pdigit->AddToTrackList(track,charge);
121   fCounter++;
122 }
123
124
125 //------------------------------------------------------------------------
126 void AliRICHMerger::Init()
127 {
128 // Initialisation
129     
130   if (fMerge && !fBgrFile) fBgrFile = InitBgr();
131 }
132
133
134
135 //------------------------------------------------------------------------
136 TFile* AliRICHMerger::InitBgr()
137 {
138 // Initialise background event
139   TFile *file = new TFile(fFnBgr);
140 // add error checking later
141   printf("\n AliRICHMerger has opened %s file with background event \n", fFnBgr);
142   fHitsBgr     = new TClonesArray("AliRICHHit",1000);
143   fSDigitsBgr  = new TClonesArray("AliRICHSDigit",1000);
144   return file;
145 }
146
147 //------------------------------------------------------------------------
148 void AliRICHMerger::Digitise(Int_t nev, Int_t flag)
149 {
150
151   // keep galice.root for signal and name differently the file for 
152   // background when add! otherwise the track info for signal will be lost !
153
154   Int_t particle;
155
156   AliRICHChamber*       iChamber;
157   AliSegmentation*  segmentation;
158   
159   fList = new TObjArray;
160     
161   AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
162   fHitMap= new AliHitMap* [kNCH];
163     
164   if (fMerge ) {
165     fBgrFile->cd();
166     //
167     // Get Hits Tree header from file
168     if(fHitsBgr) fHitsBgr->Clear();
169     if(fSDigitsBgr) fSDigitsBgr->Clear();
170     if(fTrH1) delete fTrH1;
171     fTrH1 = 0;
172       
173     char treeName[20];
174     sprintf(treeName,"TreeH%d",fEvNrBgr);
175     fTrH1 = (TTree*)gDirectory->Get(treeName);
176     if (!fTrH1) {
177       printf("ERROR: cannot find Hits Tree for event:%d\n",fEvNrBgr);
178     }
179     //
180     // Set branch addresses
181     TBranch *branch;
182     char branchname[20];
183     sprintf(branchname,"%s",pRICH->GetName());
184     if (fTrH1 && fHitsBgr) {
185       branch = fTrH1->GetBranch(branchname);
186       if (branch) branch->SetAddress(&fHitsBgr);
187     }
188     if (fTrH1 && fSDigitsBgr) {
189       branch = fTrH1->GetBranch("RICHSDigits");
190       if (branch) branch->SetAddress(&fSDigitsBgr);
191     }
192   }
193     
194   for (Int_t i =0; i<kNCH; i++) {
195     iChamber= &(pRICH->Chamber(i));
196     segmentation=iChamber->GetSegmentationModel(1);
197     fHitMap[i] = new AliRICHHitMapA1(segmentation, fList);
198   }
199   //
200   //   Loop over tracks
201   //
202     
203   fSignal = kTRUE;
204   fCounter = 0;
205   TTree *treeH = pRICH->TreeH();
206   Int_t ntracks =(Int_t) treeH->GetEntries();
207   for (fTrack=0; fTrack<ntracks; fTrack++) {
208     gAlice->ResetHits();
209     treeH->GetEvent(fTrack);
210     //
211     //   Loop over hits
212     for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1); 
213         mHit;
214         mHit=(AliRICHHit*)pRICH->NextHit()) 
215     {
216           
217       fNch = mHit->Chamber()-1;  // chamber number
218       Int_t trackIndex = mHit->Track();
219       if (fNch >= kNCH) {
220         cerr<<"AliRICHMerger: chamber nr. fNch out of range: "<<fNch<<endl;
221         cerr<<"               track: "<<fTrack<<endl;
222         continue;
223       }
224       iChamber = &(pRICH->Chamber(fNch));
225           
226
227 // 
228 // If flag is set, create digits only for some particles
229 //
230       TParticle *current = (TParticle*)gAlice->Particle(trackIndex);
231       if (current->GetPdgCode() >= 50000050)
232       {
233         TParticle *motherofcurrent = (TParticle*)gAlice->Particle(current->GetFirstMother());
234         particle = motherofcurrent->GetPdgCode();
235       }
236       else
237       {
238         particle = current->GetPdgCode();
239       }
240
241           
242       //printf("Flag:%d\n",flag);
243       //printf("Track:%d\n",track);
244       //printf("Particle:%d\n",particle);
245           
246       Int_t digitise=1;
247           
248       if (flag == 1) 
249         if(TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
250           digitise=0;
251           
252       if (flag == 2)
253         if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 
254            || TMath::Abs(particle)==311)
255           digitise=0;
256           
257       if (flag == 3 && TMath::Abs(particle)==2212)
258         digitise=0;
259           
260       if (flag == 4 && TMath::Abs(particle)==13)
261         digitise=0;
262           
263       if (flag == 5 && TMath::Abs(particle)==11)
264         digitise=0;
265           
266       if (flag == 6 && TMath::Abs(particle)==2112)
267         digitise=0;
268           
269           
270       //printf ("Particle: %d, Flag: %d, Digitise: %d\n",particle,flag,digitise); 
271           
272           
273       if (digitise)
274       {
275               
276         //
277         // Loop over pad hits
278         for (AliRICHSDigit* mPad=
279                (AliRICHSDigit*)pRICH->FirstPad(mHit,pRICH->SDigits());
280              mPad;
281              mPad=(AliRICHSDigit*)pRICH->NextPad(pRICH->SDigits()))
282         {
283           Int_t iqpad    = mPad->QPad();       // charge per pad
284           fDigits[0]=mPad->PadX();
285           fDigits[1]=mPad->PadY();
286           fDigits[2]=iqpad;
287           fDigits[3]=iqpad;
288           fDigits[4]=mPad->HitNumber();
289                   
290           // build the list of fired pads and update the info
291           if (Exists(mPad)) {
292             Update(mPad);
293           } else {
294             CreateNew(mPad);
295           }
296         } //end loop over clusters
297       }// track type condition
298     } // hit loop
299   } // track loop
300     
301   // open the file with background
302     
303   if (fMerge) {
304     fSignal = kFALSE;
305     ntracks = (Int_t)fTrH1->GetEntries();
306     //
307     //   Loop over tracks
308     //
309     for (fTrack = 0; fTrack < ntracks; fTrack++) {
310         
311       if (fHitsBgr)       fHitsBgr->Clear();
312       if (fSDigitsBgr)    fSDigitsBgr->Clear();
313         
314       fTrH1->GetEvent(fTrack);
315       //
316       //   Loop over hits
317       AliRICHHit* mHit;
318       for(Int_t i = 0; i < fHitsBgr->GetEntriesFast(); ++i) 
319       { 
320         mHit   = (AliRICHHit*) (*fHitsBgr)[i];
321         fNch   = mHit->Chamber()-1;  // chamber number
322         iChamber = &(pRICH->Chamber(fNch));
323
324         //
325         // Loop over pad hits
326         for (AliRICHSDigit* mPad =
327                (AliRICHSDigit*)pRICH->FirstPad(mHit,fSDigitsBgr);
328              mPad;
329              mPad = (AliRICHSDigit*)pRICH->NextPad(fSDigitsBgr))
330         {
331           fDigits[0] = mPad->PadX(); 
332           fDigits[1] = mPad->PadY();  
333           fDigits[2] = mPad->QPad();
334           fDigits[3] = 0;
335           fDigits[4] = -1; // set hit number to -1 for bgr
336                 
337           // build the list of fired pads and update the info
338           if (!Exists(mPad)) {
339             CreateNew(mPad);
340           } else {
341             Update(mPad);
342           } //  end if !Exists
343         } //end loop over clusters
344       } // hit loop
345     } // track loop
346       
347     TTree *treeK = gAlice->TreeK();
348     TFile *file = NULL;
349       
350     if (treeK) file = treeK->GetCurrentFile();
351     file->cd();
352   } // if fMerge
353     
354   Int_t tracks[kMAXTRACKSPERRICHDIGIT];
355   Int_t charges[kMAXTRACKSPERRICHDIGIT];
356   Int_t nentries=fList->GetEntriesFast();
357     
358   for (Int_t nent=0;nent<nentries;nent++) {
359     AliRICHTransientDigit *transDigit=(AliRICHTransientDigit*)fList->At(nent);
360     if (transDigit==0) continue; 
361     Int_t ich=transDigit->GetChamber();
362     Int_t q=transDigit->Signal(); 
363     iChamber=&(pRICH->Chamber(ich));
364     AliRICHResponse * response=iChamber->GetResponseModel();
365     Int_t adcmax= (Int_t) response->MaxAdc();
366       
367       
368     // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
369     //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY()));
370     Int_t pedestal = iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY());
371
372     //printf("Pedestal:%d\n",pedestal);
373     //Int_t pedestal=0;
374     Float_t treshold = (pedestal + 4*2.2);
375       
376     Float_t meanNoise = gRandom->Gaus(2.2, 0.3);
377     Float_t noise     = gRandom->Gaus(0, meanNoise);
378     q+=(Int_t)(noise + pedestal);
379     //q+=(Int_t)(noise);
380     //          magic number to be parametrised !!! 
381     if ( q <= treshold) 
382     {
383       q = q - pedestal;
384       continue;
385     }
386     q = q - pedestal;
387     if ( q >= adcmax) q=adcmax;
388     fDigits[0]=transDigit->PadX();
389     fDigits[1]=transDigit->PadY();
390     fDigits[2]=q;
391     fDigits[3]=transDigit->Physics();
392     fDigits[4]=transDigit->Hit();
393
394     Int_t nptracks = transDigit->GetNTracks();  
395
396     // this was changed to accomodate the real number of tracks
397     if (nptracks > kMAXTRACKSPERRICHDIGIT) {
398       printf("Attention - tracks > 10 %d\n",nptracks);
399       nptracks=kMAXTRACKSPERRICHDIGIT;
400     }
401     if (nptracks > 2) {
402       printf("Attention - tracks > 2  %d \n",nptracks);
403       //printf("cat,ich,ix,iy,q %d %d %d %d %d \n",
404       //icat,ich,digits[0],digits[1],q);
405     }
406     for (Int_t tr=0;tr<nptracks;tr++) {
407       tracks[tr]=transDigit->GetTrack(tr);
408       charges[tr]=transDigit->GetCharge(tr);
409     }      //end loop over list of tracks for one pad
410     if (nptracks < kMAXTRACKSPERRICHDIGIT ) {
411       for (Int_t t=nptracks; t<kMAXTRACKSPERRICHDIGIT; t++) {
412         tracks[t]=0;
413         charges[t]=0;
414       }
415     }
416     //write file
417     //if (ich==2)
418     //fprintf(points,"%4d,      %4d,      %4d\n",digits[0],digits[1],digits[2]);
419       
420     // fill digits
421     pRICH->AddDigits(ich,tracks,charges,fDigits);
422   }     
423   gAlice->TreeD()->Fill();
424
425   fList->Delete();
426   for(Int_t ii=0;ii<kNCH;++ii) {
427     if (fHitMap[ii]) {
428       delete fHitMap[ii];
429       fHitMap[ii]=0;
430     }
431   }
432     
433   TClonesArray *richDigits;
434   for (Int_t k=0;k<kNCH;k++) {
435     richDigits = pRICH->DigitsAddress(k);
436     Int_t ndigit=richDigits->GetEntriesFast();
437     printf ("Chamber %d digits %d \n",k,ndigit);
438   }
439   pRICH->ResetDigits(); /// ??? should it be here???
440   gAlice->TreeD()->Write(0,TObject::kOverwrite);
441   // reset tree
442   //    gAlice->TreeD()->Reset();
443   //    delete fList; // deleted in dtor
444   // gObjectTable->Print();
445
446 }