removed possibility to choose channels for stability
[u/mrichter/AliRoot.git] / T0 / AliT0Reconstructor.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$ */
17 /*********************************************************************
18  *  T0 reconstruction and filling ESD
19  *  - reconstruct mean time (interation time) 
20  *  - vertex position
21  *  -  multiplicity
22  ********************************************************************/
23
24 #include <AliESDEvent.h>
25 #include "AliLog.h"
26 #include "AliT0RecPoint.h"
27 #include "AliRawReader.h"
28 #include "AliT0RawReader.h"
29 #include "AliT0digit.h"
30 #include "AliT0Reconstructor.h"
31 #include "AliT0Parameters.h"
32 #include "AliT0Calibrator.h"
33 #include "AliESDfriend.h"
34 #include "AliESDTZEROfriend.h"
35 #include "AliLog.h"
36 #include "AliCDBEntry.h" 
37 #include "AliCDBManager.h"
38 #include "AliCTPTimeParams.h"
39 #include "AliLHCClockPhase.h"
40 #include "AliT0CalibSeasonTimeShift.h"
41 #include "AliESDRun.h"
42
43 #include <TArrayI.h>
44 #include <TGraph.h>
45 #include <TMath.h>
46 #include <Riostream.h>
47
48 ClassImp(AliT0Reconstructor)
49
50   AliT0Reconstructor:: AliT0Reconstructor(): AliReconstructor(),
51                                              fdZonA(0),
52                                              fdZonC(0),
53                                              fZposition(0),
54                                              fParam(NULL),
55                                              fAmpLEDrec(),
56                                              fQTC(0),
57                                              fAmpLED(0),
58                                              fCalib(),
59                                              fLatencyHPTDC(9000),
60                                              fLatencyL1(0),
61                                              fLatencyL1A(0),
62                                              fLatencyL1C(0),
63                                              fGRPdelays(0),
64                                              fTimeMeanShift(0x0),
65                                              fTimeSigmaShift(0x0),
66                                              fESDTZEROfriend(NULL)
67
68 {
69   for (Int_t i=0; i<24; i++)  fTime0vertex[i] =0;
70
71   //constructor
72   AliCDBEntry *entry = AliCDBManager::Instance()->Get("GRP/CTP/CTPtiming");
73   if (!entry) AliFatal("CTP timing parameters are not found in OCDB !");
74   AliCTPTimeParams *ctpParams = (AliCTPTimeParams*)entry->GetObject();
75   Float_t l1Delay = (Float_t)ctpParams->GetDelayL1L0()*25.0;
76
77   AliCDBEntry *entry1 = AliCDBManager::Instance()->Get("GRP/CTP/TimeAlign");
78   if (!entry1) AliFatal("CTP time-alignment is not found in OCDB !");
79   AliCTPTimeParams *ctpTimeAlign = (AliCTPTimeParams*)entry1->GetObject();
80   l1Delay += ((Float_t)ctpTimeAlign->GetDelayL1L0()*25.0);
81  
82   AliCDBEntry *entry4 = AliCDBManager::Instance()->Get("GRP/Calib/LHCClockPhase");
83   if (!entry4) AliFatal("LHC clock-phase shift is not found in OCDB !");
84   AliLHCClockPhase *phase = (AliLHCClockPhase*)entry4->GetObject();
85
86   fGRPdelays = l1Delay - phase->GetMeanPhase();
87
88   AliCDBEntry *entry5 = AliCDBManager::Instance()->Get("T0/Calib/TimeAdjust");
89   if (entry5) {
90     AliT0CalibSeasonTimeShift *timeshift = (AliT0CalibSeasonTimeShift*)entry5->GetObject();
91     fTimeMeanShift = timeshift->GetT0Means();
92     fTimeSigmaShift  = timeshift->GetT0Sigmas();
93    }
94   else
95     AliWarning("Time Adjust is not found in OCDB !");
96  
97   fParam = AliT0Parameters::Instance();
98   fParam->Init();
99  
100   for (Int_t i=0; i<24; i++){
101         TGraph* gr = fParam ->GetAmpLEDRec(i);
102         if (gr) fAmpLEDrec.AddAtAndExpand(gr,i) ; 
103           TGraph* gr1 = fParam ->GetAmpLED(i);
104           if (gr1) fAmpLED.AddAtAndExpand(gr1,i) ; 
105           TGraph* gr2 = fParam ->GetQTC(i);
106           if (gr2) fQTC.AddAtAndExpand(gr2,i) ;         
107           fTime0vertex[i] = fParam->GetCFD(i);
108           printf("OCDB mean CFD time %i %f \n",i, fTime0vertex[i]);
109  }
110   fLatencyL1 = fParam->GetLatencyL1();
111   fLatencyL1A = fParam->GetLatencyL1A(); 
112   fLatencyL1C = fParam->GetLatencyL1C();
113   fLatencyHPTDC = fParam->GetLatencyHPTDC();
114   AliDebug(2,Form(" LatencyL1 %f latencyL1A %f latencyL1C %f latencyHPTDC %f \n",fLatencyL1, fLatencyL1A, fLatencyL1C, fLatencyHPTDC));
115  
116   for (Int_t i=0; i<24; i++) {
117     if( fTime0vertex[i] < 500 || fTime0vertex[i] > 50000) fTime0vertex[i] =( 1000.*fLatencyHPTDC - 1000.*fLatencyL1 + 1000.*fGRPdelays)/24.4;
118  
119   }
120   
121   // fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1"));
122   //fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));
123   //here real Z position
124   fdZonC = TMath::Abs(fParam->GetZPosition("T0/C/PMT1"));
125   fdZonA = TMath::Abs(fParam->GetZPosition("T0/A/PMT15"));
126
127   fCalib = new AliT0Calibrator();
128   fESDTZEROfriend = new AliESDTZEROfriend();
129
130 }
131
132 //_____________________________________________________________________________
133 void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const
134 {
135   // T0 digits reconstruction
136   Int_t refAmp = Int_t (GetRecoParam()->GetRefAmp());
137   
138   TArrayI * timeCFD = new TArrayI(24); 
139   TArrayI * timeLED = new TArrayI(24); 
140   TArrayI * chargeQT0 = new TArrayI(24); 
141   TArrayI * chargeQT1 = new TArrayI(24); 
142
143  
144   Float_t channelWidth = fParam->GetChannelWidth() ;  
145   Float_t meanVertex = fParam->GetMeanVertex();
146   Float_t c = 0.0299792; // cm/ps
147   Double32_t vertex = 9999999;
148   Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;
149
150   
151   AliDebug(1,Form("Start DIGITS reconstruction "));
152   
153   Float_t lowAmpThreshold =  GetRecoParam()->GetLow(200);  
154   Float_t highAmpThreshold =  GetRecoParam()->GetHigh(200);  
155   Int_t badpmt = GetRecoParam()->GetRefPoint();
156
157   TBranch *brDigits=digitsTree->GetBranch("T0");
158   AliT0digit *fDigits = new AliT0digit() ;
159   if (brDigits) {
160     brDigits->SetAddress(&fDigits);
161   }else{
162     AliError(Form("EXEC Branch T0 digits not found"));
163      return;
164   }
165   
166   digitsTree->GetEvent(0);
167   digitsTree->GetEntry(0);
168   brDigits->GetEntry(0);
169   fDigits->GetTimeCFD(*timeCFD);
170   fDigits->GetTimeLED(*timeLED);
171   fDigits->GetQT0(*chargeQT0);
172   fDigits->GetQT1(*chargeQT1);
173   Int_t onlineMean =  fDigits->MeanTime();
174
175   Bool_t tr[5];
176   for (Int_t i=0; i<5; i++) tr[i]=false; 
177   
178   Double32_t besttimeA=999999;
179   Double32_t besttimeC=999999;
180   Int_t pmtBestA=99999;
181   Int_t pmtBestC=99999;
182   
183   AliT0RecPoint* frecpoints= new AliT0RecPoint ();
184   clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints);
185   
186   Float_t time[24], adc[24], adcmip[24];
187   for (Int_t ipmt=0; ipmt<24; ipmt++) {
188     if(timeCFD->At(ipmt)>0 && ipmt != badpmt) {
189      if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)  
190         adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
191       else
192         adc[ipmt] = 0;
193       
194      time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, Int_t(adc[ipmt]),  timeCFD->At(ipmt)) ;
195              
196       Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
197       //    time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl),  timeCFD->At(ipmt) ) ;
198       AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
199                        ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
200
201       Double_t ampMip =((TGraph*)fAmpLED.At(ipmt))->Eval(sl);
202       Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
203       AliDebug(5,Form("  Amlitude in MIPS LED %f ,  QTC %f in channels %f\n ",ampMip,qtMip, adc[ipmt]));
204       
205       frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
206       frecpoints->SetAmpLED(ipmt, Float_t( ampMip)); 
207       frecpoints->SetAmp(ipmt, Float_t(qtMip));
208       adcmip[ipmt]=qtMip;
209       
210     }
211     else {
212       time[ipmt] = 0;
213       adc[ipmt] = 0;
214       adcmip[ipmt] = 0;
215
216     }
217   }
218   
219   for (Int_t ipmt=0; ipmt<12; ipmt++){
220     if(time[ipmt] > 1 && ipmt != badpmt && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold) {
221       if(time[ipmt]<besttimeC){
222         besttimeC=time[ipmt]; //timeC
223         pmtBestC=ipmt;
224       }
225     }
226   }
227   for ( Int_t ipmt=12; ipmt<24; ipmt++){
228     if(time[ipmt] > 1 &&  ipmt != badpmt && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold) {
229       if(time[ipmt]<besttimeA) {
230         besttimeA=time[ipmt]; //timeA
231         pmtBestA=ipmt;}
232     }
233   }
234   if(besttimeA < 999999) {
235     frecpoints->SetTimeBestA(Int_t(besttimeA *channelWidth - fdZonA/c));
236     tr[1]=true;
237   }
238   if( besttimeC < 999999 ) {
239     frecpoints->SetTimeBestC(Int_t(besttimeC *channelWidth - fdZonA/c));
240     tr[2]=true;
241   }
242   AliDebug(5,Form(" besttimeA %f ch,  besttimeC %f ch",besttimeA, besttimeC));
243   if(besttimeA <999999 && besttimeC < 999999 ){
244     //    timeDiff = (besttimeC - besttimeA)*channelWidth;
245     timeDiff = (besttimeA - besttimeC)*channelWidth;
246     meanTime = (besttimeA + besttimeC)/2;// * channelWidth); 
247     timeclock = meanTime *channelWidth -fdZonA/c ;
248     vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
249     tr[0]=true; 
250   }
251   frecpoints->SetVertex(vertex);
252   frecpoints->SetMeanTime(meanTime);
253   frecpoints->SetT0clock(timeclock);
254   frecpoints->SetT0Trig(tr);
255
256   AliDebug(5,Form("T0 triggers %d %d %d %d %d",tr[0],tr[1],tr[2],tr[3],tr[4]));
257
258   //online mean
259   frecpoints->SetOnlineMean(Int_t(onlineMean));
260   AliDebug(10,Form("  timeDiff %f #channel,  meanTime %f #channel, vertex %f cm online mean %i timeclock %f ps",timeDiff, meanTime,vertex, Int_t(onlineMean), timeclock));
261   
262   
263
264    
265   
266   clustersTree->Fill();
267
268   delete timeCFD;
269   delete timeLED;
270   delete chargeQT0; 
271   delete chargeQT1; 
272 }
273
274
275 //_______________________________________________________________________
276
277 void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
278 {
279   // T0 raw ->
280   //
281   // reference amplitude and time ref. point from reco param
282
283   // Float_t refAmp = GetRecoParam()->GetRefAmp();
284
285   //  Int_t refPoint = 0;
286   Int_t badpmt[24];
287   //Bad channel
288   for (Int_t i=0; i<24; i++) {
289     badpmt[i] = GetRecoParam() -> GetBadChannels(i);
290   }
291   Int_t low[500], high[500];
292
293   Int_t allData[110][5];
294   
295   Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
296   Double32_t timeDiff, meanTime, timeclock;
297   timeDiff =  meanTime = timeclock = 9999999;
298   Float_t c = 29.9792458; // cm/ns
299   Double32_t vertex = 9999999;
300   Int_t onlineMean=0;
301   // Float_t meanVertex = fParam->GetMeanVertex();
302   Float_t meanVertex = 0;
303   for (Int_t i0=0; i0<24; i0++) {
304     low[i0] = Int_t(fTime0vertex[i0]) - 200;
305     high[i0] = Int_t(fTime0vertex[i0]) + 200;
306   }
307   
308   for (Int_t i0=0; i0<110; i0++)
309     {
310       for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0; 
311       //    low[i0] = Int_t (GetRecoParam()->GetLow(i0));       
312       // high[i0] = Int_t (GetRecoParam()->GetHigh(i0));
313       }
314   
315   Float_t lowAmpThreshold =  GetRecoParam()->GetAmpLowThreshold();  
316   Float_t highAmpThreshold =  GetRecoParam()->GetAmpHighThreshold(); 
317   
318   Double32_t besttimeA=9999999;
319   Double32_t besttimeC=9999999;
320   Int_t pmtBestA=99999;
321   Int_t pmtBestC=99999;
322    
323   AliT0RecPoint* frecpoints= new AliT0RecPoint ();
324   
325   recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints);
326    
327   AliDebug(10," before read data ");
328   AliT0RawReader myrawreader(rawReader);
329
330   UInt_t type =rawReader->GetType();
331
332   if (!myrawreader.Next())
333     AliDebug(1,Form(" no raw data found!!"));
334   else
335     {  
336    for (Int_t i=0; i<24; i++)
337     {
338       timeCFD[i]=0; timeLED[i]=0; chargeQT0[i]=0; chargeQT1[i]=0;
339     }
340      Int_t fBCID=Int_t (rawReader->GetBCID());
341       Int_t trmbunch= myrawreader.GetTRMBunchID();
342       AliDebug(10,Form(" CDH BC ID %i, TRM BC ID %i \n", fBCID, trmbunch ));
343  
344       if(type == 7  ) {  //only physics 
345         for (Int_t i=0; i<107; i++) {
346         for (Int_t iHit=0; iHit<5; iHit++) 
347           {
348             allData[i][iHit] = myrawreader.GetData(i,iHit);
349           }
350         }
351         Int_t ref=0;
352
353         //      if (refPoint>0) 
354         //  ref = allData[refPoint][0]-5000;
355
356         
357         Float_t channelWidth = fParam->GetChannelWidth() ;  
358         
359         //       Int_t meanT0 = fParam->GetMeanT0();
360         
361         for (Int_t in=0; in<12; in++)  
362           {
363             for (Int_t iHit=0; iHit<5; iHit++) 
364               {
365                 if(allData[in+1][iHit] > low[in] && 
366                    allData[in+1][iHit] < high[in])
367                   {
368                     timeCFD[in] = allData[in+1][iHit] ; 
369                     break;
370                   }
371               }
372             for (Int_t iHit=0; iHit<5; iHit++) 
373               {
374                 if(allData[in+1+56][iHit] > low[in] && 
375                    allData[in+1+56][iHit] < high[in])
376                   {
377                     timeCFD[in+12] = allData[in+56+1][iHit] ;
378                     break;
379                   }
380               }
381             timeLED[in+12] = allData[in+68+1][0] ;
382             timeLED[in] = allData[in+12+1][0] ;
383             AliDebug(5, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
384                               in, timeCFD[in],timeCFD[in+12],timeLED[in], 
385                              timeLED[in+12]));   
386             
387           }
388         
389         
390         for (Int_t in=0; in<12;  in++)
391           {
392             chargeQT0[in]=allData[2*in+25][0];
393             chargeQT1[in]=allData[2*in+26][0];
394             AliDebug(25, Form(" readed Raw %i %i %i",
395                               in, chargeQT0[in],chargeQT1[in]));
396           }     
397         for (Int_t in=12; in<24;  in++)
398           {
399             chargeQT0[in]=allData[2*in+57][0];
400             chargeQT1[in]=allData[2*in+58][0];
401             AliDebug(25, Form(" readed Raw %i %i %i",
402                               in, chargeQT0[in],chargeQT1[in]));
403             
404           }
405         
406         for (Int_t iHit=0; iHit<5; iHit++) 
407           {
408               if(allData[49][iHit] > low[49] && 
409                  allData[49][iHit] < high[49]){
410                 onlineMean = allData[49][iHit];
411                 break;
412               }
413           }
414         Double32_t time[24], adc[24], adcmip[24], noncalibtime[24];
415         for (Int_t ipmt=0; ipmt<24; ipmt++) {
416           if(timeCFD[ipmt] >  0 /* && badpmt[ipmt]==0*/ ){
417            //for simulated data
418              //for physics  data
419            if(( chargeQT0[ipmt] - chargeQT1[ipmt])>0)  {
420              adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
421            }
422            else
423              adc[ipmt] = 0;
424            //      time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
425            
426            time[ipmt] = fCalib-> WalkCorrection(Int_t (fTime0vertex[ipmt]), ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
427            Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
428            // time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
429            AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
430                             ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
431            Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
432            Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
433            AliDebug(10,Form("  Amlitude in MIPS LED %f ; QTC %f;  in channels %f\n ",ampMip,qtMip, adc[ipmt]));
434            //bad peak removing
435              frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
436              // frecpoints->SetTime(ipmt,Double32_t(timeCFD[ipmt]));
437              frecpoints->SetAmp(ipmt, Double32_t( qtMip)); 
438              adcmip[ipmt]=qtMip;
439              frecpoints->SetAmpLED(ipmt, Double32_t(ampMip));        
440              noncalibtime[ipmt]= Double32_t (timeCFD[ipmt]);
441          }
442          else {
443            time[ipmt] = 0;
444            adc[ipmt] = 0;
445            adcmip[ipmt] = 0;
446            noncalibtime[ipmt] = 0;
447          }
448        }
449        fESDTZEROfriend->SetT0timeCorr(noncalibtime) ;     
450        for (Int_t ipmt=0; ipmt<12; ipmt++){
451          if(time[ipmt] !=0 /*&& badpmt[ipmt]==0 */&&  adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold )
452            {
453                //              if(TMath::Abs(time[ipmt])<TMath::Abs(besttimeC)) {
454              if(time[ipmt]<besttimeC){
455                besttimeC=time[ipmt]; //timeC
456                   pmtBestC=ipmt;
457              }
458            }
459        }
460        for ( Int_t ipmt=12; ipmt<24; ipmt++)
461          {
462            if(time[ipmt] != 0 /* && badpmt[ipmt]==0*/ && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold)
463              {
464                if(time[ipmt]<besttimeA) {
465                  //            if(TMath::Abs(time[ipmt])<TMath::Abs(besttimeA)) {
466                  besttimeA=time[ipmt]; //timeA
467                  pmtBestA=ipmt;
468                }
469              }
470          }
471        if(besttimeA < 999999) 
472          frecpoints->SetTimeBestA((besttimeA * channelWidth)- 1000.*fLatencyHPTDC + 1000.*fLatencyL1A - 1000.*fGRPdelays - fTimeMeanShift[1] ); 
473        //        frecpoints->SetTimeBestA((besttimeA * channelWidth- fTimeMeanShift[1])); 
474
475        if( besttimeC < 999999 ) 
476          frecpoints->SetTimeBestC((besttimeC * channelWidth)- 1000.*fLatencyHPTDC +1000.*fLatencyL1C - 1000.*fGRPdelays - fTimeMeanShift[2]);
477        // frecpoints->SetTimeBestC((besttimeC * channelWidth - fTimeMeanShift[2]));
478        AliDebug(5,Form(" pmtA %i besttimeA %f shift A %f ps, pmtC %i besttimeC %f shiftC %f ps",
479                        pmtBestA,besttimeA, fTimeMeanShift[1],
480                        pmtBestC,  besttimeC,fTimeMeanShift[2]));
481         if(besttimeA <999999 && besttimeC < 999999 ){
482           //     timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth + fLatencyL1A - fLatencyL1C;
483           timeclock = channelWidth * Float_t( besttimeA+besttimeC)/2. - 1000.*fLatencyHPTDC + 1000.*fLatencyL1 - 1000.*fGRPdelays - fTimeMeanShift[0] ;  
484          meanTime = (besttimeA+besttimeC-2.*Float_t(ref))/2.;
485          timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth ;
486          // timeclock = channelWidth * Float_t( besttimeA+besttimeC)/2. - fTimeMeanShift[0] ;  
487          vertex =  meanVertex - c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2; 
488         }
489       }  //if phys event       
490       AliDebug(10,Form("  timeDiff %f #channel,  meanTime %f #channel, TOFmean%f  vertex %f cm meanVertex %f online mean %i \n",timeDiff, meanTime,timeclock, vertex,meanVertex, onlineMean));
491       frecpoints->SetT0clock(timeclock);
492       frecpoints->SetVertex(vertex);
493       frecpoints->SetMeanTime(meanTime);
494       frecpoints->SetOnlineMean(Int_t(onlineMean));
495         // Set triggers
496       
497       Bool_t tr[5];
498       Int_t trchan[5]= {50,51,52,55,56};
499       for (Int_t i=0; i<5; i++) tr[i]=false; 
500       for (Int_t itr=0; itr<5; itr++) {
501         for (Int_t iHit=0; iHit<5; iHit++) 
502           {
503             Int_t trr=trchan[itr];
504             if( allData[trr][iHit] > 0) tr[itr]=true;
505           }
506       }
507       frecpoints->SetT0Trig(tr);
508    
509       //Set MPD
510       if(allData[53][0]>0 && allData[54][0]) 
511         frecpoints->SetMultA(allData[53][0]-allData[54][0]);
512         if(allData[105][0]>0 && allData[106][0]) 
513           frecpoints->SetMultC(allData[105][0]-allData[106][0]);
514         
515         
516     } // if (else )raw data
517   recTree->Fill();
518   if(frecpoints) delete frecpoints;
519 }
520   
521   
522   //____________________________________________________________
523   
524   void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
525   {
526
527   /***************************************************
528   Resonstruct digits to vertex position
529   ****************************************************/
530   
531   AliDebug(1,Form("Start FillESD T0"));
532   if(!pESD) {
533     AliError("No ESD Event");
534     return;
535   }
536   pESD ->SetT0spread(fTimeSigmaShift);
537  
538
539   Float_t channelWidth = fParam->GetChannelWidth() ;  
540   Float_t c = 0.0299792458; // cm/ps
541   Float_t currentVertex=0, shift=0;
542   Int_t ncont=-1;
543   const AliESDVertex* vertex = pESD->GetPrimaryVertex();
544   if (!vertex)        vertex = pESD->GetPrimaryVertexSPD();
545   if (!vertex)        vertex = pESD->GetPrimaryVertexTPC();
546   if (!vertex)        vertex = pESD->GetVertex();
547
548   if (vertex) {
549     AliDebug(2, Form("Got %s (%s) from ESD: %f", 
550                     vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
551     currentVertex = vertex->GetZ();
552     
553     ncont = vertex->GetNContributors();
554     //  cout<<"@@ spdver "<<spdver<<" ncont "<<ncont<<endl;
555     if(ncont>0 ) {
556       shift = currentVertex/c;
557     }
558   }
559   TTree *treeR = clustersTree;
560   
561   AliT0RecPoint* frecpoints= new AliT0RecPoint ();
562   if (!frecpoints) {
563     AliError("Reconstruct Fill ESD >> no recpoints found");
564     return;
565   }
566   
567   AliDebug(1,Form("Start FillESD T0"));
568   TBranch *brRec = treeR->GetBranch("T0");
569   if (brRec) {
570     brRec->SetAddress(&frecpoints);
571   }else{
572     AliError(Form("EXEC Branch T0 rec not found"));
573     return;
574   } 
575   
576   brRec->GetEntry(0);
577   Double32_t amp[24], time[24], ampQTC[24], timecorr[24];  
578   Double32_t* tcorr;
579   for(Int_t i=0; i<24; i++) 
580     amp[i]=time[i]=ampQTC[i]=timecorr[i]=0;
581
582
583   Double32_t timeClock[3];
584   Double32_t zPosition = frecpoints -> GetVertex();
585   Double32_t timeStart = frecpoints -> GetMeanTime();
586   timeClock[0] = frecpoints -> GetT0clock() ;
587   timeClock[1] = frecpoints -> GetBestTimeA() + shift;
588   timeClock[2] = frecpoints -> GetBestTimeC() - shift;
589
590   for ( Int_t i=0; i<24; i++) {
591     time[i] =  frecpoints -> GetTime(i); // ps to ns
592     //    if ( time[i] >1) {
593     if ( time[i] != 0) {
594       ampQTC[i] = frecpoints -> GetAmp(i);
595       amp[i] = frecpoints -> AmpLED(i);
596       AliDebug(1,Form("T0: %i  time %f  ampQTC %f ampLED %f \n", i, time[i], ampQTC[i], amp[i]));
597    }
598   }
599   Int_t trig= frecpoints ->GetT0Trig();
600   pESD->SetT0Trig(trig);
601   
602   pESD->SetT0zVertex(zPosition); //vertex Z position 
603
604   Double32_t multA=frecpoints ->GetMultA();
605   Double32_t multC=frecpoints ->GetMultC();
606   //  pESD->SetT0MultC(multC);        // multiplicity Cside 
607   //  pESD->SetT0MultA(multA);        // multiplicity Aside 
608   pESD->SetT0(multA); // for backward compatubility
609   pESD->SetT0clock(multC); // for backward compatubility
610
611   for(Int_t i=0; i<3; i++) 
612     pESD->SetT0TOF(i,timeClock[i]);   // interaction time (ns) 
613   pESD->SetT0time(time);         // best TOF on each PMT 
614   pESD->SetT0amplitude(ampQTC);     // number of particles(MIPs) on each PMT
615   
616   AliDebug(1,Form("T0: SPDshift %f Vertex %f (T0A+T0C)/2 %f #channels T0signal %f ns OrA %f ns OrC %f T0trig %i\n",shift, zPosition, timeStart, timeClock[0], timeClock[1], timeClock[2], trig));
617   
618   if (pESD) {
619     
620     AliESDfriend *fr = (AliESDfriend*)pESD->FindListObject("AliESDfriend");
621     if (fr) {
622       AliDebug(1, Form("Writing TZERO friend data to ESD tree"));
623
624       //     if (ncont>2) {
625         tcorr = fESDTZEROfriend->GetT0timeCorr();
626         for ( Int_t i=0; i<24; i++) {
627           if(i<12 && time[i]>1) timecorr[i] = tcorr[i] -  shift/channelWidth;
628           if(i>11 && time[i]>1) timecorr[i] = tcorr[i] +  shift/channelWidth;
629           if(time[i]>1)   AliDebug(10,Form("T0 friend : %i time %f  ampQTC %f ampLED %f \n", i, timecorr[i], ampQTC[i], amp[i]));
630         }
631         fESDTZEROfriend->SetT0timeCorr( timecorr) ;
632         fESDTZEROfriend->SetT0ampLEDminCFD(amp);
633         fESDTZEROfriend->SetT0ampQTC(ampQTC);
634         fr->SetTZEROfriend(fESDTZEROfriend);
635         //      }//
636     }
637   }
638      
639
640
641 } // vertex in 3 sigma
642
643
644
645
646
647