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