Protection against missing digit in a digits array in added (P.Hristov)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawDigiProducer.cxx
1 /**************************************************************************
2  * Copyright(c) 2007, 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 //This class produces PHOS digits of one event
19 //using AliPHOSRawDecoder. 
20 //
21 //   For example:
22 //   TClonesArray *digits = new TClonesArray("AliPHOSDigit",100);
23 //   AliRawReader* rawReader = new AliRawReaderDate("2006run2211.raw");
24 //   AliPHOSRawDecoder dc(rawReader);
25 //   while (rawReader->NextEvent()) {
26 //     AliPHOSRawDigiProducer producer;
27 //     producer.MakeDigits(digits,&dc);
28 //   }
29
30 // Author: Boris Polichtchouk
31
32 // --- ROOT system ---
33 #include "TClonesArray.h"
34
35 // --- AliRoot header files ---
36 #include "AliPHOSRawDigiProducer.h"
37 #include "AliPHOSRawDecoder.h"
38 #include "AliPHOSGeometry.h"
39 #include "AliPHOSDigit.h"
40 #include "AliPHOSRecoParam.h"
41 #include "AliPHOSCalibData.h"
42 #include "AliLog.h"
43
44 ClassImp(AliPHOSRawDigiProducer)
45
46 AliPHOSCalibData * AliPHOSRawDigiProducer::fgCalibData  = 0 ; 
47
48 //--------------------------------------------------------------------------------------
49 AliPHOSRawDigiProducer::AliPHOSRawDigiProducer():TObject(),
50   fEmcMinE(0.),fCpvMinE(0.),fEmcCrystals(0),fGeom(0){
51
52 }
53 //--------------------------------------------------------------------------------------
54 AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRecoParam* parEmc, const AliPHOSRecoParam* parCpv):TObject(),
55   fEmcMinE(0.),fCpvMinE(0.),fEmcCrystals(0),fGeom(0){
56
57   if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
58   if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
59
60   fEmcMinE = parEmc->GetMinE();
61   fCpvMinE = parCpv->GetMinE();
62
63   fGeom=AliPHOSGeometry::GetInstance() ;
64   if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP");
65
66   fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ;
67
68   GetCalibrationParameters() ; 
69 }
70 //--------------------------------------------------------------------------------------                       
71 AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp):TObject(),
72   fEmcMinE(0.),fCpvMinE(0.),fEmcCrystals(0),fGeom(0){                                                          
73
74   fEmcMinE = dp.fEmcMinE ;
75   fCpvMinE = dp.fCpvMinE ;
76   fEmcCrystals = dp.fEmcCrystals ;
77   fGeom = dp.fGeom ;
78 }
79 //--------------------------------------------------------------------------------------
80 AliPHOSRawDigiProducer& AliPHOSRawDigiProducer::operator= (const AliPHOSRawDigiProducer &dp){
81   if(&dp == this) return *this;
82
83   fEmcMinE = dp.fEmcMinE ;
84   fCpvMinE = dp.fCpvMinE ;
85   fEmcCrystals = dp.fEmcCrystals ;
86   fGeom = dp.fGeom ;
87   return  *this;
88
89  
90 //--------------------------------------------------------------------------------------
91 void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawDecoder* decoder) 
92 {
93   //Makes the job.
94   //TClonesArray *digits and raw data decoder should be provided by calling function.
95
96   digits->Clear();
97  
98   Int_t    iDigit   = 0 ;
99   Int_t relId[4], absId =0;
100
101   const Double_t baseLine=1. ; //Minimal energy of digit in ADC ch. 
102   const Double_t highLowDiff=2. ; //Maximal difference between High and Low channels in LG adc channels 
103
104   //Temporary array for LowGain digits
105   TClonesArray tmpLG("AliPHOSDigit",10000) ;
106   Int_t ilgDigit=0 ;
107   
108   while (decoder->NextDigit()) {
109
110     Double_t energy=decoder->GetEnergy() ; 
111     if(energy<=baseLine) //in ADC channels
112       continue ;
113
114     Bool_t lowGainFlag = decoder->IsLowGain();
115
116     relId[0] = decoder->GetModule();
117     relId[1] = 0;
118     relId[2] = decoder->GetRow();
119     relId[3] = decoder->GetColumn();
120     fGeom->RelToAbsNumbering(relId, absId);
121
122     Double_t time = decoder->GetTime() ;
123     time = CalibrateT(time,relId,lowGainFlag) ;
124
125     energy = CalibrateE(energy,relId,lowGainFlag) ;
126
127     if(energy <= 0.) 
128        continue;
129
130     if(lowGainFlag){
131       new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
132       ilgDigit++ ; 
133     }
134     else{ 
135       if(decoder->IsOverflow()) //Keep this digit to replace it by Low Gain later.
136                                 //If there is no LogGain it wil be removed by cut on Min E
137         new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time);
138       else
139         new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time);
140       iDigit++;
141     }
142   }
143
144   //Now scan created LG and HG digits and keep only those which appeared in both lists 
145   //replace energy of HighGain digits only if there is overflow
146   //negative energy (overflow)
147   digits->Sort() ;
148   tmpLG.Sort() ;
149   Int_t iLG = 0;
150   Int_t nLG1 = tmpLG.GetEntriesFast()-1 ;
151   for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) { 
152     AliPHOSDigit * digHG = dynamic_cast<AliPHOSDigit*>(digits->At(iDig)) ;
153     if (!digHG) continue;
154     AliPHOSDigit * digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
155     while(iLG<nLG1 && digHG->GetId()> digLG->GetId()){
156       iLG++ ;
157       digLG = dynamic_cast<AliPHOSDigit*>(tmpLG.At(iLG)) ;
158     }
159     if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair
160       if(digHG->GetEnergy()<0.){ //This is overflow in HG
161         digHG->SetTime(digLG->GetTime()) ;
162         digHG->SetEnergy(digLG->GetEnergy()) ;
163       } 
164       else{ //Make approximate comparison of HG and LG energies
165         Double_t de = (digHG->GetEnergy()-digLG->GetEnergy())/16./0.005 ; //aproximate difference in LG ADC ch. 
166         if(TMath::Abs(de)>highLowDiff){ //too strong difference, remove digit
167           digits->RemoveAt(iDig) ;
168         }
169       }
170     }
171     else{ //no pair - remove
172       // temporary fix for dead LG channels
173       Int_t absId=digHG->GetId() ;
174       Int_t relId[4] ;
175       fGeom->AbsToRelNumbering(absId,relId) ;
176       if(relId[2]%2==1 && relId[3]%16==4) 
177         continue ;
178       digits->RemoveAt(iDig) ;                                                                                                            
179     }
180   }
181
182   CleanDigits(digits) ;
183
184 }
185 //____________________________________________________________________________
186 Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain)
187 {
188   // Convert EMC measured amplitude into real energy.
189   // Calibration parameters are taken from calibration data base for raw data,
190   // or from digitizer parameters for simulated data.                        
191   if(fgCalibData){ 
192     Int_t   module = relId[0];  
193     Int_t   column = relId[3];
194     Int_t   row    = relId[2];
195     if(relId[1]==0) { // this is EMC 
196       if(isLowGain){
197         amp*= fgCalibData->GetHighLowRatioEmc(module,column,row);
198       }
199       amp *= fgCalibData->GetADCchannelEmc(module,column,row);                                                          
200       return amp ;         
201     }         
202   }          
203   return 0;        
204 }
205 //____________________________________________________________________________
206 Double_t AliPHOSRawDigiProducer::CalibrateT(Double_t time, Int_t * relId, Bool_t /* isLowGain */)
207 {
208   //Calibrate time
209   if(fgCalibData){
210     Int_t   module = relId[0];
211     Int_t   column = relId[3];
212     Int_t   row    = relId[2];
213     if(relId[1]==0) { // this is EMC
214       time += fgCalibData->GetTimeShiftEmc(module,column,row);                   
215       return time ;             
216     }
217   }
218  
219   return -999.;
220 }
221 //____________________________________________________________________________
222 void AliPHOSRawDigiProducer::CleanDigits(TClonesArray * digits)
223 {
224   // remove digits with amplitudes below threshold.
225   // remove digits in bad channels
226   // sort digits with icreasing AbsId
227   
228   //remove digits in bad map and below threshold
229   Bool_t isBadMap = 0 ;
230   if(fgCalibData->GetNumOfEmcBadChannels()){
231     isBadMap=1 ;
232   }
233   
234   for(Int_t i=0; i<digits->GetEntriesFast(); i++){
235     AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
236     if(!digit)
237       continue  ;
238     if ( (IsInEMC(digit) && digit->GetEnergy() < fEmcMinE) ||
239          (IsInCPV(digit) && digit->GetEnergy() < fCpvMinE) ){
240       digits->RemoveAt(i) ;
241       continue ;
242     }
243     if(isBadMap){ //check bad map now
244       Int_t relid[4] ;
245       fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
246       if(fgCalibData->IsBadChannelEmc(relid[0],relid[2],relid[3])){
247         digits->RemoveAt(i) ;
248       }
249     }
250   }
251
252   //Compress, sort and set indexes
253   digits->Compress() ;
254 //  digits->Sort(); already sorted earlier
255   for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) { 
256     AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ; 
257     digit->SetIndexInList(i) ;     
258   }
259 }
260 //____________________________________________________________________________
261 Bool_t AliPHOSRawDigiProducer::IsInEMC(AliPHOSDigit * digit) const
262 {
263   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
264   return digit->GetId() <= fEmcCrystals ;
265
266 }
267
268 //____________________________________________________________________________
269 Bool_t AliPHOSRawDigiProducer::IsInCPV(AliPHOSDigit * digit) const
270 {
271   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
272   return digit->GetId() > fEmcCrystals ;
273 }
274 //____________________________________________________________________________
275 void AliPHOSRawDigiProducer::GetCalibrationParameters() 
276 {
277   // Set calibration parameters:
278   // if calibration database exists, they are read from database,
279   // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
280   //
281   // It is a user responsilibity to open CDB before reconstruction, for example: 
282   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
283
284   if (!fgCalibData){
285     fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
286   }
287   if (fgCalibData->GetCalibDataEmc() == 0)
288     AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
289   if (fgCalibData->GetCalibDataCpv() == 0)
290     AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
291 }