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