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