uuhhhhaaa - what did I change? Added the run number and fixed coding conventions ;-)
[u/mrichter/AliRoot.git] / TOF / AliTOFcalib.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 /*
17 $Log$
18 Revision 1.9  2006/04/20 22:30:50  hristov
19 Coding conventions (Annalisa)
20
21 Revision 1.8  2006/04/16 22:29:05  hristov
22 Coding conventions (Annalisa)
23
24 Revision 1.7  2006/04/16 20:12:46  hristov
25 Removing memory leak in case of cached CDB entries
26
27 Revision 1.6  2006/04/11 15:28:32  hristov
28 Checks on cache status before deleting calibration objects (A.Colla)
29
30 Revision 1.5  2006/04/05 08:35:38  hristov
31 Coding conventions (S.Arcelli, C.Zampolli)
32
33 Revision 1.4  2006/03/31 11:26:46  arcelli
34  changing CDB Ids according to standard convention
35
36 Revision 1.3  2006/03/28 14:57:02  arcelli
37 updates to handle new V5 geometry & some re-arrangements
38
39 Revision 1.2  2006/02/13 17:22:26  arcelli
40 just Fixing Log info
41
42 Revision 1.1  2006/02/13 16:10:48  arcelli
43 Add classes for TOF Calibration (C.Zampolli)
44
45 author: Chiara Zampolli, zampolli@bo.infn.it
46 */  
47
48 ///////////////////////////////////////////////////////////////////////////////
49 //                                                                           //
50 // class for TOF calibration                                                 //
51 //                                                                           //
52 ///////////////////////////////////////////////////////////////////////////////
53
54 #include "TF1.h"
55 #include "TFile.h"
56 #include "TH1F.h"
57 #include "TH2F.h"
58 #include "TList.h"
59 #include "TROOT.h"
60 #include "TStyle.h"
61
62 #include "AliCDBEntry.h"
63 #include "AliCDBId.h"
64 #include "AliCDBManager.h"
65 #include "AliCDBMetaData.h"
66 #include "AliESDtrack.h"
67 #include "AliESD.h"
68 #include "AliLog.h"
69
70 #include "AliTOFCal.h"
71 #include "AliTOFcalibESD.h"
72 #include "AliTOFcalib.h"
73 #include "AliTOFChannel.h"
74 #include "AliTOFGeometryV5.h"
75 #include "AliTOFGeometry.h"
76
77 extern TROOT *gROOT;
78 extern TStyle *gStyle;
79
80 ClassImp(AliTOFcalib)
81
82 const Int_t AliTOFcalib::fgkchannel = 5000;
83 //_______________________________________________________________________
84 AliTOFcalib::AliTOFcalib():
85   TTask("AliTOFcalib",""),
86   fNChannels(-1),
87   fNSector(-1),
88   fNPlate(-1),
89   fNStripA(-1),
90   fNStripB(-1),
91   fNStripC(-1),
92   fNpadZ(-1),
93   fNpadX(-1),
94   fNevents(0),
95   fESDsel(0x0),
96   fArrayToT(0x0),
97   fArrayTime(0x0),
98   fTOFCal(0x0),
99   fTOFSimCal(0x0),
100   fTOFSimToT(0x0)
101
102   //TOF Calibration Class ctor
103   fArrayToT = 0x0;
104   fArrayTime = 0x0;
105   fESDsel = 0x0;
106   AliTOFGeometry *geom=new AliTOFGeometryV5();
107   AliInfo("V5 TOF Geometry is taken as the default");
108   fNSector = geom->NSectors();
109   fNPlate  = geom->NPlates();
110   fNStripA = geom->NStripA();
111   fNStripB = geom->NStripB();
112   fNStripC = geom->NStripC();
113   fNpadZ = geom->NpadZ();
114   fNpadX = geom->NpadX();
115   fNChannels = fNSector*(2*(fNStripC+fNStripB)+fNStripA)*fNpadZ*fNpadX; //generalized version
116   fTOFCal = 0x0;
117   fTOFSimCal = 0x0;
118   fTOFSimToT = 0x0;
119   delete geom;
120 }
121 //_______________________________________________________________________
122 AliTOFcalib::AliTOFcalib(AliTOFGeometry *geom):
123   TTask("AliTOFcalib",""),
124   fNChannels(-1),
125   fNSector(-1),
126   fNPlate(-1),
127   fNStripA(-1),
128   fNStripB(-1),
129   fNStripC(-1),
130   fNpadZ(-1),
131   fNpadX(-1),
132   fNevents(0),
133   fESDsel(0x0),
134   fArrayToT(0x0),
135   fArrayTime(0x0),
136   fTOFCal(0x0),
137   fTOFSimCal(0x0),
138   fTOFSimToT(0x0)
139
140   //TOF Calibration Class ctor, taking the TOF geometry as input
141   fArrayToT = 0x0;
142   fArrayTime = 0x0;
143   fESDsel = 0x0;
144   fNSector = geom->NSectors();
145   fNPlate  = geom->NPlates();
146   fNStripA = geom->NStripA();
147   fNStripB = geom->NStripB();
148   fNStripC = geom->NStripC();
149   fNpadZ = geom->NpadZ();
150   fNpadX = geom->NpadX();
151   fNChannels = fNSector*(2*(fNStripC+fNStripB)+fNStripA)*fNpadZ*fNpadX; //generalized version
152   fTOFCal = 0x0;
153   fTOFSimCal = 0x0;
154   fTOFSimToT = 0x0;
155 }
156 //____________________________________________________________________________ 
157
158 AliTOFcalib::AliTOFcalib(const AliTOFcalib & calib):
159   TTask("AliTOFcalib",""),
160   fNChannels(-1),
161   fNSector(-1),
162   fNPlate(-1),
163   fNStripA(-1),
164   fNStripB(-1),
165   fNStripC(-1),
166   fNpadZ(-1),
167   fNpadX(-1),
168   fNevents(0),
169   fESDsel(0x0),
170   fArrayToT(0x0),
171   fArrayTime(0x0),
172   fTOFCal(0x0),
173   fTOFSimCal(0x0),
174   fTOFSimToT(0x0)
175 {
176   //TOF Calibration Class copy ctor
177   fNSector = calib.fNSector;
178   fNPlate = calib.fNPlate;
179   fNStripA = calib.fNStripA;
180   fNStripB = calib.fNStripB;
181   fNStripC = calib.fNStripC;
182   fNpadZ = calib.fNpadZ;
183   fNpadX = calib.fNpadX;
184   fNChannels = calib.fNChannels;
185   fArrayToT = calib.fArrayToT;
186   fArrayTime = calib.fArrayTime;
187   fTOFCal=calib.fTOFCal;
188   fTOFSimCal = calib.fTOFSimCal;
189   fTOFSimToT=calib.fTOFSimToT;
190 }
191
192 //____________________________________________________________________________ 
193
194 AliTOFcalib& AliTOFcalib::operator=(const AliTOFcalib &calib)
195 {
196   //TOF Calibration Class assignment operator
197   this->fNSector = calib.fNSector;
198   this->fNPlate = calib.fNPlate;
199   this->fNStripA = calib.fNStripA;
200   this->fNStripB = calib.fNStripB;
201   this->fNStripC = calib.fNStripC;
202   this->fNpadZ = calib.fNpadZ;
203   this->fNpadX = calib.fNpadX;
204   this->fNChannels = calib.fNChannels;
205   this->fArrayToT = calib.fArrayToT;
206   this->fArrayTime = calib.fArrayTime;
207   this->fTOFCal=calib.fTOFCal;
208   this->fTOFSimCal = calib.fTOFSimCal;
209   this->fTOFSimToT=calib.fTOFSimToT;
210   return *this;
211 }
212
213 //____________________________________________________________________________ 
214
215 AliTOFcalib::~AliTOFcalib()
216 {
217   //TOF Calibration Class dtor
218   delete fArrayToT;
219   delete fArrayTime;
220
221   if(!(AliCDBManager::Instance()->GetCacheFlag())){ // CDB objects must NOT be deleted if cache is active!
222         delete fTOFCal;
223         delete fTOFSimCal;
224   }
225   
226   delete fESDsel;
227 }
228 //__________________________________________________________________________
229
230 TF1* AliTOFcalib::SetFitFunctions(TH1F *histo)
231 {
232   //Define Fit Functions for Slewing Correction
233   TF1 * fpol[3];
234   const Int_t knbins = histo->GetNbinsX();
235   Float_t delta = histo->GetBinWidth(1);  //all the bins have the same width
236   Double_t max = histo->GetBinLowEdge(knbins)+delta;
237   max = 15;
238   fpol[0]=new TF1("poly3","pol3",5,max);
239   fpol[1]=new TF1("poly4","pol4",5,max);
240   fpol[2]=new TF1("poly5","pol5",5,max);
241   char npoly[10];
242   Double_t chi[3]={1E6,1E6,1E6};
243   Int_t ndf[3]={-1,-1,-1};
244   Double_t nchi[3]={1E6,1E6,1E6};
245   Double_t bestchi=1E6;
246   TF1 * fGold=0x0;
247   Int_t nonzero =0;
248   Int_t numberOfpar =0;
249   for (Int_t j=0; j<knbins; j++){
250     if (histo->GetBinContent(j)!=0) {
251       nonzero++;
252     }
253   }
254   Int_t norderfit = 0;
255   if (nonzero<=4) {
256     AliError(" Too few points in the histo. No fit performed.");
257     return 0x0;
258   }
259   else if (nonzero<=5) {
260     norderfit = 3;
261     AliInfo(" Only 3rd order polynomial fit possible.");
262   }
263   else if (nonzero<=6) {
264     norderfit = 4;
265     AliInfo(" Only 3rd and 4th order polynomial fit possible.");
266   }
267   else {
268     norderfit = 5;
269     AliInfo(" All 3rd, 4th and 5th order polynomial fit possible.");
270   }
271   for (Int_t ifun=norderfit-3;ifun<norderfit-2;ifun++){
272     sprintf(npoly,"poly%i",ifun+3);
273     histo->Fit(npoly, "ERN", " ", 5.,14.);
274     chi[ifun] = fpol[ifun]->GetChisquare();
275     ndf[ifun] = fpol[ifun]->GetNDF();
276     nchi[ifun] = (Double_t)chi[ifun]/ndf[ifun];
277     if (nchi[ifun]<bestchi) {
278       bestchi=nchi[ifun];
279       fGold = fpol[ifun];
280       numberOfpar = fGold->GetNpar();
281     }
282   }
283   fGold=fpol[2];  //Gold fit function set to pol5 in any case
284   histo->Fit(fGold,"ER " ," ",5.,15.);
285   return fGold;
286 }
287 //____________________________________________________________________________
288
289 void AliTOFcalib::SelectESD(AliESD *event) 
290 {
291   //track selection for Calibration
292   Float_t lowerMomBound=0.8; // [GeV/c] default value Pb-Pb
293   Float_t upperMomBound=1.8 ; // [GeV/c] default value Pb-Pb
294   Int_t ntrk =0;
295   Int_t ngoodtrkfinalToT = 0;
296   ntrk=event->GetNumberOfTracks();
297   fESDsel = new TObjArray(ntrk);
298   fESDsel->SetOwner();
299   TObjArray  uCdatatemp(ntrk);
300   Int_t ngoodtrk = 0;
301   Int_t ngoodtrkfinal = 0;
302   Float_t mintime =1E6;
303   for (Int_t itrk=0; itrk<ntrk; itrk++) {
304     AliESDtrack *t=event->GetTrack(itrk);
305     //track selection: reconstrution to TOF:
306     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
307       continue;
308     }
309     //IsStartedTimeIntegral
310     if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
311       continue;
312     }
313     Double_t time=t->GetTOFsignal();    
314     time*=1.E-3; // tof given in nanoseconds
315     if(time <= mintime)mintime=time;
316     Double_t mom=t->GetP();
317     if (!(mom<=upperMomBound && mom>=lowerMomBound))continue;
318     UInt_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
319     if(assignedTOFcluster==0){ // not matched
320       continue;
321     }
322     AliTOFcalibESD *unc = new AliTOFcalibESD;
323     unc->CopyFromAliESD(t);
324     Double_t c1[15]; 
325     unc->GetExternalCovariance(c1);
326     uCdatatemp.Add(unc);
327     ngoodtrk++;
328   }
329   for (Int_t i = 0; i < ngoodtrk ; i ++){
330     AliTOFcalibESD *unc = (AliTOFcalibESD*)uCdatatemp.At(i);
331     if((unc->GetTOFsignal()-mintime*1.E3)<5.E3){
332       fESDsel->Add(unc);
333       ngoodtrkfinal++;
334       ngoodtrkfinalToT++;
335     }
336   }
337   fESDsel->Sort();
338 }
339 //_____________________________________________________________________________
340
341 void AliTOFcalib::CombESDId()
342 {
343   //track PID for calibration
344   Float_t t0offset=0;
345   Float_t loffset=0;
346   Int_t   ntracksinset=6;
347   Float_t exptof[6][3];
348   Float_t momentum[6]={0.,0.,0.,0.,0.,0.};
349   Int_t   assparticle[6]={3,3,3,3,3,3};
350   Float_t massarray[3]={0.13957,0.493677,0.9382723};
351   Float_t timeofflight[6]={0.,0.,0.,0.,0.,0.};
352   Float_t beta[6]={0.,0.,0.,0.,0.,0.};
353   Float_t texp[6]={0.,0.,0.,0.,0.,0.};
354   Float_t sqMomError[6]={0.,0.,0.,0.,0.,0.};
355   Float_t sqTrackError[6]={0.,0.,0.,0.,0.,0.};
356   Float_t tracktoflen[6]={0.,0.,0.,0.,0.,0.};
357   Float_t timeResolution   = 0.90e-10; // 90 ps by default      
358   Float_t timeresolutioninns=timeResolution*(1.e+9); // convert in [ns]
359   Float_t timezero[6]={0.,0.,0.,0.,0.,0.};
360   Float_t weightedtimezero[6]={0.,0.,0.,0.,0.,0.};
361   Float_t besttimezero[6]={0.,0.,0.,0.,0.,0.};
362   Float_t bestchisquare[6]={0.,0.,0.,0.,0.,0.};
363   Float_t bestweightedtimezero[6]={0.,0.,0.,0.,0.,0.};
364   Float_t bestsqTrackError[6]={0.,0.,0.,0.,0.,0.};
365
366   Int_t nelements = fESDsel->GetEntries();
367   Int_t nset= (Int_t)(nelements/ntracksinset);
368   for (Int_t i=0; i< nset; i++) {   
369
370     AliTOFcalibESD **unc=new AliTOFcalibESD*[ntracksinset];
371     for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
372       Int_t index = itrk+i*ntracksinset;
373       AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(index);
374       unc[itrk]=element;
375     }
376     
377     for (Int_t j=0; j<ntracksinset; j++) {
378       AliTOFcalibESD *element=unc[j];
379       Double_t mom=element->GetP();
380       Double_t time=element->GetTOFsignal()*1.E-3; // in ns     
381       Double_t exptime[10]; 
382       element->GetIntegratedTimes(exptime);
383       Double_t toflen=element->GetIntegratedLength()/100.;  // in m
384       timeofflight[j]=time+t0offset;
385       tracktoflen[j]=toflen+loffset;
386       exptof[j][0]=exptime[2]*1.E-3+0.005;
387       exptof[j][1]=exptime[3]*1.E-3+0.005;
388       exptof[j][2]=exptime[4]*1.E-3+0.005;
389       momentum[j]=mom;
390     }
391     Float_t t0best=999.;
392     Float_t et0best=999.;
393     Float_t chisquarebest=999.;
394     for (Int_t i1=0; i1<3;i1++) {
395       beta[0]=momentum[0]/sqrt(massarray[i1]*massarray[i1]+momentum[0]*momentum[0]);
396       texp[0]=exptof[0][i1];
397       for (Int_t i2=0; i2<3;i2++) { 
398         beta[1]=momentum[1]/sqrt(massarray[i2]*massarray[i2]+momentum[1]*momentum[1]);
399         texp[1]=exptof[1][i2];
400         for (Int_t i3=0; i3<3;i3++) {
401           beta[2]=momentum[2]/sqrt(massarray[i3]*massarray[i3]+momentum[2]*momentum[2]);
402           texp[2]=exptof[2][i3];
403           for (Int_t i4=0; i4<3;i4++) {
404             beta[3]=momentum[3]/sqrt(massarray[i4]*massarray[i4]+momentum[3]*momentum[3]);
405             texp[3]=exptof[3][i4];
406             
407             for (Int_t i5=0; i5<3;i5++) {
408               beta[4]=momentum[4]/sqrt(massarray[i5]*massarray[i5]+momentum[4]*momentum[4]);
409               texp[4]=exptof[4][i5];
410               for (Int_t i6=0; i6<3;i6++) {
411                 beta[5]=momentum[5]/sqrt(massarray[i6]*massarray[i6]+momentum[5]*momentum[5]);
412                 texp[5]=exptof[5][i6];
413         
414                 Float_t sumAllweights=0.;
415                 Float_t meantzero=0.;
416                 Float_t emeantzero=0.;
417                 
418                 for (Int_t itz=0; itz<ntracksinset;itz++) {
419                   sqMomError[itz]=
420                     ((1.-beta[itz]*beta[itz])*0.025)*
421                     ((1.-beta[itz]*beta[itz])*0.025)*
422                     (tracktoflen[itz]/
423                      (0.299792*beta[itz]))*
424                     (tracktoflen[itz]/
425                      (0.299792*beta[itz])); 
426                   sqTrackError[itz]=
427                     (timeresolutioninns*
428                      timeresolutioninns
429                      +sqMomError[itz]); 
430                   
431                   timezero[itz]=texp[itz]-timeofflight[itz];                
432                   weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
433                   sumAllweights+=1./sqTrackError[itz];
434                   meantzero+=weightedtimezero[itz];
435                   
436                 } // end loop for (Int_t itz=0; itz<15;itz++)
437                 
438                 meantzero=meantzero/sumAllweights; // it is given in [ns]
439                 emeantzero=sqrt(1./sumAllweights); // it is given in [ns]
440                 
441                 // calculate chisquare
442                 
443                 Float_t chisquare=0.;           
444                 for (Int_t icsq=0; icsq<ntracksinset;icsq++) {
445                   chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
446                 } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
447                 //              cout << " chisquare " << chisquare << endl;
448                 
449                 Int_t npion=0;
450                 if(i1==0)npion++;
451                 if(i2==0)npion++;
452                 if(i3==0)npion++;
453                 if(i4==0)npion++;
454                 if(i5==0)npion++;
455                 if(i6==0)npion++;
456                 
457                 if(chisquare<=chisquarebest  && ((Float_t) npion/ ((Float_t) ntracksinset)>0.3)){
458                   //  if(chisquare<=chisquarebest){
459                   
460                   for(Int_t iqsq = 0; iqsq<ntracksinset; iqsq++) {
461                     bestsqTrackError[iqsq]=sqTrackError[iqsq]; 
462                     besttimezero[iqsq]=timezero[iqsq]; 
463                     bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
464                     bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq]; 
465                   }
466                   
467                   assparticle[0]=i1;
468                   assparticle[1]=i2;
469                   assparticle[2]=i3;
470                   assparticle[3]=i4;
471                   assparticle[4]=i5;
472                   assparticle[5]=i6;
473                   
474                   chisquarebest=chisquare;
475                   t0best=meantzero;
476                   et0best=emeantzero;
477                 } // close if(dummychisquare<=chisquare)
478               } // end loop on i6
479             } // end loop on i5
480           } // end loop on i4
481         } // end loop on i3
482       } // end loop on i2
483     } // end loop on i1
484
485
486     Float_t confLevel=999;
487     if(chisquarebest<999.){
488       Double_t dblechisquare=(Double_t)chisquarebest;
489       confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinset-1); 
490     }
491     // assume they are all pions for fake sets
492     if(confLevel<0.01 || confLevel==999. ){
493       for (Int_t itrk=0; itrk<ntracksinset; itrk++)assparticle[itrk]=0;
494     }
495     for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
496       Int_t index = itrk+i*ntracksinset;
497       AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(index);
498       element->SetCombID(assparticle[itrk]);
499     }
500   }
501 }
502
503 //_____________________________________________________________________________
504
505 void AliTOFcalib::CalibrateESD(){
506   //Calibrate selected ESD times
507   Int_t nelements = fESDsel->GetEntries();
508   Int_t *number=new Int_t[fNChannels];
509   fArrayToT = new AliTOFArray(fNChannels);
510   fArrayTime = new AliTOFArray(fNChannels);
511   for (Int_t i=0; i<fNChannels; i++){
512     number[i]=0;
513     fArrayToT->AddArray(i, new TArrayF(fgkchannel));
514     TArrayF * parrToT = fArrayToT->GetArray(i);
515     TArrayF & refaToT = * parrToT;
516     fArrayTime->AddArray(i, new TArrayF(fgkchannel));
517     TArrayF * parrTime = fArrayToT->GetArray(i);
518     TArrayF & refaTime = * parrTime;
519     for (Int_t j = 0;j<AliTOFcalib::fgkchannel;j++){
520       refaToT[j]=0.;      //ToT[i][j]=j;
521       refaTime[j]=0.;      //Time[i][j]=j;
522     }
523   }
524   
525   for (Int_t i=0; i< nelements; i++) {
526     AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(i);
527     Int_t ipid = element->GetCombID();
528     Double_t etime = 0;   //expected time
529     Double_t expTime[10]; 
530     element->GetIntegratedTimes(expTime);
531     if (ipid == 0) etime = expTime[2]*1E-3; //ns
532     else if (ipid == 1) etime = expTime[3]*1E-3; //ns
533     else if (ipid == 2) etime = expTime[4]*1E-3; //ns
534     else AliError("No pid from combinatorial algo for this track");
535     Double_t mtime = (Double_t)element->GetTOFsignal()*1E-3;  //measured time
536     Double_t mToT = (Double_t) element->GetToT();  //measured ToT, ns
537     //select the correspondent channel with its simulated ToT spectrum
538     //summing up everything, index = 0 for all channels:
539     Int_t index = element->GetTOFCalChannel();
540     Int_t index2 = number[index];
541     TArrayF * parrToT = fArrayToT->GetArray(index);
542     TArrayF & refaToT = * parrToT;
543     refaToT[index2] = (Float_t)mToT;
544     TArrayF * parrTime = fArrayTime->GetArray(index);
545     TArrayF & refaTime = * parrTime;
546     refaTime[index2] = (Float_t)(mtime-etime);
547     number[index]++;
548   }
549
550   for (Int_t i=0;i<1;i++){
551     TH1F * hProf = Profile(i);
552     TF1* fGold = SetFitFunctions(hProf);
553     Int_t nfpar = fGold->GetNpar();
554     Float_t par[6];    
555     for(Int_t kk=0;kk<6;kk++){
556       par[kk]=0;
557     }
558     for (Int_t kk = 0; kk< nfpar; kk++){
559       par[kk]=fGold->GetParameter(kk);
560     }
561     if (!fTOFCal) {
562       AliTOFGeometry *geom=new AliTOFGeometryV5();
563       fTOFCal = new AliTOFCal(geom);
564       fTOFCal->CreateArray();
565       delete geom;
566     }
567     AliTOFChannel * calChannel = fTOFCal->GetChannel(i);
568     calChannel->SetSlewPar(par);
569   }
570   delete[] number;
571 }
572
573 //___________________________________________________________________________
574
575 TH1F* AliTOFcalib::Profile(Int_t ich)
576 {
577   //Prepare histograms for Slewing Correction
578   const Int_t knbinToT = 650;
579   Int_t nbinTime = 400;
580   Float_t minTime = -10.5; //ns
581   Float_t maxTime = 10.5; //ns
582   Float_t minToT = 7.5; //ns
583   Float_t maxToT = 40.; //ns
584   Float_t deltaToT = (maxToT-minToT)/knbinToT;
585   Double_t mTime[knbinToT+1],mToT[knbinToT+1],meanTime[knbinToT+1], meanTime2[knbinToT+1],vToT[knbinToT+1], vToT2[knbinToT+1],meanToT[knbinToT+1],meanToT2[knbinToT+1],vTime[knbinToT+1],vTime2[knbinToT+1],xlow[knbinToT+1],sigmaTime[knbinToT+1];
586   Int_t n[knbinToT+1], nentrx[knbinToT+1];
587   Double_t sigmaToT[knbinToT+1];
588   for (Int_t i = 0; i < knbinToT+1 ; i++){
589     mTime[i]=0;
590     mToT[i]=0;
591     n[i]=0;
592     meanTime[i]=0;
593     meanTime2[i]=0;
594     vToT[i]=0;
595     vToT2[i]=0;
596     meanToT[i]=0;
597     meanToT2[i]=0;
598     vTime[i]=0;
599     vTime2[i]=0;
600     xlow[i]=0;
601     sigmaTime[i]=0;
602     sigmaToT[i]=0;
603     n[i]=0;
604     nentrx[i]=0;
605   }
606   TH2F* hSlewing = new TH2F("hSlewing", "hSlewing", knbinToT, minToT, maxToT, nbinTime, minTime, maxTime);
607   TArrayF * parrToT = fArrayToT->GetArray(ich);
608   TArrayF & refaToT = * parrToT;
609   TArrayF * parrTime = fArrayTime->GetArray(ich);
610   TArrayF & refaTime = * parrTime;
611   for (Int_t j = 0; j < AliTOFcalib::fgkchannel; j++){
612     if (refaToT[j] == 0) continue; 
613     Int_t nx = (Int_t)((refaToT[j]-minToT)/deltaToT)+1;
614     if ((refaToT[j] != 0) && (refaTime[j] != 0)){
615       vTime[nx]+=refaTime[j];
616       vTime2[nx]+=(refaTime[j])*(refaTime[j]);
617       vToT[nx]+=refaToT[j];
618       vToT2[nx]+=refaToT[j]*refaToT[j];
619       nentrx[nx]++;
620       hSlewing->Fill(refaToT[j],refaTime[j]);
621     }
622   }
623   Int_t nbinsToT=hSlewing->GetNbinsX();
624   if (nbinsToT != knbinToT) {
625     AliError("Profile :: incompatible numbers of bins");
626     return 0x0;
627   }
628
629   Int_t usefulBins=0;
630   TH1F *histo = new TH1F("histo", "1D Time vs ToT", nbinsToT, minToT, maxToT);
631   for (Int_t i=1;i<=nbinsToT;i++){
632     if (nentrx[i]!=0){
633     n[usefulBins]+=nentrx[i];
634     if (n[usefulBins]==0 && i == nbinsToT) {
635       break;
636     }
637     meanTime[usefulBins]+=vTime[i];
638     meanTime2[usefulBins]+=vTime2[i];
639     meanToT[usefulBins]+=vToT[i];
640     meanToT2[usefulBins]+=vToT2[i];
641     if (n[usefulBins]<20 && i!=nbinsToT) continue; 
642     mTime[usefulBins]=meanTime[usefulBins]/n[usefulBins];
643     mToT[usefulBins]=meanToT[usefulBins]/n[usefulBins];
644     sigmaTime[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
645                                    *(meanTime2[usefulBins]-meanTime[usefulBins]
646                                      *meanTime[usefulBins]/n[usefulBins]));
647     if ((1./n[usefulBins]/n[usefulBins]
648          *(meanToT2[usefulBins]-meanToT[usefulBins]
649            *meanToT[usefulBins]/n[usefulBins]))< 0) {
650       AliError(" too small radical" );
651       sigmaToT[usefulBins]=0;
652     }
653     else{       
654       sigmaToT[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
655                                      *(meanToT2[usefulBins]-meanToT[usefulBins]
656                                        *meanToT[usefulBins]/n[usefulBins]));
657     }
658     usefulBins++;
659     }
660   }
661   for (Int_t i=0;i<usefulBins;i++){
662     Int_t binN = (Int_t)((mToT[i]-minToT)/deltaToT)+1;
663     histo->Fill(mToT[i],mTime[i]);
664     histo->SetBinError(binN,sigmaTime[i]);
665   } 
666   return histo;
667 }
668 //_____________________________________________________________________________
669
670 void AliTOFcalib::CorrectESDTime()
671 {
672   //Calculate the corrected TOF time
673   Int_t nelements = fESDsel->GetEntries();
674   for (Int_t i=0; i< nelements; i++) {
675     AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(i);
676     Int_t index = element->GetTOFCalChannel();
677     Float_t tToT = element->GetToT();
678     //select the correspondent channel with its simulated ToT spectrum
679     //summing up everything, index = 0 for all channels:
680     Int_t ipid = element->GetCombID();
681     Double_t etime = 0;   //expected time
682     Double_t expTime[10]; 
683     element->GetIntegratedTimes(expTime);
684     if (ipid == 0) etime = expTime[2]*1E-3; //ns
685     else if (ipid == 1) etime = expTime[3]*1E-3; //ns
686     else if (ipid == 2) etime = expTime[4]*1E-3; //ns
687     Float_t par[6];
688     if (!fTOFCal) {
689       AliTOFGeometry *geom=new AliTOFGeometryV5();
690       fTOFCal = new AliTOFCal(geom);
691       fTOFCal->CreateArray();
692       delete geom;
693     }
694     AliTOFChannel * calChannel = fTOFCal->GetChannel(index);
695     for (Int_t j = 0; j<6; j++){
696       par[j]=calChannel->GetSlewPar(j);
697     }
698     Float_t timeCorr=0;
699     timeCorr= par[0]+par[1]*tToT+par[2]*tToT*tToT+par[3]*tToT*tToT*tToT+par[4]*tToT*tToT*tToT*tToT+par[5]*tToT*tToT*tToT*tToT*tToT;
700   }
701 }
702 //_____________________________________________________________________________
703
704 void AliTOFcalib::CorrectESDTime(AliESD *event){
705   //Calculate the corrected TOF time
706
707   Int_t ntrk =0;
708   ntrk=event->GetNumberOfTracks();
709   for (Int_t itrk=0; itrk<ntrk; itrk++) {
710     AliESDtrack *t=event->GetTrack(itrk);
711     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
712       continue;
713     }
714     //IsStartedTimeIntegral
715     if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
716       continue;
717     }
718     UInt_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
719     if(assignedTOFcluster==0){ // not matched
720       continue;
721     }
722     Int_t index = t->GetTOFCalChannel();
723     if (!fTOFCal) {
724       AliTOFGeometry *geom=new AliTOFGeometryV5();
725       fTOFCal = new AliTOFCal(geom);
726       fTOFCal->CreateArray();
727       delete geom;
728     }
729     AliTOFChannel * calChannel = fTOFCal->GetChannel(index);
730     Float_t par[6];
731     for (Int_t j = 0; j<6; j++){
732       par[j]=calChannel->GetSlewPar(j);
733     }
734     Float_t tToT = t->GetTOFsignalToT();
735     Float_t timeCorr =0; 
736     timeCorr=par[0]+par[1]*tToT+par[2]*tToT*tToT+par[3]*tToT*tToT*tToT+par[4]*tToT*tToT*tToT*tToT+par[5]*tToT*tToT*tToT*tToT*tToT;
737   }
738 }
739 //_____________________________________________________________________________
740
741 void AliTOFcalib::WriteParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun)
742 {
743   //Write calibration parameters to the CDB
744   AliCDBManager *man = AliCDBManager::Instance();
745   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
746   Char_t *sel1 = "Par" ;
747   Char_t  out[100];
748   sprintf(out,"%s/%s",sel,sel1); 
749   AliCDBId id(out,minrun,maxrun);
750   AliCDBMetaData *md = new AliCDBMetaData();
751   md->SetResponsible("Chiara Zampolli");
752   if (!fTOFCal) {
753     AliTOFGeometry *geom=new AliTOFGeometryV5();
754     fTOFCal = new AliTOFCal(geom);
755     fTOFCal->CreateArray();
756     delete geom;
757   }
758   man->Put(fTOFCal,id,md);
759   delete md;
760 }
761 //_____________________________________________________________________________
762
763 void AliTOFcalib::WriteParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun, AliTOFCal *cal){
764   //Write calibration parameters to the CDB
765   AliCDBManager *man = AliCDBManager::Instance();
766   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
767   Char_t *sel1 = "Par" ;
768   Char_t  out[100];
769   sprintf(out,"%s/%s",sel,sel1); 
770   AliCDBId id(out,minrun,maxrun);
771   AliCDBMetaData *md = new AliCDBMetaData();
772   md->SetResponsible("Chiara Zampolli");
773   man->Put(cal,id,md);
774   delete md;
775 }
776 //_____________________________________________________________________________
777
778 void AliTOFcalib::ReadParFromCDB(Char_t *sel, Int_t nrun)
779 {
780   //Read calibration parameters from the CDB
781   AliCDBManager *man = AliCDBManager::Instance();
782   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
783   Char_t *sel1 = "Par" ;
784   Char_t  out[100];
785   sprintf(out,"%s/%s",sel,sel1); 
786   AliCDBEntry *entry = man->Get(out,nrun);
787   AliTOFCal *cal =(AliTOFCal*)entry->GetObject();
788   fTOFCal = cal;
789 }
790 //_____________________________________________________________________________
791 void AliTOFcalib::WriteSimParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun)
792 {
793   //Write Sim miscalibration parameters to the CDB
794
795
796   //for the time being, only one spectrum is used
797   TFile *spFile = new TFile("$ALICE_ROOT/TOF/data/spectrum.root","read");
798   TH1F * hToT;
799   // Retrieve ToT Spectrum
800   spFile->GetObject("ToT",hToT);
801
802   fTOFSimToT=hToT;
803   
804   // Retrieve Time over TOT dependence 
805   
806   TH1F * h = (TH1F*)spFile->Get("TimeToTFit");
807   TList * list = (TList*)h->GetListOfFunctions();
808   TF1* f = (TF1*)list->At(0);
809   Float_t par[6] = {0,0,0,0,0,0};
810   Int_t npar=f->GetNpar();
811   for (Int_t ipar=0;ipar<npar;ipar++){
812     par[ipar]=f->GetParameter(ipar);
813   }
814   if (!fTOFSimCal) {
815     AliTOFGeometry *geom=new AliTOFGeometryV5();
816     fTOFSimCal = new AliTOFCal(geom);
817     fTOFSimCal->CreateArray();
818     delete geom;
819   }
820   for(Int_t iTOFch=0; iTOFch<fTOFSimCal->NPads();iTOFch++){
821     AliTOFChannel * calChannel = fTOFSimCal->GetChannel(iTOFch);
822     calChannel->SetSlewPar(par);
823   }
824
825   // Store them in the CDB
826
827   AliCDBManager *man = AliCDBManager::Instance();
828   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
829   AliCDBMetaData *md = new AliCDBMetaData();
830   md->SetResponsible("Chiara Zampolli");
831   Char_t *sel1 = "SimPar" ;
832   Char_t  out[100];
833   sprintf(out,"%s/%s",sel,sel1); 
834   AliCDBId id1(out,minrun,maxrun);
835   man->Put(fTOFSimCal,id1,md);
836   Char_t *sel2 = "SimHisto" ;
837   sprintf(out,"%s/%s",sel,sel2); 
838   AliCDBId id2(out,minrun,maxrun);
839   man->Put(fTOFSimToT,id2,md);
840   delete md;
841 }
842
843 //_____________________________________________________________________________
844 void AliTOFcalib::WriteSimParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun, AliTOFCal *cal, TH1F * histo){
845   //Write Sim miscalibration parameters to the CDB
846
847   fTOFSimToT=histo;
848   fTOFSimCal=cal;  
849   AliCDBManager *man = AliCDBManager::Instance();
850   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
851   AliCDBMetaData *md = new AliCDBMetaData();
852   md->SetResponsible("Chiara Zampolli");
853   Char_t *sel1 = "SimPar" ;
854   Char_t  out[100];
855   sprintf(out,"%s/%s",sel,sel1); 
856   AliCDBId id1(out,minrun,maxrun);
857   man->Put(fTOFSimCal,id1,md);
858   Char_t *sel2 = "SimHisto" ;
859   sprintf(out,"%s/%s",sel,sel2); 
860   AliCDBId id2(out,minrun,maxrun);
861   man->Put(fTOFSimToT,id2,md);
862   delete md;
863 }
864 //_____________________________________________________________________________
865 void AliTOFcalib::ReadSimParFromCDB(Char_t *sel, Int_t nrun)
866 {
867   //Read miscalibration parameters from the CDB
868   AliCDBManager *man = AliCDBManager::Instance();
869   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
870   Char_t *sel1 = "SimPar" ;
871   Char_t  out[100];
872   sprintf(out,"%s/%s",sel,sel1); 
873   AliCDBEntry *entry1 = man->Get(out,nrun);
874   AliTOFCal *cal =(AliTOFCal*)entry1->GetObject();
875   fTOFSimCal=cal;
876   Char_t *sel2 = "SimHisto" ;
877   sprintf(out,"%s/%s",sel,sel2); 
878   AliCDBEntry *entry2 = man->Get(out,nrun);
879   TH1F *histo =(TH1F*)entry2->GetObject();
880   fTOFSimToT=histo;
881 }
882 //_____________________________________________________________________________
883
884 Int_t AliTOFcalib::GetIndex(Int_t *detId)
885 {
886   //Retrieve calibration channel index 
887   Int_t isector = detId[0];
888   if (isector >= fNSector)
889     AliError(Form("Wrong sector number in TOF (%d) !",isector));
890   Int_t iplate = detId[1];
891   if (iplate >= fNPlate)
892     AliError(Form("Wrong plate number in TOF (%d) !",iplate));
893   Int_t istrip = detId[2];
894   Int_t ipadz = detId[3];
895   Int_t ipadx = detId[4];
896   Int_t stripOffset = 0;
897   switch (iplate) {
898   case 0:
899     stripOffset = 0;
900     break;
901   case 1:
902     stripOffset = fNStripC;
903     break;
904   case 2:
905     stripOffset = fNStripC+fNStripB;
906     break;
907   case 3:
908     stripOffset = fNStripC+fNStripB+fNStripA;
909     break;
910   case 4:
911     stripOffset = fNStripC+fNStripB+fNStripA+fNStripB;
912     break;
913   default:
914     AliError(Form("Wrong plate number in TOF (%d) !",iplate));
915     break;
916   };
917
918   Int_t idet = ((2*(fNStripC+fNStripB)+fNStripA)*fNpadZ*fNpadX)*isector +
919                (stripOffset*fNpadZ*fNpadX)+
920                (fNpadZ*fNpadX)*istrip+
921                (fNpadX)*ipadz+
922                 ipadx;
923   return idet;
924 }
925