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