PatRec adapted to new IO
[u/mrichter/AliRoot.git] / RICH / AliRICHDigitizer.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
17 //Piotr.Skowronski@cern.ch :
18 //Corrections applied in order to compile (only) with new I/O and folder structure
19 //To be implemented correctly by responsible
20
21 #include <Riostream.h> 
22
23 #include <TTree.h> 
24 #include <TObjArray.h>
25 #include <TFile.h>
26 #include <TDirectory.h>
27 #include <TParticle.h>
28
29 #include "AliRunLoader.h"
30 #include "AliLoader.h"
31
32 #include "AliRICHDigitizer.h"
33 #include "AliRICHChamber.h"
34 #include "AliHitMap.h"
35 #include "AliRICHHitMapA1.h"
36 #include "AliRICH.h"
37 #include "AliRICHSDigit.h"
38 #include "AliRICHDigit.h"
39 #include "AliRICHTransientDigit.h"
40 #include "AliRun.h"
41 #include "AliPDG.h"
42 #include "AliRunDigitizer.h"
43
44 ClassImp(AliRICHDigitizer)
45
46 //__________________________________________________________________________________________________
47 AliRICHDigitizer::AliRICHDigitizer()
48 {//default ctor - don't use it   
49   fHits = 0;
50   fSDigits = 0;
51   fHitMap = 0;
52   fTDList = 0;
53 }//default ctor
54 //__________________________________________________________________________________________________
55 AliRICHDigitizer::AliRICHDigitizer(AliRunDigitizer *pManager) 
56                  :AliDigitizer(pManager)
57 {//main ctor which should be used
58   if(fManager->GetDebug())Info("main ctor","Start.");
59   fHits = 0;
60   fSDigits = 0;
61   fHitMap = 0;
62   fTDList = 0;
63   fDebug  = pManager->GetDebug(); 
64 }//main ctor
65 //__________________________________________________________________________________________________
66 AliRICHDigitizer::~AliRICHDigitizer()
67 {//dtor
68   if(fManager->GetDebug())Info("dtor","Start.");
69   
70   if(fHits)    {fHits->Delete();   delete fHits;}
71   if(fSDigits) {fSDigits->Delete();delete fSDigits;}
72   if(fTDList)  {fTDList->Delete(); delete fTDList;}
73   for(Int_t i=0; i<kNCH; i++) delete fHitMap[i];    delete [] fHitMap;
74 }//dtor
75 //__________________________________________________________________________________________________
76 Bool_t AliRICHDigitizer::Exists(const AliRICHSDigit *p) 
77 {
78   return (fHitMap[fNch]->TestHit(p->PadX(),p->PadY()));
79 }//Exists
80 //__________________________________________________________________________________________________
81 void AliRICHDigitizer::Update(AliRICHSDigit *padhit)
82 {
83   AliRICHTransientDigit *pdigit = 
84     static_cast<AliRICHTransientDigit*>(
85       fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
86
87   // update charge
88   //
89   Int_t iqpad    = Int_t(padhit->QPad()); // charge per pad
90   pdigit->AddSignal(iqpad);
91   pdigit->AddPhysicsSignal(iqpad);            
92
93   // update list of tracks
94   //
95   Int_t track, charge;
96   track = fTrack+fMask;
97   if (fSignal) {
98     charge = iqpad;
99   } else {
100     charge = kBgTag;
101   }
102   pdigit->UpdateTrackList(track,charge);
103 }//Update()
104 //__________________________________________________________________________________________________
105 void AliRICHDigitizer::CreateNew(AliRICHSDigit *padhit)
106 {
107   fTDList->AddAtAndExpand(
108     new AliRICHTransientDigit(fNch,fDigits),fCounter);
109   fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
110
111   AliRICHTransientDigit* pdigit = 
112     static_cast<AliRICHTransientDigit*>(fTDList->Last());
113   // list of tracks
114   Int_t track, charge;
115   if (fSignal) {
116     track = fTrack;
117     charge = padhit->QPad();
118   } else {
119     track = kBgTag;
120     charge = kBgTag;
121   }
122   pdigit->AddToTrackList(track,charge);
123   fCounter++;
124 }//CreateNew()
125 //__________________________________________________________________________________________________
126 Bool_t AliRICHDigitizer::Init()
127 {// Initialisation
128   if(GetDebug())Info("Init","Start.");
129   fHits     = new TClonesArray("AliRICHhit",1000);
130   fSDigits  = new TClonesArray("AliRICHSDigit",1000);
131   return kTRUE;
132 }//Init()
133 //__________________________________________________________________________________________________
134 void AliRICHDigitizer::Exec(Option_t* option)
135 {
136   if(GetDebug())Info("Exec","Start with option=%s",option);
137
138   AliRICHChamber*       iChamber;
139   AliSegmentation*  segmentation;
140  
141   AliRunLoader *pInAL, *pOutAL;//in and out Run loaders
142   AliLoader    *pInRL, *pOutRL;//in and out RICH loaders
143  
144   pOutAL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
145   pOutRL = pOutAL->GetLoader("RICHLoader");
146
147   fTDList = new TObjArray;
148   
149      
150   AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
151   
152   if(!pOutRL->TreeD()) pOutRL->MakeTree("D");  pRICH->MakeBranch("D");
153   
154   
155   fHitMap= new AliHitMap* [kNCH];
156         
157   for (Int_t i =0; i<kNCH; i++) {
158     iChamber= &(pRICH->Chamber(i));
159     segmentation=iChamber->GetSegmentationModel();
160     fHitMap[i] = new AliRICHHitMapA1(segmentation, fTDList);
161   }
162
163 // Loop over files to digitize
164   fSignal = kTRUE;
165   fCounter = 0;
166   for (Int_t inputFile=0;inputFile<fManager->GetNinputs();inputFile++){//files loop
167     pInAL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
168     pInRL = pInAL->GetLoader("RICHLoader");
169     
170 // Connect RICH branches
171
172     if (inputFile > 0 ) fSignal = kFALSE;
173     TBranch *branchHits = 0;
174     TBranch *branchSDigits = 0;
175
176     pInRL->LoadHits("READ");
177     
178     TTree *treeH = pInRL->TreeH();
179     if(pInRL->TreeH()) {
180       branchHits = treeH->GetBranch("RICH");
181       if(branchHits){
182         fHits->Clear();
183         branchHits->SetAddress(&fHits);}
184       else
185         Error("Exec","branch RICH was not found");
186     }
187
188     if(treeH) {
189       branchSDigits = treeH->GetBranch("RICHSpecials");
190       if(branchSDigits) 
191         branchSDigits->SetAddress(&fSDigits);
192       else
193         Error("exec","branch RICHSpecials was not found");
194     }
195     
196     for (fTrack=0; fTrack<treeH->GetEntries(); fTrack++) {//prims loop
197       fHits->Clear();                   fSDigits->Clear();  
198       fHits->Print();
199
200       branchHits->GetEntry(fTrack);     branchSDigits->GetEntry(fTrack);
201       for(Int_t i=0;i<fHits->GetEntriesFast();++i){//hits loop
202         AliRICHhit* pHit = static_cast<AliRICHhit*>(fHits->At(i));
203         fNch = pHit->Chamber()-1;  // chamber number
204         if (fNch >= kNCH) {
205           cerr<<"AliRICHDigitizer: chamber nr. fNch out of range: "<<fNch<<endl;
206           cerr<<"               track: "<<fTrack<<endl;
207           continue;
208         }
209         iChamber = &(pRICH->Chamber(fNch));
210           
211         for (AliRICHSDigit* mPad=FirstPad(pHit,fSDigits);mPad;mPad=NextPad(fSDigits)){//specials loop
212           Int_t iqpad    = mPad->QPad();       // charge per pad
213           fDigits[0]=mPad->PadX();
214           fDigits[1]=mPad->PadY();
215           fDigits[2]=iqpad;
216           fDigits[3]=iqpad;
217           fDigits[4]=mPad->HitNumber();
218           if (Exists(mPad))       // build the list of fired pads and update the info
219             Update(mPad);
220           else
221             CreateNew(mPad);
222         }//specials loop
223       }//hits loop
224     }//prims loop
225   }//files loop
226   if(GetDebug()) Info("Exec","END OF FILE LOOP");
227
228   Int_t tracks[kMAXTRACKSPERRICHDIGIT];
229   Int_t charges[kMAXTRACKSPERRICHDIGIT];
230   Int_t nentries=fTDList->GetEntriesFast();
231     
232   for (Int_t nent=0;nent<nentries;nent++){//transient digits loop
233     AliRICHTransientDigit *transDigit=(AliRICHTransientDigit*)fTDList->At(nent);
234     if (transDigit==0) continue; 
235     Int_t ich=transDigit->GetChamber();
236     Int_t q=transDigit->Signal(); 
237     iChamber=&(pRICH->Chamber(ich));
238     AliRICHResponse * response=iChamber->GetResponseModel();
239     Int_t adcmax= (Int_t) response->MaxAdc();
240       
241       
242     // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
243     //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY()));
244
245 // tmp change
246 //    Int_t pedestal = iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY());
247     Int_t pedestal = 0;
248
249     Float_t treshold = (pedestal + 4*2.2);
250       
251     Float_t meanNoise = gRandom->Gaus(2.2, 0.3);
252     Float_t noise     = gRandom->Gaus(0, meanNoise);
253     q+=(Int_t)(noise + pedestal);
254     //q+=(Int_t)(noise);
255     //          magic number to be parametrised !!! 
256     if ( q <= treshold) 
257     {
258       q = q - pedestal;
259       continue;
260     }
261     q = q - pedestal;
262     if ( q >= adcmax) q=adcmax;
263     fDigits[0]=transDigit->PadX();
264     fDigits[1]=transDigit->PadY();
265     fDigits[2]=q;
266     fDigits[3]=transDigit->Physics();
267     fDigits[4]=transDigit->Hit();
268
269     Int_t nptracks = transDigit->GetNTracks();  
270
271     
272     if(nptracks>kMAXTRACKSPERRICHDIGIT){// this was changed to accomodate the real number of tracks
273       printf("Attention - tracks > 10 %d\n",nptracks);
274       nptracks=kMAXTRACKSPERRICHDIGIT;
275     }
276     if (nptracks > 2) {
277       printf("Attention - tracks > 2  %d \n",nptracks);
278     }
279     for (Int_t tr=0;tr<nptracks;tr++) {
280       tracks[tr]=transDigit->GetTrack(tr);
281       charges[tr]=transDigit->GetCharge(tr);
282     }      //end loop over list of tracks for one pad
283     if (nptracks < kMAXTRACKSPERRICHDIGIT ) {
284       for (Int_t t=nptracks; t<kMAXTRACKSPERRICHDIGIT; t++) {
285       tracks[t]=0;
286       charges[t]=0;
287       }
288     }
289     pRICH->AddDigitOld(ich+1,tracks,charges,fDigits);//OLD
290   }//transient digits loop      
291   pOutRL->TreeD()->Fill();
292
293   //pRICH->ResetDigits();
294   fTDList->Delete();    // or fTDList->Clear(); ???
295   for(Int_t ii=0;ii<kNCH;++ii) {
296     if (fHitMap[ii]) {
297       delete fHitMap[ii];
298       fHitMap[ii]=0;
299     }
300   }
301     
302   TClonesArray *richDigits;
303   for (Int_t k=0;k<kNCH;k++) {
304     richDigits = pRICH->DigitsAddress(k);
305     Int_t ndigit=richDigits->GetEntriesFast();
306     printf ("Chamber %d digits %d \n",k+1,ndigit);
307   }
308   pRICH->ResetDigits(); /// ??? should it be here???
309   
310   pOutRL->WriteDigits("OVERWRITE");
311
312   delete [] fHitMap;
313   delete fTDList;
314   pOutRL->UnloadHits();
315   pOutRL->UnloadDigits();
316 //  if (fHits)    fHits->Delete();
317 //  if (fSDigits) fSDigits->Delete();
318   if(GetDebug())Info("Exec","Stop.");
319 }//Exec()
320 //__________________________________________________________________________________________________
321
322 static Int_t sMaxIterPad=0;    // Static variables for the pad-hit iterator routines
323 static Int_t sCurIterPad=0;
324
325 //__________________________________________________________________________________________________
326 AliRICHSDigit* AliRICHDigitizer::FirstPad(AliRICHhit*  hit,TClonesArray *clusters ) 
327 {// Initialise the pad iterator Return the address of the first sdigit for hit
328     TClonesArray *theClusters = clusters;
329     Int_t nclust = theClusters->GetEntriesFast();
330     if (nclust && hit->PHlast() > 0) {
331         sMaxIterPad=Int_t(hit->PHlast());
332         sCurIterPad=Int_t(hit->PHfirst());
333         return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
334     } else {
335         return 0;
336     }
337     
338 }//FirstPad
339 //__________________________________________________________________________________________________
340 AliRICHSDigit* AliRICHDigitizer::NextPad(TClonesArray *clusters) 
341 {// Iterates over pads
342   
343     sCurIterPad++;
344     if (sCurIterPad <= sMaxIterPad) {
345         return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
346     } else {
347         return 0;
348     }
349 }//NextPad
350 //__________________________________________________________________________________________________