Bug fix in the order of the Ds cuts (Sadhana, Francesco)
[u/mrichter/AliRoot.git] / T0 / AliT0Reconstructor.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /* $Id$ */
18 /*********************************************************************
19  *  T0 reconstruction and filling ESD
20  *  - reconstruct mean time (interation time) 
21  *  - vertex position
22  *  -  multiplicity
23  ********************************************************************/
24
25 #include <AliESDEvent.h>
26 #include "AliLog.h"
27 #include "AliT0RecPoint.h"
28 #include "AliRawReader.h"
29 #include "AliT0RawReader.h"
30 #include "AliT0digit.h"
31 #include "AliT0Reconstructor.h"
32 #include "AliT0Parameters.h"
33 #include "AliT0Calibrator.h"
34 #include "AliESDfriend.h"
35 #include "AliESDTZEROfriend.h"
36 #include "AliLog.h"
37
38 #include <TArrayI.h>
39 #include <TGraph.h>
40 #include <TMath.h>
41 #include <Riostream.h>
42
43 ClassImp(AliT0Reconstructor)
44
45   AliT0Reconstructor:: AliT0Reconstructor(): AliReconstructor(),
46                                              fdZonA(0),
47                                              fdZonC(0),
48                                              fZposition(0),
49                                              fParam(NULL),
50                                              fAmpLEDrec(),
51                                              fQTC(0),
52                                              fAmpLED(0),
53                                              fCalib(),
54                                              fLatencyHPTDC(9000),
55                                              fLatencyL1(0),
56                                              fLatencyL1A(0),
57                                              fLatencyL1C(0),
58                                              fESDTZEROfriend(NULL)
59
60 {
61   //constructor
62
63   fParam = AliT0Parameters::Instance();
64   fParam->Init();
65  
66   for (Int_t i=0; i<24; i++){
67         TGraph* gr = fParam ->GetAmpLEDRec(i);
68         if (gr) fAmpLEDrec.AddAtAndExpand(gr,i) ; 
69           TGraph* gr1 = fParam ->GetAmpLED(i);
70           if (gr1) fAmpLED.AddAtAndExpand(gr1,i) ; 
71           TGraph* gr2 = fParam ->GetQTC(i);
72           if (gr2) fQTC.AddAtAndExpand(gr2,i) ;         
73   }
74
75   fLatencyL1 = fParam->GetLatencyL1();
76   fLatencyL1A = fParam->GetLatencyL1A();
77   fLatencyL1C = fParam->GetLatencyL1C();
78   fLatencyHPTDC = fParam->GetLatencyHPTDC();
79   AliDebug(10,Form(" LatencyL1 %f latencyL1A %f latencyL1C %f latencyHPTDC %f \n",fLatencyL1, fLatencyL1A, fLatencyL1C, fLatencyHPTDC));
80   
81   // fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1"));
82   //fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));
83
84
85   fCalib = new AliT0Calibrator();
86   fESDTZEROfriend = new AliESDTZEROfriend();
87
88 }
89
90 //_____________________________________________________________________________
91 void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const
92   
93 {
94   // T0 digits reconstruction
95   Int_t refAmp = GetRecoParam()->GetRefAmp();
96   
97   TArrayI * timeCFD = new TArrayI(24); 
98   TArrayI * timeLED = new TArrayI(24); 
99   TArrayI * chargeQT0 = new TArrayI(24); 
100   TArrayI * chargeQT1 = new TArrayI(24); 
101
102  
103   Float_t channelWidth = fParam->GetChannelWidth() ;  
104   Float_t meanVertex = fParam->GetMeanVertex();
105   Float_t c = 0.0299792; // cm/ps
106   Double32_t vertex = 9999999;
107   Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;
108
109   
110   AliDebug(1,Form("Start DIGITS reconstruction "));
111   
112
113   TBranch *brDigits=digitsTree->GetBranch("T0");
114   AliT0digit *fDigits = new AliT0digit() ;
115   if (brDigits) {
116     brDigits->SetAddress(&fDigits);
117   }else{
118     AliError(Form("EXEC Branch T0 digits not found"));
119      return;
120   }
121   
122   digitsTree->GetEvent(0);
123   digitsTree->GetEntry(0);
124   brDigits->GetEntry(0);
125   fDigits->GetTimeCFD(*timeCFD);
126   fDigits->GetTimeLED(*timeLED);
127   fDigits->GetQT0(*chargeQT0);
128   fDigits->GetQT1(*chargeQT1);
129   Int_t onlineMean =  fDigits->MeanTime();
130
131   Bool_t tr[5];
132   for (Int_t i=0; i<5; i++) tr[i]=false; 
133   
134   Double32_t besttimeA=999999;
135   Double32_t besttimeC=999999;
136   Int_t pmtBestA=99999;
137   Int_t pmtBestC=99999;
138   
139   AliT0RecPoint* frecpoints= new AliT0RecPoint ();
140   clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
141   
142   Float_t time[24], adc[24];
143   for (Int_t ipmt=0; ipmt<24; ipmt++) {
144     if(timeCFD->At(ipmt)>0 ){
145      if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)  
146         adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
147       else
148         adc[ipmt] = 0;
149       
150      // time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt],  timeCFD->At(ipmt)) ;
151              
152       Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
153       time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl),  timeCFD->At(ipmt) ) ;
154       AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
155                        ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
156
157       Double_t ampMip =((TGraph*)fAmpLED.At(ipmt))->Eval(sl);
158       Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
159       AliDebug(10,Form("  Amlitude in MIPS LED %f ,  QTC %f in channels %i\n ",ampMip,qtMip, adc[ipmt]));
160       
161       frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
162       frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam 
163       frecpoints->SetAmpLED(ipmt, Float_t(qtMip));
164       
165     }
166     else {
167       time[ipmt] = 0;
168       adc[ipmt] = 0;
169     }
170   }
171   
172   for (Int_t ipmt=0; ipmt<12; ipmt++){
173     if(time[ipmt] > 1 ) {
174       if(time[ipmt]<besttimeC){
175         besttimeC=time[ipmt]; //timeC
176         pmtBestC=ipmt;
177       }
178     }
179   }
180   for ( Int_t ipmt=12; ipmt<24; ipmt++){
181     if(time[ipmt] > 1) {
182       if(time[ipmt]<besttimeA) {
183         besttimeA=time[ipmt]; //timeA
184         pmtBestA=ipmt;}
185     }
186   }
187   if(besttimeA < 999999) {
188     frecpoints->SetTimeBestA(Int_t(besttimeA *channelWidth));
189     tr[1]=true;
190   }
191   if( besttimeC < 999999 ) {
192     frecpoints->SetTimeBestC(Int_t(besttimeC *channelWidth));
193     tr[2]=true;
194   }
195   AliDebug(10,Form(" besttimeA %f ch,  besttimeC %f ch",besttimeA, besttimeC));
196   if(besttimeA <999999 && besttimeC < 999999 ){
197     //    timeDiff = (besttimeC - besttimeA)*channelWidth;
198     timeDiff = (besttimeA - besttimeC)*channelWidth;
199     meanTime = (besttimeA + besttimeC)/2;// * channelWidth); 
200     timeclock = meanTime *channelWidth ;
201     vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
202     tr[0]=true; 
203   }
204   frecpoints->SetVertex(vertex);
205   frecpoints->SetMeanTime(meanTime);
206   frecpoints->SetT0clock(timeclock);
207   frecpoints->SetT0Trig(tr);
208
209   AliInfo(Form("T0 triggers %d %d %d %d %d",tr[0],tr[1],tr[2],tr[3],tr[4]));
210
211   //online mean
212   frecpoints->SetOnlineMean(Int_t(onlineMean));
213   AliDebug(10,Form("  timeDiff %i #channel,  meanTime %i #channel, vertex %f cm online mean %i timeclock %i ps",timeDiff, meanTime,vertex, Int_t(onlineMean), timeclock));
214   
215   
216
217    
218   
219   clustersTree->Fill();
220
221   delete timeCFD;
222   delete timeLED;
223   delete chargeQT0; 
224   delete chargeQT1; 
225 }
226
227
228 //_______________________________________________________________________
229
230 void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
231 {
232   // T0 raw ->
233   //
234   // reference amplitude and time ref. point from reco param
235
236   Int_t refAmp = GetRecoParam()->GetRefAmp();
237   Int_t refPoint = GetRecoParam()->GetRefPoint();
238
239   Int_t allData[110][5];
240   
241   Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
242   Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;
243   Float_t c = 29.9792458; // cm/ns
244   Double32_t vertex = 9999999;
245   Int_t onlineMean=0;
246   // Float_t meanVertex = fParam->GetMeanVertex();
247   Float_t meanVertex = 0;
248   for (Int_t i0=0; i0<105; i0++)
249     {
250       for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;   
251     }
252    
253   Double32_t besttimeA=9999999;
254   Double32_t besttimeC=9999999;
255   Int_t pmtBestA=99999;
256   Int_t pmtBestC=99999;
257    
258   AliT0RecPoint* frecpoints= new AliT0RecPoint ();
259   
260   recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
261    
262   AliDebug(10," before read data ");
263   AliT0RawReader myrawreader(rawReader);
264
265   UInt_t type =rawReader->GetType();
266
267   if (!myrawreader.Next())
268     AliDebug(1,Form(" no raw data found!!"));
269   else
270     {  
271       if(type == 7) {  //only physics 
272       for (Int_t i=0; i<105; i++) {
273         for (Int_t iHit=0; iHit<5; iHit++) 
274           {
275             allData[i][iHit] = myrawreader.GetData(i,iHit);
276           }
277       }
278       Int_t ref=0;
279       if (refPoint>0) 
280       ref = allData[refPoint][0]-5000;
281
282       Float_t channelWidth = fParam->GetChannelWidth() ;  
283       
284       //       Int_t meanT0 = fParam->GetMeanT0();
285        
286           
287           for (Int_t in=0; in<12; in++)  
288             {
289               timeCFD[in] = allData[in+1][0] ;
290               timeCFD[in+12] = allData[in+56+1][0] ;
291               timeLED[in] = allData[in+12+1][0] ;
292               timeLED[in+12] = allData[in+68+1][0] ;
293               AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
294                                 in, timeCFD[in],timeCFD[in+12],timeLED[in], 
295                                 timeLED[in+12]));   
296             }
297           
298           for (Int_t in=0; in<12;  in++)
299             {
300               chargeQT0[in]=allData[2*in+25][0];
301               chargeQT1[in]=allData[2*in+26][0];
302             }
303           
304            for (Int_t in=12; in<24;  in++)
305              {
306                chargeQT0[in]=allData[2*in+57][0];
307                chargeQT1[in]=allData[2*in+58][0];
308              }
309            
310            //    } //cosmic with physics event
311        for (Int_t in=0; in<24; in++)  
312          AliDebug(10, Form(" readed Raw %i %i %i %i %i",
313                            in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
314         onlineMean = allData[49][0];       
315        
316         Double32_t time[24], adc[24],  noncalibtime[24];
317        for (Int_t ipmt=0; ipmt<24; ipmt++) {
318          if(timeCFD[ipmt]>0 && timeLED[ipmt]>0){
319            //for simulated data
320              //for physics  data
321            if(( chargeQT0[ipmt] - chargeQT1[ipmt])>0)  {
322              adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
323            }
324            else
325              adc[ipmt] = 0;
326            
327
328            //      time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD[ipmt] ) ;
329            
330            Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
331            time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
332            AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
333                             ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
334            Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
335            Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
336            AliDebug(10,Form("  Amlitude in MIPS LED %f ; QTC %f;  in channels %i\n ",ampMip,qtMip, adc[ipmt]));
337            //bad peak removing
338            if(sl<540) {
339              frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
340              // frecpoints->SetTime(ipmt,Double32_t(timeCFD[ipmt]));
341              frecpoints->SetAmpLED(ipmt, Double32_t( qtMip)); //for cosmic &pp beam 
342              frecpoints->SetAmp(ipmt, Double32_t(ampMip));           
343              noncalibtime[ipmt]= Double32_t (timeCFD[ipmt]);
344            }
345          }
346          else {
347            time[ipmt] = 0;
348            adc[ipmt] = 0;
349            noncalibtime[ipmt] = 0;
350          }
351        }
352        fESDTZEROfriend->SetT0timeCorr(noncalibtime) ;     
353        for (Int_t ipmt=0; ipmt<12; ipmt++){
354          if(time[ipmt] > 1 
355             && (timeLED[ipmt] - timeCFD[ipmt])<540 )
356            {
357              if(time[ipmt]<besttimeC){
358                besttimeC=time[ipmt]; //timeC
359                pmtBestC=ipmt;
360              }
361            }
362        }
363        for ( Int_t ipmt=12; ipmt<24; ipmt++)
364          {
365            if(time[ipmt] > 1 
366               && (timeLED[ipmt] - timeCFD[ipmt])<540 )
367              {
368                if(time[ipmt]<besttimeA) {
369                  besttimeA=time[ipmt]; //timeA
370                  pmtBestA=ipmt;
371                }
372              }
373          }
374        if(besttimeA < 999999) 
375          frecpoints->SetTimeBestA(besttimeA * channelWidth - 1000.*fLatencyHPTDC + 1000.*fLatencyL1A);
376        if( besttimeC < 999999 ) 
377          frecpoints->SetTimeBestC(besttimeC * channelWidth - 1000.*fLatencyHPTDC +1000.*fLatencyL1C);
378        AliDebug(10,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f ps",
379                        pmtBestA,besttimeA, pmtBestC,  besttimeC));
380         if(besttimeA <999999 && besttimeC < 999999 ){
381          timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth + fLatencyL1A - fLatencyL1C;
382          timeclock = channelWidth * Float_t( besttimeA+besttimeC)/2. - 1000.*fLatencyHPTDC + 1000.*fLatencyL1;  
383          meanTime = (besttimeA+besttimeC-2.*Float_t(ref))/2.;
384          vertex =  meanVertex - c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2; 
385         }
386       }  //if phys event       
387       AliDebug(5,Form("  timeDiff %f #channel,  meanTime %f #channel, TOFmean%f  vertex %f cm meanVertex %f online mean %i \n",timeDiff, meanTime,timeclock, vertex,meanVertex, onlineMean));
388       frecpoints->SetT0clock(timeclock);
389       frecpoints->SetVertex(vertex);
390       frecpoints->SetMeanTime(meanTime);
391       frecpoints->SetOnlineMean(Int_t(onlineMean));
392         // Set triggers
393       
394       Bool_t tr[5];
395       Int_t trchan[5]= {50,51,52,55,56};
396       for (Int_t i=0; i<5; i++) tr[i]=false; 
397       for (Int_t itr=0; itr<5; itr++) {
398         if(allData[trchan[itr]][0]>0) tr[itr]=true;
399         frecpoints->SetT0Trig(tr);
400       }
401     } // if (else )raw data
402   recTree->Fill();
403   if(frecpoints) delete frecpoints;
404 }
405   
406   
407   //____________________________________________________________
408   
409   void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
410   {
411
412   /***************************************************
413   Resonstruct digits to vertex position
414   ****************************************************/
415   
416   AliDebug(1,Form("Start FillESD T0"));
417   Float_t channelWidth = fParam->GetChannelWidth() ;  
418   Float_t c = 0.0299792458; // cm/ps
419   Float_t currentVertex=0, shift=0;
420   Int_t ncont=0;
421   const AliESDVertex* vertex = pESD->GetPrimaryVertex();
422   if (!vertex)        vertex = pESD->GetPrimaryVertexSPD();
423   if (!vertex)        vertex = pESD->GetPrimaryVertexTPC();
424   if (!vertex)        vertex = pESD->GetVertex();
425
426   if (vertex) {
427     AliDebug(2, Form("Got %s (%s) from ESD: %f", 
428                     vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
429     currentVertex = vertex->GetZ();
430     
431     ncont = vertex->GetNContributors();
432     // cout<<" spdver "<<spdver<<" ncont "<<ncont<<endl;
433     if(ncont>2 ) {
434       shift = currentVertex/c;
435       //          cout<<" vertex shif "<<shift<<" vertex "<<spdver<<" IsFromVertexer3D  "<<fverSPD->IsFromVertexer3D()<<endl;
436     }
437   }
438   TTree *treeR = clustersTree;
439   
440   AliT0RecPoint* frecpoints= new AliT0RecPoint ();
441   if (!frecpoints) {
442     AliError("Reconstruct Fill ESD >> no recpoints found");
443     return;
444   }
445   
446   AliDebug(1,Form("Start FillESD T0"));
447   TBranch *brRec = treeR->GetBranch("T0");
448   if (brRec) {
449     brRec->SetAddress(&frecpoints);
450   }else{
451     AliError(Form("EXEC Branch T0 rec not found"));
452     return;
453   } 
454   
455   brRec->GetEntry(0);
456   Double32_t amp[24], time[24], ampQTC[24], timecorr[24];  
457   Double32_t* tcorr;
458   for(Int_t i=0; i<24; i++) 
459     amp[i]=time[i]=ampQTC[i]=timecorr[i]=0;
460
461   Double32_t timeClock[3];
462   Double32_t zPosition = frecpoints -> GetVertex();
463   Double32_t timeStart = frecpoints -> GetMeanTime();
464   timeClock[0] = frecpoints -> GetT0clock() ;
465   timeClock[1] = frecpoints -> GetBestTimeA() + shift;
466   timeClock[2] = frecpoints -> GetBestTimeC() - shift;
467   for ( Int_t i=0; i<24; i++) {
468     time[i] =  frecpoints -> GetTime(i); // ps to ns
469     if ( time[i] >1) {
470       amp[i] = frecpoints -> GetAmp(i);
471       ampQTC[i] = frecpoints -> AmpLED(i);
472     }
473   }
474   Int_t trig= frecpoints ->GetT0Trig();
475   pESD->SetT0Trig(trig);
476   
477   pESD->SetT0zVertex(zPosition); //vertex Z position 
478   pESD->SetT0(timeStart);        // interaction time 
479   for(Int_t i=0; i<3; i++) 
480     pESD->SetT0TOF(i,timeClock[i]);   // interaction time (ns) 
481   pESD->SetT0time(time);         // best TOF on each PMT 
482   pESD->SetT0amplitude(amp);     // number of particles(MIPs) on each PMT
483   
484   AliDebug(1,Form("T0: Vertex %f (T0A+T0C)/2 %f #channels T0signal %f ns OrA %f ns OrC %f T0trig %i\n",zPosition, timeStart, timeClock[0], timeClock[1], timeClock[2], trig));
485   
486   if (pESD) {
487     
488     AliESDfriend *fr = (AliESDfriend*)pESD->FindListObject("AliESDfriend");
489     if (fr) {
490       AliDebug(1, Form("Writing TZERO friend data to ESD tree"));
491
492       if (ncont>2) {
493         tcorr = fESDTZEROfriend->GetT0timeCorr();
494         for ( Int_t i=0; i<24; i++) {
495           timecorr[i]=tcorr[i];
496           if(i<12 && time[i]>1) timecorr[i] -=  shift/channelWidth;
497           if(i>11 && time[i]>1) timecorr[i] +=  shift/channelWidth;
498         }
499         fESDTZEROfriend->SetT0timeCorr( timecorr) ;
500         fESDTZEROfriend->SetT0ampLEDminCFD(amp);
501         fESDTZEROfriend->SetT0ampQTC(ampQTC);
502         fr->SetTZEROfriend(fESDTZEROfriend);
503       }//
504     }
505   }
506      
507
508
509 } // vertex in 3 sigma
510
511
512
513
514
515