add ADDigits branch to ADDigits tree in every event
[u/mrichter/AliRoot.git] / AD / ADrec / AliADReconstructor.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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: AliADReconstructor.cxx 20956 2007-09-26 14:22:18Z mrodrigu $ */
17 //////////////////////////////////////////////////////////////////////////////
18 //                                                                          //
19 //  Class for AD reconstruction                                         //
20 //////////////////////////////////////////////////////////////////////////////
21 #include <TParameter.h>
22
23 #include "AliRawReader.h"
24 #include "AliGRPObject.h"
25 #include "AliCDBManager.h"
26 #include "AliCDBStorage.h"
27 #include "AliCDBEntry.h"
28 #include "AliESDEvent.h"
29 #include "AliRunInfo.h"
30 #include "AliCTPTimeParams.h"
31 #include "AliLHCClockPhase.h"
32
33 #include "AliADReconstructor.h"
34 #include "AliADdigit.h"
35 #include "AliESDAD.h"
36 #include "AliESDADfriend.h"
37 #include "AliADConst.h"
38 #include "AliADCalibData.h"
39 #include "AliADRawStream.h"
40
41 ClassImp(AliADReconstructor)
42 //_____________________________________________________________________________
43 AliADReconstructor:: AliADReconstructor():
44   AliReconstructor(),
45   fESDAD(0x0),
46   fESD(0x0),
47   fESDADfriend(0x0),
48   fCalibData(NULL),
49   fDigitsArray(0),
50   fCollisionMode(0),
51   fBeamEnergy(0.)
52
53 {
54   fCalibData = GetCalibData();
55   // Default constructor  
56   // Get calibration data
57
58 }
59
60 //_____________________________________________________________________________
61 AliADReconstructor& AliADReconstructor::operator = 
62   (const AliADReconstructor& /*reconstructor*/)
63 {
64 // assignment operator
65
66   Fatal("operator =", "assignment operator not implemented");
67   return *this;
68 }
69
70 //_____________________________________________________________________________
71 AliADReconstructor::~AliADReconstructor()
72 {
73 // destructor
74   delete fESDAD;
75   delete fESDADfriend;
76   delete fDigitsArray;
77 }
78
79 //_____________________________________________________________________________
80 void AliADReconstructor::Init()
81 {
82 // initializer
83     fESDAD  = new AliESDAD;
84     fESDADfriend  = new AliESDADfriend;
85 }
86
87 //_____________________________________________________________________________
88 void AliADReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
89 {
90 // converts RAW to digits 
91
92   if (!digitsTree) {
93     AliError("No digits tree!");
94     return;
95   }
96
97   if (!fDigitsArray) fDigitsArray = new TClonesArray("AliADdigit", 16);
98   digitsTree->Branch("ADDigit", &fDigitsArray);
99
100   rawReader->Reset();
101   AliADRawStream rawStream(rawReader);
102   if (rawStream.Next()) { 
103
104     for(Int_t iChannel=0; iChannel < 16; ++iChannel) {
105       Int_t offlineCh = kOfflineChannel[iChannel];
106       // ADC charge samples
107       Short_t chargeADC[kNClocks];
108       for(Int_t iClock=0; iClock < kNClocks; ++iClock) {
109         chargeADC[iClock] = rawStream.GetPedestal(iChannel,iClock);
110       }
111       // Integrator flag
112       Bool_t integrator = rawStream.GetIntegratorFlag(iChannel,kNClocks/2);
113       Bool_t BBflag = rawStream.GetBBFlag(iChannel,kNClocks/2); 
114       Bool_t BGflag = rawStream.GetBGFlag(iChannel,kNClocks/2);
115    
116       // HPTDC data (leading time and width)
117       Int_t board = AliADCalibData::GetBoardNumber(offlineCh);
118       Float_t time = rawStream.GetTime(iChannel)*fCalibData->GetTimeResolution(board);
119       Float_t width = rawStream.GetWidth(iChannel)*fCalibData->GetWidthResolution(board);
120       // Add a digit
121       if(!fCalibData->IsChannelDead(iChannel)){
122           new ((*fDigitsArray)[fDigitsArray->GetEntriesFast()]) AliADdigit(offlineCh, time, width,integrator, chargeADC, BBflag, BGflag);
123       }
124       
125       fESDADfriend->SetTriggerInputs(rawStream.GetTriggerInputs());
126       fESDADfriend->SetTriggerInputsMask(rawStream.GetTriggerInputsMask());
127     
128       fESDADfriend->SetBBScalers(offlineCh,rawStream.GetBBScalers(iChannel));
129       fESDADfriend->SetBGScalers(offlineCh,rawStream.GetBGScalers(iChannel));
130       for (Int_t iEv = 0; iEv < kNClocks; iEv++) {
131           fESDADfriend->SetBBFlag(offlineCh,iEv,rawStream.GetBBFlag(iChannel,iEv));
132           fESDADfriend->SetBGFlag(offlineCh,iEv,rawStream.GetBGFlag(iChannel,iEv));
133       }
134     }
135     for(Int_t iScaler = 0; iScaler < AliESDADfriend::kNScalers; iScaler++)
136       fESDADfriend->SetTriggerScalers(iScaler,rawStream.GetTriggerScalers(iScaler));
137
138     digitsTree->Fill();
139   }
140
141   fDigitsArray->Clear();
142
143 }
144
145 //_____________________________________________________________________________
146 void AliADReconstructor::FillESD(TTree* digitsTree, TTree* /*clustersTree*/,AliESDEvent* esd) const
147 {
148
149   printf("Running AD Reconstruction \n");
150
151   if (!digitsTree) {
152       AliError("No digits tree!");
153       return;
154   }
155
156   TBranch* digitBranch = digitsTree->GetBranch("ADDigit");
157   digitBranch->SetAddress(&fDigitsArray);
158
159   Float_t   mult[16];  
160   Float_t    adc[16]; 
161   Float_t   time[16]; 
162   Float_t  width[16];
163   Bool_t aBBflag[16];
164   Bool_t aBGflag[16];
165    
166   for (Int_t i=0; i<16; i++){
167        adc[i]    = 0.0;
168        mult[i]   = 0.0;
169        time[i]   = kInvalidTime;
170        width[i]  = 0.0;
171        aBBflag[i] = kFALSE;
172        aBGflag[i] = kFALSE;
173   }
174
175   Int_t nEntries = (Int_t)digitsTree->GetEntries();
176   for (Int_t e=0; e<nEntries; e++) {
177     digitsTree->GetEvent(e);
178
179     Int_t nDigits = fDigitsArray->GetEntriesFast();
180     
181     for (Int_t d=0; d<nDigits; d++) {    
182         AliADdigit* digit = (AliADdigit*) fDigitsArray->At(d);      
183         Int_t  pmNumber = digit->PMNumber();
184
185         // Pedestal retrieval and suppression
186         Bool_t integrator = digit->Integrator();
187         Float_t maxadc = 0;
188         Int_t imax = -1;
189         Float_t adcPedSub[kNClocks];
190         for(Int_t iClock=0; iClock < kNClocks; ++iClock) {
191           Short_t charge = digit->ChargeADC(iClock);
192           Bool_t iIntegrator = (iClock%2 == 0) ? integrator : !integrator;
193           Int_t k = pmNumber + 16*iIntegrator;
194           adcPedSub[iClock] = (Float_t)charge - fCalibData->GetPedestal(k);
195           //if(adcPedSub[iClock] <= GetRecoParam()->GetNSigmaPed()*fCalibData->GetSigma(k)) {
196           if(adcPedSub[iClock] <= 2*fCalibData->GetSigma(k)) {
197             adcPedSub[iClock] = 0;
198             continue;
199           }
200           //if(iClock < GetRecoParam()->GetStartClock() || iClock > GetRecoParam()->GetEndClock()) continue;
201           if(iClock < 8 || iClock > 12) continue;
202           if(adcPedSub[iClock] > maxadc) {
203             maxadc = adcPedSub[iClock];
204             imax   = iClock;
205           }
206         }
207
208         if (imax != -1) {
209           //Int_t start = imax - GetRecoParam()->GetNPreClocks();
210           Int_t start = imax - 2;
211           if (start < 0) start = 0;
212           //Int_t end = imax + GetRecoParam()->GetNPostClocks();
213           Int_t end = imax + 1;
214           if (end > 20) end = 20;
215           for(Int_t iClock = start; iClock <= end; iClock++) {
216             adc[pmNumber] += adcPedSub[iClock];
217           }
218         }
219
220         // HPTDC leading time and width
221         // Correction for slewing and various time delays
222         //time[pmNumber]  =  CorrectLeadingTime(pmNumber,digit->Time(),adc[pmNumber]);
223         time[pmNumber]  =  digit->Time();
224         width[pmNumber] =  digit->Width();
225         aBBflag[pmNumber] = digit->GetBBflag();
226         aBGflag[pmNumber] = digit->GetBGflag();
227
228         if (adc[pmNumber] > 0) {
229           AliDebug(1,Form("PM = %d ADC = %.2f (%.2f) TDC %.2f (%.2f)   Int %d (%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d)    %.2f %.2f   %.2f %.2f    %d %d",pmNumber, adc[pmNumber],
230                        digit->ChargeADC(11)+digit->ChargeADC(10)+digit->ChargeADC(9)+digit->ChargeADC(8)+
231                        digit->ChargeADC(7)+digit->ChargeADC(6)+digit->ChargeADC(5)+digit->ChargeADC(4)-
232                        4.*fCalibData->GetPedestal(pmNumber)-4.*fCalibData->GetPedestal(pmNumber+16),
233                           digit->Time(),time[pmNumber],
234                           integrator,
235                           digit->ChargeADC(0),digit->ChargeADC(1),digit->ChargeADC(2),digit->ChargeADC(3),digit->ChargeADC(4),digit->ChargeADC(5),digit->ChargeADC(6),digit->ChargeADC(7),
236                           digit->ChargeADC(8),digit->ChargeADC(9),digit->ChargeADC(10),
237                           digit->ChargeADC(11),digit->ChargeADC(12),
238                           digit->ChargeADC(13),digit->ChargeADC(14),digit->ChargeADC(15),digit->ChargeADC(16),digit->ChargeADC(17),digit->ChargeADC(18),digit->ChargeADC(19),digit->ChargeADC(20),
239                           fCalibData->GetPedestal(pmNumber),fCalibData->GetSigma(pmNumber),
240                           fCalibData->GetPedestal(pmNumber+16),fCalibData->GetSigma(pmNumber+16),
241                           aBBflag[pmNumber],aBGflag[pmNumber]));
242             };
243
244         // Fill ESD friend object
245         for (Int_t iEv = 0; iEv < kNClocks; iEv++) {
246           fESDADfriend->SetPedestal(pmNumber,iEv,(Float_t)digit->ChargeADC(iEv));
247           fESDADfriend->SetIntegratorFlag(pmNumber,iEv,(iEv%2 == 0) ? integrator : !integrator);
248         }
249         fESDADfriend->SetTime(pmNumber,digit->Time());
250         fESDADfriend->SetWidth(pmNumber,digit->Width());
251
252     } // end of loop over digits
253   } // end of loop over events in digits tree
254          
255   //fESDAD->SetBit(AliESDAD::kCorrectedLeadingTime,kTRUE);
256   //fESDAD->SetMultiplicity(mult);
257   fESDAD->SetADC(adc);
258   fESDAD->SetTime(time);
259   fESDAD->SetWidth(width);
260   fESDAD->SetBit(AliESDAD::kOnlineBitsFilled,kTRUE);
261   fESDAD->SetBBFlag(aBBflag);
262   fESDAD->SetBGFlag(aBGflag);
263   //fESDAD->SetBit(AliESDAD::kCorrectedForSaturation,kTRUE);
264
265   /*/ now fill the AD decision and channel flags
266   {
267     AliADDecision offlineDecision;
268     offlineDecision.SetRecoParam(GetRecoParam());
269     offlineDecision.FillDecisions(fESDAD, fCalibData, fTimeSlewing);
270   }/*/
271
272   if (esd) { 
273      AliDebug(1, Form("Writing AD data to ESD tree"));
274      esd->SetADData(fESDAD);
275  
276      AliESDfriend *fr = (AliESDfriend*)esd->FindListObject("AliESDfriend");
277      if (fr) {
278         AliDebug(1, Form("Writing AD friend data to ESD tree"));
279         fr->SetADfriend(fESDADfriend);
280     }
281   }
282
283   fDigitsArray->Clear();
284
285 }
286
287 //_____________________________________________________________________________
288 AliADCalibData* AliADReconstructor::GetCalibData() const
289 {
290   AliCDBManager *man = AliCDBManager::Instance();
291
292   AliCDBEntry *entry=0;
293
294   entry = man->Get("AD/Calib/Data");
295   if(!entry){
296     AliWarning("Load of calibration data from default storage failed!");
297     AliWarning("Calibration data will be loaded from local storage ($ALICE_ROOT)");
298         
299     man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
300     man->SetRun(1);
301     entry = man->Get("AD/Calib/Data");
302   }
303   // Retrieval of data in directory AD/Calib/Data:
304
305   AliADCalibData *calibdata = 0;
306
307   if (entry) calibdata = (AliADCalibData*) entry->GetObject();
308   if (!calibdata)  AliFatal("No calibration data from calibration database !");
309
310   return calibdata;
311 }
312
313 //_____________________________________________________________________________
314 AliCDBStorage* AliADReconstructor::SetStorage(const char *uri) 
315 {
316 // Sets the storage  
317
318   Bool_t deleteManager = kFALSE;
319   
320   AliCDBManager *manager = AliCDBManager::Instance();
321   AliCDBStorage *defstorage = manager->GetDefaultStorage();
322   
323   if(!defstorage || !(defstorage->Contains("AD"))){ 
324      AliWarning("No default storage set or default storage doesn't contain AD!");
325      manager->SetDefaultStorage(uri);
326      deleteManager = kTRUE;
327   }
328  
329   AliCDBStorage *storage = manager->GetDefaultStorage();
330
331   if(deleteManager){
332      AliCDBManager::Instance()->UnsetDefaultStorage();
333      defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
334   }
335
336   return storage; 
337 }
338
339 //____________________________________________________________________________
340 void AliADReconstructor::GetCollisionMode()
341 {
342   // Retrieval of collision mode 
343
344   TString beamType = GetRunInfo()->GetBeamType();
345   if(beamType==AliGRPObject::GetInvalidString()){
346      AliError("AD cannot retrieve beam type");
347      return;
348   }
349
350   if( (beamType.CompareTo("P-P") ==0)  || (beamType.CompareTo("p-p") ==0) ){
351     fCollisionMode=0;
352   }
353   else if( (beamType.CompareTo("Pb-Pb") ==0)  || (beamType.CompareTo("A-A") ==0) ){
354     fCollisionMode=1;
355   }
356     
357   fBeamEnergy = GetRunInfo()->GetBeamEnergy();
358   if(fBeamEnergy==AliGRPObject::GetInvalidFloat()) {
359      AliError("Missing value for the beam energy ! Using 0");
360      fBeamEnergy = 0.;
361   }
362   
363   AliDebug(1,Form("\n ++++++ Beam type and collision mode retrieved as %s %d @ %1.3f GeV ++++++\n\n",beamType.Data(), fCollisionMode, fBeamEnergy));
364
365 }