Task for TOF Calibration, optimized for pp runs
[u/mrichter/AliRoot.git] / TOF / AliTOFCalibTask.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 $Log
17 */
18
19 // task to perform TOF calibration
20 // optimized for pp runs
21 // expect a maximum of 100 entries per channel
22 // different ways to calibrate are implemented
23 // C. Zampolli 
24
25 #include <TChain.h>
26 #include <TObject.h> 
27 #include <TCanvas.h>
28 #include <TROOT.h>
29 #include <TSystem.h>
30 #include <TProfile.h>
31 #include <TF1.h>
32 #include <TTree.h>
33 #include <TFile.h>
34 #include <TH1F.h>
35 #include <TH1I.h>
36 #include <TH1D.h>
37 #include <TH2F.h>
38   
39 #include "AliTOFCalibTask.h" 
40 #include "AliTOFChannelTask.h" 
41 #include "AliESD.h" 
42 #include "AliESDtrack.h" 
43 #include "AliLog.h"
44 //______________________________________________________________________________
45 AliTOFCalibTask::AliTOFCalibTask(const char *name) : 
46   AliAnalysisTask(name,""),  
47   fdir(0),
48   fChain(0),
49   fESD(0),
50   fMinEntries(0),
51   fIch(0),
52   fToT(0),
53   fTime(0),
54   fExpTimePi(0),
55   fExpTimeKa(0),
56   fExpTimePr(0),
57   ftree(0x0),
58   fMinTime(0),
59   fbigarray(0x0),
60   findexarray(0x0),
61   fnESD(0),
62   fnESDselected(0),
63   fnESDkTOFout(0),
64   fnESDkTIME(0),
65   fnESDassTOFcl(0),
66   fnESDTIMEcut(0),
67   fnESDTRDcut(0),
68   fhToT(0),
69   fhTime(0),
70   fhExpTimePi(0),
71   fhExpTimeKa(0),
72   fhExpTimePr(0),
73   fhPID(0),
74   fhch(0),
75   fOutputContainer(0)
76 {
77   // Constructor.
78   // Input slot #0 works with an Ntuple
79
80   DefineInput(0, TChain::Class());
81   // Output slot #0 writes into a TH1 container
82   DefineOutput(0,  TObjArray::Class()) ; 
83   fdir = gDirectory->GetPath();
84 }
85
86 //______________________________________________________________________________
87 AliTOFCalibTask::AliTOFCalibTask(const AliTOFCalibTask &calibtask) : 
88   AliAnalysisTask("AliTOFCalibTask",""),  
89   fdir(0),
90   fChain(0),
91   fESD(0),
92   fMinEntries(0),
93   fIch(0),
94   fToT(0),
95   fTime(0),
96   fExpTimePi(0),
97   fExpTimeKa(0),
98   fExpTimePr(0),
99   ftree(0x0),
100   fMinTime(0),
101   fbigarray(0x0),
102   findexarray(0x0),
103   fnESD(0),
104   fnESDselected(0),
105   fnESDkTOFout(0),
106   fnESDkTIME(0),
107   fnESDassTOFcl(0),
108   fnESDTIMEcut(0),
109   fnESDTRDcut(0),
110   fhToT(0),
111   fhTime(0),
112   fhExpTimePi(0),
113   fhExpTimeKa(0),
114   fhExpTimePr(0),
115   fhPID(0),
116   fhch(0),
117   fOutputContainer(0)
118 {
119   // Copy Constructor.
120   fdir=calibtask.fdir;
121   fChain=calibtask.fChain;
122   fESD=calibtask.fESD;
123   fMinEntries=calibtask.fMinEntries;
124   fIch=calibtask.fIch;
125   fToT=calibtask.fToT;
126   fTime=calibtask.fTime;
127   fExpTimePi=calibtask.fExpTimePi;
128   fExpTimeKa=calibtask.fExpTimeKa;
129   fExpTimePr=calibtask.fExpTimePr;
130   ftree=calibtask.ftree;
131   fMinTime=calibtask.fMinTime;
132   fbigarray=calibtask.fbigarray;
133   findexarray=calibtask.findexarray;
134   fnESD=calibtask.fnESD;
135   fnESDselected=calibtask.fnESDselected;
136   fnESDkTOFout=calibtask.fnESDkTOFout;
137   fnESDkTIME=calibtask.fnESDkTIME;
138   fnESDassTOFcl=calibtask.fnESDassTOFcl;
139   fnESDTIMEcut=calibtask.fnESDTIMEcut;
140   fnESDTRDcut=calibtask.fnESDTRDcut;
141   fhToT=calibtask.fhToT;
142   fhTime=calibtask.fhTime;
143   fhExpTimePi=calibtask.fhExpTimePi;
144   fhExpTimeKa=calibtask.fhExpTimeKa;
145   fhExpTimePr=calibtask.fhExpTimePr;
146   fhPID=calibtask.fhPID;
147   fhch=calibtask.fhch;
148   fOutputContainer=calibtask.fOutputContainer; 
149 }
150 //______________________________________________________________________________
151 AliTOFCalibTask:: ~AliTOFCalibTask() 
152 {
153   // destructor
154
155   AliInfo("TOF Calib Task: Deleting");
156   delete[] fbigarray;
157   delete findexarray;
158   delete fOutputContainer;
159   delete fhToT;
160   delete fhTime;
161   delete fhExpTimePi; 
162   delete fhExpTimeKa; 
163   delete fhExpTimePr; 
164   delete fhPID; 
165   delete fhch;
166 }
167 //______________________________________________________________________________
168 AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)  
169
170    //assignment operator
171   this->fdir=calibtask.fdir;
172   this->fChain=calibtask.fChain;
173   this->fESD=calibtask.fESD;
174   this->fMinEntries=calibtask.fMinEntries;
175   this->fIch=calibtask.fIch;
176   this->fToT=calibtask.fToT;
177   this->fTime=calibtask.fTime;
178   this->fExpTimePi=calibtask.fExpTimePi;
179   this->fExpTimeKa=calibtask.fExpTimeKa;
180   this->fExpTimePr=calibtask.fExpTimePr;
181   this->ftree=calibtask.ftree;
182   this->fMinTime=calibtask.fMinTime;
183   this->fbigarray=calibtask.fbigarray;
184   this->findexarray=calibtask.findexarray;
185   this->fnESD=calibtask.fnESD;
186   this->fnESDselected=calibtask.fnESDselected;
187   this->fnESDkTOFout=calibtask.fnESDkTOFout;
188   this->fnESDkTIME=calibtask.fnESDkTIME;
189   this->fnESDassTOFcl=calibtask.fnESDassTOFcl;
190   this->fnESDTIMEcut=calibtask.fnESDTIMEcut;
191   this->fnESDTRDcut=calibtask.fnESDTRDcut;
192   this->fhToT=calibtask.fhToT;
193   this->fhTime=calibtask.fhTime;
194   this->fhExpTimePi=calibtask.fhExpTimePi;
195   this->fhExpTimeKa=calibtask.fhExpTimeKa;
196   this->fhExpTimePr=calibtask.fhExpTimePr;
197   this->fOutputContainer=calibtask.fOutputContainer; 
198   this->fhPID=calibtask.fhPID;
199   this->fhch=calibtask.fhch;
200   return *this;
201 }
202 //--------------------------------------------------------------------------
203 void AliTOFCalibTask::BookHistos(){
204
205   //booking histos
206
207   AliInfo(Form("*** Booking Histograms %s", GetName())) ; 
208
209   fhToT=
210     new TH1F("hToT", " ToT distribution (ns) ", 400, 0, 40);
211   fhTime=
212     new TH1F("hTime", " Time distribution (ns) ", 400, 0, 40);
213   fhExpTimePi=
214     new TH1F("hExpTimePi", " Exp Time distribution, Pi (ns) ", 400, 0, 40);
215   fhExpTimeKa=
216     new TH1F("hExpTimeKa", " Exp Time distribution, Ka (ns) ", 400, 0, 40);
217   fhExpTimePr=
218     new TH1F("hExpTimePr", " Exp Time distribution, Pr (ns) ", 400, 0, 40);
219   fhPID=
220     new TH1I("hPID", " Combinatorial PID identity ", 3, 0, 3);
221   fhch=
222     new TH1D("hch", " TOF channel ", TOFCHANNELS, 0, TOFCHANNELS);
223
224   //create the putput container
225   fOutputContainer = new TObjArray(7) ; 
226   fOutputContainer->SetName(GetName()) ; 
227
228   fOutputContainer->AddAt(fhToT,             0) ; 
229   fOutputContainer->AddAt(fhTime,            1) ; 
230   fOutputContainer->AddAt(fhExpTimePi,       2) ; 
231   fOutputContainer->AddAt(fhExpTimeKa,       3) ; 
232   fOutputContainer->AddAt(fhExpTimePr,       4) ; 
233   fOutputContainer->AddAt(fhPID,             5) ; 
234   fOutputContainer->AddAt(fhch,              6) ; 
235
236 }
237 //----------------------------------------------------------------------------
238 void AliTOFCalibTask::DrawHistos(){
239
240   // drawing output histos
241
242   AliInfo(Form("*** Drawing Histograms %s", GetName())) ; 
243
244   TCanvas * canvasToTTime = new TCanvas("canvasToTTime", " ToT and Time ",400, 30, 550, 630) ;
245   canvasToTTime->Divide(1,2);
246   canvasToTTime->cd(1);
247   fhToT->SetLineColor(4);
248   fhToT->GetXaxis()->SetTitle("ToT (ns)");
249   fhToT->Draw("hist");
250   canvasToTTime->cd(2);
251   fhTime->SetLineColor(4);
252   fhTime->GetXaxis()->SetTitle("Time (ns)");
253   fhTime->Draw("hist");
254   canvasToTTime->Update();
255   canvasToTTime->Print("ToTTime.gif");
256
257   TCanvas * canvasExpTime = new TCanvas("canvasExpTime", " Expected Times ",400, 30, 550, 630) ;
258   canvasExpTime->Divide(1,3);
259   canvasExpTime->cd(1);
260   fhExpTimePi->SetLineColor(4);
261   fhExpTimePi->GetXaxis()->SetTitle("Exp Time (ns), #pi");
262   fhExpTimePi->Draw("hist");
263   canvasExpTime->cd(2);
264   fhExpTimeKa->SetLineColor(4);
265   fhExpTimeKa->GetXaxis()->SetTitle("Exp Time (ns), K");
266   fhExpTimeKa->Draw("hist");
267   canvasExpTime->cd(3);
268   fhExpTimePr->SetLineColor(4);
269   fhExpTimePr->GetXaxis()->SetTitle("Exp Time (ns), p");
270   fhExpTimePr->Draw("hist");
271
272   canvasExpTime->Print("ExpTime.gif");
273
274   TCanvas * canvasPID = new TCanvas("canvasPID", " Combinatorial PID ",400, 30, 550, 400);
275   fhPID->GetXaxis()->SetTitle("Comb PID");
276   fhPID->GetXaxis()->SetBinLabel(1,"#pi");
277   fhPID->GetXaxis()->SetBinLabel(2,"K");
278   fhPID->GetXaxis()->SetBinLabel(3,"p");
279   fhPID->Draw("hist");
280
281   canvasPID->Print("PID.gif");
282
283   TCanvas * canvasrndch = new TCanvas("canvasrndch", " TOF channel ",400, 30, 550, 400);
284   fhch->GetXaxis()->SetTitle("TOF ch");
285   fhch->Draw("hist");
286   Float_t meanTOFch = 0;
287   for (Int_t ibin=0;ibin<TOFCHANNELS;ibin++){
288     meanTOFch+=(Float_t)fhch->GetBinContent(ibin+1);
289   }
290
291   meanTOFch/=TOFCHANNELS;
292   AliDebug(1,Form(" Mean number of tracks/channel = %f ",meanTOFch));
293
294   canvasrndch->Print("rndch.gif");
295
296   char line[1024] ; 
297   sprintf(line, ".!tar -zcvf %s.tar.gz *.gif", GetName()) ; 
298   gROOT->ProcessLine(line);
299   sprintf(line, ".!rm -fR *.gif"); 
300   gROOT->ProcessLine(line);
301   AliInfo(Form("*** TOF Calib Task: plots saved in %s.tar.gz...\n", GetName())) ;
302 }
303
304 //______________________________________________________________________________
305 void AliTOFCalibTask::ConnectInputData(const Option_t*)
306 {
307   // Initialization of branch container and histograms 
308     
309   //  AliLog::SetClassDebugLevel("AliTOFCalibTask",1);
310   AliInfo(Form("*** Initialization of %s", GetName())) ; 
311   
312   // Get input data
313   fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
314   if (!fChain) {
315     AliError(Form("Input 0 for %s not found\n", GetName()));
316     return ;
317   }
318   
319   // One should first check if the branch address was taken by some other task
320   char ** address = (char **)GetBranchAddress(0, "ESD");
321   if (address) {
322     fESD = (AliESD*)(*address);
323   } else {
324     fESD = new AliESD();
325     fESD->ReadFromTree(fChain) ;  
326   }
327   fbigarray = new Float_t*[TOFCHANNELS];
328   findexarray = new Int_t[TOFCHANNELS];
329   for (Int_t i =0;i<TOFCHANNELS;i++){
330     fbigarray[i] = new Float_t[CHENTRIES];
331     findexarray[i]=0;
332     for (Int_t j =0;j<CHENTRIES;j++){
333       fbigarray[i][j]=-1;
334     }
335   }
336
337   BookHistos();
338
339 }
340 //________________________________________________________________________
341 void AliTOFCalibTask::CreateOutputObjects()
342 {
343 // Create histograms
344 }
345
346 //______________________________________________________________________________
347 void AliTOFCalibTask::Exec(Option_t * opt) 
348 {
349
350   // main 
351
352   AliInfo(Form("*** Executing %s", GetName())) ; 
353
354 //******* The loop over events -----------------------------------------------
355
356   // Processing of one event
357   Long64_t entry = fChain->GetReadEntry() ;  
358   if (!fESD) {
359     AliError("fESD is not connected to the input!") ; 
360     return ; 
361   }
362
363   if ( !((entry-1)%100) ) 
364     AliDebug(1,Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ; 
365   
366   // ************************  TOF *************************************
367
368   fMinTime=22E3;   //ns; not used
369   Int_t ntrk = fESD->GetNumberOfTracks() ;
370   fnESD+=ntrk;
371   Int_t nselected = 0;
372   Int_t itr = -1;
373   while ( ntrk-- ) {
374     itr++;
375     AliESDtrack * t = fESD->GetTrack(ntrk) ;
376     //selecting only good quality tracks
377     if (!Select(t)) continue;
378     nselected++;
379     Int_t ich = Int_t(t->GetTOFCalChannel()); 
380     fhch->Fill(ich);
381     AliDebug(2,Form(" ESD in channel %i, filling array %i",t->GetTOFCalChannel(),ich));
382     Float_t tot = t->GetTOFsignalToT();
383     Float_t time = t->GetTOFsignalRaw();
384     AliDebug(2,Form(" track # %i in channel %i, time = %f \n",ntrk,ich,time)); 
385     Double_t expTime[10]; 
386     t->GetIntegratedTimes(expTime);
387     Float_t expTimePi = expTime[2]*1.E-3;
388     Float_t expTimeKa = expTime[3]*1.E-3;
389     Float_t expTimePr = expTime[4]*1.E-3;
390     if (findexarray[ich]==(Int_t)(CHENTRIES/NIDX)) {
391       AliInfo(Form("too many tracks in channel %i, not storing any more...",ich));
392       continue;
393     } 
394     findexarray[ich]++;
395     AliDebug(2,Form("tracks in channel %i = %i, storing... ",ich, findexarray[ich] ));
396     Int_t ientry=(findexarray[ich]-1)*NIDX;
397     fbigarray[ich][ientry+DELTAIDXTOT]=tot;  //in ns
398     fbigarray[ich][ientry+DELTAIDXTIME]=time*1E-3; // in ns
399     fbigarray[ich][ientry+DELTAIDXEXTIMEPI]=expTimePi;
400     fbigarray[ich][ientry+DELTAIDXEXTIMEKA]=expTimeKa;
401     fbigarray[ich][ientry+DELTAIDXEXTIMEPR]=expTimePr;
402     fhToT->Fill(fbigarray[ich][ientry+DELTAIDXTOT]);
403     fhTime->Fill(fbigarray[ich][ientry+DELTAIDXTIME]);
404     fhExpTimePi->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEPI]);
405     fhExpTimeKa->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEKA]);
406     fhExpTimePr->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEPR]);
407     AliDebug(2, Form("track = %i, tot = %f, time = %f, and Exp time in TOF: pi = %f, K = %f, p = %f",itr, fbigarray[ich][ientry+DELTAIDXTOT], fbigarray[ich][ientry+DELTAIDXTIME], expTimePi,expTimeKa,expTimePr));
408   }
409   fnESDselected+=nselected;
410
411   PostData(0, fOutputContainer);  
412 }
413
414 //_____________________________________________________________________________
415 void AliTOFCalibTask::Terminate(Option_t *)
416 {
417   // Processing when the event loop is ended
418   
419   // some plots
420
421   TH1::AddDirectory(0);
422   AliInfo("TOF Calib Task: End of events loop");
423   AliInfo(Form(" Number of analyzed ESD tracks: %i\n",fnESD));
424   AliInfo(Form(" Number of selected ESD tracks: %i\n",fnESDselected));
425   AliInfo(Form(" Number of ESD tracks with kTOFout: %i\n",fnESDkTOFout));
426   AliInfo(Form(" Number of ESD tracks with kTIME: %i\n",fnESDkTIME));
427   AliInfo(Form(" Number of ESD tracks with TRDcut: %i\n",fnESDTRDcut));
428   AliInfo(Form(" Number of ESD tracks with TIMEcut: %i\n",fnESDTIMEcut));
429   AliInfo(Form(" Number of ESD tracks with assTOFcl: %i\n",fnESDassTOFcl));
430
431   for (Int_t i = 0;i<TOFCHANNELS;i++){
432     Int_t size=findexarray[i]*NIDX;
433     AliDebug(2, Form(" entries %i in channel %i ",findexarray[i],i));
434     if (findexarray[i]<6) {
435       AliDebug(1, Form(" not enough statistics for combined PID for channel %i, putting all the tracks as if they were pions",i));
436       continue;
437     }
438     if (!CombPID(&fbigarray[i][0], size)) AliError("ERROR!!!!ERROR!!!");
439   }
440   DrawHistos();
441
442   // saving data in a tree
443   
444   TDirectory *dir = gDirectory;
445   TFile *outFile = 0;
446   TTree * tree = 0x0;
447   Bool_t isThere=kFALSE;
448   const char *dirname = "./";
449   TString filename= "TOFCalib.root";
450
451   if((gSystem->FindFile(dirname,filename))!=NULL){
452     isThere=kTRUE;
453   }
454   
455   if(!isThere){
456     AliInfo("File with tree for Calibration not there, creating it");
457     outFile = new TFile( "TOFCalib.root","RECREATE");
458     outFile->cd();
459     tree = new TTree("T", "Tree for TOF calibration");
460     Float_t p[CHENTRIESSMALL];
461     Int_t nentries;
462     tree->Branch("nentries",&nentries,"nentries/I");
463     tree->Branch("TOFentries",p,"TOFentries[nentries]/F");
464     for (Int_t i=0;i<TOFCHANNELS;i++){
465       nentries=findexarray[i]*(NIDXSMALL); // when filling small array, 
466                                            // only first 3 floats taken 
467                                            // into account
468       for (Int_t j=0; j<findexarray[i];j++){
469         for (Int_t k=0; k<NIDXSMALL;k++){
470           Int_t index1= j*NIDXSMALL+k;   // index in small array
471           Int_t index2=j*NIDX+k;         // index in big array
472           p[index1]=fbigarray[i][index2];
473         }
474       }
475       tree->Fill();
476     }
477     tree->Write("",TObject::kOverwrite);
478     outFile->Close();
479     //delete outFile;
480     //outFile=0x0;
481     //delete tree;
482     //tree=0x0;
483   }
484   else{
485     AliInfo("File with tree for Calibration already there, updating it");
486     outFile = new TFile( "TOFCalib.root","READ");
487     outFile->cd();
488     tree=(TTree*)outFile->Get("T");
489     Float_t p[CHENTRIESSMALL];
490     Float_t ptemp[MAXCHENTRIESSMALL];
491     Int_t nentries;
492     Int_t nentries1;
493     Int_t nentriestemp;
494     tree->SetBranchAddress("nentries",&nentries);
495     tree->SetBranchAddress("TOFentries",p);
496     TFile * outFile1 = new TFile( "TOFCalib.root","RECREATE");
497     outFile1->cd();
498     TTree *treeTemp = new TTree("T","Tree for TOF calibration");
499     treeTemp->Branch("nentries",&nentriestemp,"nentries/I");
500     treeTemp->Branch("TOFentries",ptemp,"TOFentries[nentries]/F");
501     for (Int_t i=0;i<TOFCHANNELS;i++){
502       tree->GetEntry(i);
503       nentries1=nentries;
504       if (nentries+findexarray[i]*(NIDXSMALL) > MAXCHENTRIESSMALL){
505         AliDebug(2, Form("findexarray[%i] = %i, nentries = %i",i,findexarray[i],nentries)); 
506         AliInfo(Form("Too many entries in channel %i, stopping at %i  entries, no storing any more...",i,MAXCHENTRIESSMALL/NIDXSMALL));
507         findexarray[i]=(MAXCHENTRIESSMALL-nentries)/NIDXSMALL;
508       }
509       for (Int_t kk = 0;kk<nentries;kk++){
510         ptemp[kk]=p[kk];
511       }
512       nentries+=findexarray[i]*(NIDXSMALL); // when filling small array, 
513                                             // only first 3 floats taken 
514                                             // into account 
515       nentriestemp=nentries; 
516       for (Int_t j=0; j<findexarray[i];j++){
517         for (Int_t k=0; k<NIDXSMALL;k++){
518           Int_t index1= j*NIDXSMALL+k;   // index in small array
519           Int_t index2=j*NIDX+k;   // index in big array
520           ptemp[index1+nentries1]=fbigarray[i][index2];
521         }
522       }
523       treeTemp->Fill();
524     }
525     treeTemp->SetName("T");
526     treeTemp->SetTitle("Tree for TOF calibration");
527     outFile->Close();
528     treeTemp->Write("",TObject::kOverwrite);
529     outFile1->Close();
530   }
531   dir->cd();
532   Int_t calibrationStatus = 0;
533   Int_t ch[2]={3,1003};
534   calibrationStatus = Calibrate(2,ch,"save");
535   //calibrationStatus = Calibrate(3,4,"save");
536   //calibrationStatus = CalibrateFromProfile(3,"save");
537   //calibrationStatus = Calibrate(3,"save");
538   //  calibrationStatus = Calibrate("");
539   AliInfo(Form("Calibration Status = %i, bye.... \n",calibrationStatus));
540   
541 }
542 //----------------------------------------------------------------------------
543 Int_t AliTOFCalibTask::Calibrate(Int_t ichmin, Int_t ichmax, Option_t *optionSave, Option_t *optionFit){
544
545   // calibrating summing more than one channels
546   // computing calibration parameters
547   // Returning codes:
548   // 0 -> everything was ok
549   // 1 -> no tree for calibration found
550   // 2 -> not enough statistics to perform calibration
551   // 3 -> problems with arrays
552   
553   TH1::AddDirectory(0);
554
555   AliInfo(Form("*** Calibrating Histograms %s, summing more channels, from channel %i, to channel %i, storing Calib Pars in channel %i", GetName(),ichmin,ichmax,ichmin)) ; 
556   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
557   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
558   const char *dirname = "./";
559   TString filename= "TOFCalib.root";
560   
561   if((gSystem->FindFile(dirname,filename))==NULL){
562     AliInfo("No file with tree for calibration found! exiting....");
563     return 1;
564   }
565
566   TFile *file = new TFile("TOFCalib.root","READ");
567   TTree *tree = (TTree*)file->Get("T");
568   Float_t p[MAXCHENTRIESSMALL];
569   Int_t nentries;
570   tree->SetBranchAddress("nentries",&nentries);
571   tree->SetBranchAddress("TOFentries",p);
572
573   Float_t nentriesmean =0;
574   for (Int_t i=ichmin; i<ichmax; i++){
575     tree->GetEntry(i);
576     Int_t ntracks=nentries/3;
577     nentriesmean+=ntracks;
578   }
579
580   //  nentriesmean/=(ichmax-ichmin);
581   if (nentriesmean < MEANENTRIES) {
582     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
583     return 2;
584   }
585
586   //filling ToT and Time arrays
587
588   Int_t nbinToT = 100;  // ToT bin width in Profile = 48.8 ps 
589   Float_t minToT = 0;   // ns
590   Float_t maxToT = 4.88;  // ns
591   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
592
593   TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
594   TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
595   Int_t ntracksTot = 0;
596   Int_t ntracks = 0;
597   Double_t binsProfile[101]; // sized larger than necessary, the correct 
598                              // dim being set in the booking of the profile
599   Int_t nusefulbins=0;
600   Float_t meantime=0;
601   for (Int_t i = ichmin;i<ichmax;i++){
602     tree->GetEntry(i);
603     ntracksTot+=nentries/3;
604     ntracks=nentries/3;
605     AliDebug(2,Form("channel %i, nentries = %i, ntracks = %i",i,nentries, ntracks));
606     for (Int_t j=0;j<ntracks;j++){
607       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
608       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
609       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
610       Float_t tot = p[idxexToT];
611       hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
612       meantime+=p[idxexTime]-p[idxexExTime];
613       hToT->Fill(tot);
614     }
615   }
616   nusefulbins = FindBins(hToT,&binsProfile[0]);
617   meantime/=ntracksTot;
618   AliDebug(2, Form("meantime = %f",meantime));
619   
620   for (Int_t j=1;j<=nusefulbins;j++) {
621     AliDebug(2,Form(" summing channels from %i to %i, nusefulbins = %i, bin %i = %f",ichmin,ichmax,nusefulbins,j,binsProfile[j])); 
622   }
623
624   TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G");  // CHECK THE BUILD OPTION, PLEASE!!!!!!
625   TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
626   for (Int_t i=ichmin; i<ichmax; i++){
627     tree->GetEntry(i);
628     ntracks=nentries/3;
629     for (Int_t j=0;j<ntracks;j++){
630       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
631       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
632       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
633       Float_t tot = p[idxexToT];
634       Float_t time = p[idxexTime]-p[idxexExTime];
635       AliDebug (2, Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
636       hSlewingProf->Fill(tot,time);  // if meantime is not used, the fill may be moved in the loop above
637       htimetot->Fill(tot,time-meantime);  // if meantime is not used, the fill may be moved in the loop above
638     }
639   }
640   hSlewingProf->Fit("pol5",optionFit,"",0,4);
641   TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
642   Float_t par[6];    
643   for(Int_t kk=0;kk<6;kk++){
644     par[kk]=calibfunc->GetParameter(kk);
645     AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
646   }
647
648   if(strstr(optionSave,"save")){
649     TFile * fileProf = new TFile("TOFCalibSave.root","recreate");
650     fileProf->cd(); 
651     TString profName=Form("Profile%06i_%06i",ichmin,ichmax);
652     TString timeTotName=Form("TimeTot%06i_%06i",ichmin,ichmax);
653     TString totName=Form("Tot%06i_%06i",ichmin,ichmax);
654     TString deltaName=Form("Delta%06i_%06i",ichmin,ichmax);
655     hSlewingProf->Write(profName);
656     htimetot->Write(timeTotName);
657     hToT->Write(totName);
658     hdeltaTime->Write(deltaName);
659     fileProf->Close();
660     delete fileProf;
661     fileProf=0x0;
662   }
663
664   delete hToT;
665   hToT=0x0;
666   delete hSlewingProf;
667   hSlewingProf=0x0;
668   delete htimetot;
669   htimetot=0x0;
670   delete hdeltaTime;
671   hdeltaTime=0x0;
672   file->Close();
673   delete file;
674   file=0x0;
675
676   AliTOFChannelTask * calChannel = new AliTOFChannelTask();
677   calChannel->SetSlewPar(par);
678   // saving parameters in chmin
679   arrayCal->AddAt(calChannel,ichmin);
680   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
681   arrayCal->Write("array",TObject::kSingleKey);
682   filecalib->Close();
683   delete filecalib;
684   filecalib = 0x0;
685   delete arrayCal;
686   arrayCal = 0x0;
687   delete calChannel;
688   calChannel =0x0;
689   return 0;
690 }
691 //----------------------------------------------------------------------------
692 Int_t AliTOFCalibTask::Calibrate(Int_t i, Option_t *optionSave, Option_t *optionFit){
693
694   // computing calibration parameters for channel i
695   // Returning codes:
696   // 0 -> everything was ok
697   // 1 -> no tree for calibration found
698   // 2 -> not enough statistics to perform calibration
699   // 3 -> problems with arrays
700
701   TH1::AddDirectory(0);
702   
703   AliInfo(Form("*** Calibrating Histograms (one channel) %s", GetName())) ; 
704   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
705   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
706   const char *dirname = "./";
707   TString filename= "TOFCalib.root";
708   
709   if((gSystem->FindFile(dirname,filename))==NULL){
710     AliInfo("No file with tree for calibration found! exiting....");
711     return 1;
712   }
713
714   TFile *file = new TFile("TOFCalib.root","READ");
715   TTree *tree = (TTree*)file->Get("T");
716   Float_t p[MAXCHENTRIESSMALL];
717   Int_t nentries;
718   tree->SetBranchAddress("nentries",&nentries);
719   tree->SetBranchAddress("TOFentries",p);
720
721   tree->GetEntry(i);
722   Int_t nentriesmean=nentries/3;
723
724   if (nentriesmean < MEANENTRIES) {  
725     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
726     return 2;
727   }
728
729   //filling ToT and Time arrays
730
731   Int_t nbinToT = 100;  // ToT bin width in Profile = 48.8 ps 
732   Float_t minToT = 0;   // ns
733   Float_t maxToT = 4.88;  // ns
734   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
735
736   TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
737   TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
738   Int_t ntracks = 0;
739   Double_t binsProfile[101]; // sized larger than necessary, the correct 
740                              // dim being set in the booking of the profile
741   Int_t nusefulbins=0;
742   Float_t meantime=0;
743   ntracks=nentries/3;
744   AliDebug(2,Form("channel %i, nentries = %i, ntracks = %i",i ,nentries, ntracks));
745   for (Int_t j=0;j<ntracks;j++){
746     Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
747     Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
748     Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
749     Float_t tot = p[idxexToT];
750     meantime+=p[idxexTime]-p[idxexExTime];
751     hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
752     hToT->Fill(tot);
753   }
754   
755   nusefulbins = FindBins(hToT,&binsProfile[0]);
756   meantime/=ntracks;
757   AliDebug(2,Form("meantime = %f",meantime));
758   
759   for (Int_t j=1;j<=nusefulbins;j++) {
760     AliDebug(2,Form(" channel %i, nusefulbins = %i, bin %i = %f",i,nusefulbins,j,binsProfile[j])); 
761   }
762
763   TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G");  // CHECK THE BUILD OPTION, PLEASE!!!!!!
764   TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
765   tree->GetEntry(i);
766   ntracks=nentries/3;
767   for (Int_t j=0;j<ntracks;j++){
768     Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
769     Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
770     Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
771     Float_t tot = p[idxexToT];
772     Float_t time = p[idxexTime]-p[idxexExTime];
773     AliDebug (2,Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
774     hSlewingProf->Fill(tot,time);  // if meantime is not used, the fill may be moved in the loop above
775     htimetot->Fill(tot,time-meantime);  // if meantime is not used, the fill may be moved in the loop above
776   }
777   hSlewingProf->Fit("pol5",optionFit,"",0,4);
778   TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
779   Float_t par[6];    
780   for(Int_t kk=0;kk<6;kk++){
781     par[kk]=calibfunc->GetParameter(kk);
782     AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
783   }
784
785   file->Close();
786   delete file;
787   file=0x0;
788
789   if(strstr(optionSave,"save")){
790     TFile * fileProf = new TFile("TOFCalibSave.root","recreate");
791     fileProf->cd();   
792     TString profName=Form("Profile%06i",i);
793     TString timeTotName=Form("TimeTot%06i",i);
794     TString totName=Form("Tot%06i",i);
795     TString deltaName=Form("Delta%06i",i);
796     hSlewingProf->Write(profName);
797     htimetot->Write(timeTotName);
798     hToT->Write(totName);
799     hdeltaTime->Write(deltaName);
800     fileProf->Close();
801     delete fileProf;
802     fileProf=0x0;
803   }
804
805   delete hToT;
806   hToT=0x0; 
807   delete hSlewingProf;
808   hSlewingProf=0x0;
809   delete htimetot;
810   htimetot=0x0;
811   delete hdeltaTime;
812   hdeltaTime=0x0;
813
814   AliTOFChannelTask * calChannel = new AliTOFChannelTask();
815   calChannel->SetSlewPar(par);
816   arrayCal->AddAt(calChannel,i);
817   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
818   arrayCal->Write("array",TObject::kSingleKey);
819   filecalib->Close();
820   delete filecalib;
821   filecalib=0x0;
822   delete arrayCal;
823   arrayCal = 0x0;
824   delete calChannel;
825   calChannel =0x0;
826   
827   return 0;
828 }
829 //----------------------------------------------------------------------------
830 Int_t AliTOFCalibTask::Calibrate(Int_t nch, Int_t *ch, Option_t *optionSave, Option_t *optionFit){
831
832   // calibrating an array of channels
833   // computing calibration parameters
834   // Returning codes:
835   // 0 -> everything was ok
836   // 1 -> no tree for calibration found
837   // 2 -> not enough statistics to perform calibration
838   // 3 -> problems with arrays
839   
840   TH1::AddDirectory(0);
841
842   AliInfo(Form("*** Calibrating Histograms %s, number of channels = %i", GetName(),nch)) ; 
843   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
844   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
845   for (Int_t ich=0; ich<nch; ich++){
846     Int_t i = ch[ich];
847     AliInfo(Form("Calibrating channel = %i",i )) ; 
848   }
849   const char *dirname = "./";
850   TString filename= "TOFCalib.root";
851   
852   if((gSystem->FindFile(dirname,filename))==NULL){
853     AliInfo("No file with tree for calibration found! exiting....");
854     return 1;
855   }
856
857   TFile *file = new TFile("TOFCalib.root","READ");
858   TTree *tree = (TTree*)file->Get("T");
859   Float_t p[MAXCHENTRIESSMALL];
860   Int_t nentries;
861   tree->SetBranchAddress("nentries",&nentries);
862   tree->SetBranchAddress("TOFentries",p);
863
864   Float_t nentriesmean =0;
865   for (Int_t ich=0; ich<nch; ich++){
866     Int_t i = ch[ich];
867     tree->GetEntry(i);
868     Int_t ntracks=nentries/3;
869     nentriesmean+=ntracks;
870   }
871
872   nentriesmean/=nch;
873   if (nentriesmean < MEANENTRIES) { 
874     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
875     return 2;
876   }
877
878   //filling ToT and Time arrays
879
880   Int_t nbinToT = 100;  // ToT bin width in Profile = 48.8 ps 
881   Float_t minToT = 0;   // ns
882   Float_t maxToT = 4.88;  // ns
883   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
884   arrayCal->SetOwner();
885   TFile * fileProf=0x0;
886   if(strstr(optionSave,"save")){
887     fileProf = new TFile("TOFCalibSave.root","recreate");
888   }
889   for (Int_t ich=0; ich<nch; ich++) {
890     TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
891     TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
892     Double_t binsProfile[101]; // sized larger than necessary, the correct 
893                               // dim being set in the booking of the profile
894     TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
895     Int_t ntracksTot = 0;
896     Int_t ntracks = 0;
897     Int_t nusefulbins=0;
898     Float_t meantime=0;
899     Int_t i = ch[ich];
900     AliDebug(2,Form("Calibrating channel %i",i));
901     tree->GetEntry(i);
902     ntracksTot+=nentries/3;
903     ntracks=nentries/3;
904     if (ntracks < MEANENTRIES) {
905       AliInfo(Form(" Too small mean number of entires in channel %i (number of tracks = %f), not calibrating channel and continuing.....",i,ntracks));
906       continue;
907     }
908     for (Int_t j=0;j<ntracks;j++){
909       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
910       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
911       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
912       Float_t tot = p[idxexToT];
913       hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
914       meantime+=p[idxexTime]-p[idxexExTime];
915       hToT->Fill(tot);
916     }
917
918     nusefulbins = FindBins(hToT,&binsProfile[0]);
919     meantime/=ntracksTot;
920     for (Int_t j=1;j<=nusefulbins;j++) {
921       AliDebug(2,Form(" channel %i, nusefulbins = %i, bin %i = %f",i,nusefulbins,j,binsProfile[j])); 
922     }
923
924     TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G");  // CHECK THE BUILD OPTION, PLEASE!!!!!!
925     tree->GetEntry(i);
926     ntracks=nentries/3;
927     for (Int_t j=0;j<ntracks;j++){
928       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
929       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
930       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
931       Float_t tot = p[idxexToT];
932       Float_t time = p[idxexTime]-p[idxexExTime];
933       AliDebug(2,Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
934       hSlewingProf->Fill(tot,time);  // if meantime is not used, the fill may be moved in the loop above
935       htimetot->Fill(tot,time-meantime);  // if meantime is not used, the fill may be moved in the loop above
936     }
937     
938     hSlewingProf->Fit("pol5",optionFit,"",1,4);
939     TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
940     Float_t par[6];    
941     for(Int_t kk=0;kk<6;kk++){
942       par[kk]=calibfunc->GetParameter(kk);
943       AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
944     }
945     
946     if(strstr(optionSave,"save") && fileProf){
947       TString profName=Form("Profile%06i",i);
948       TString timeTotName=Form("TimeTot%06i",i);
949       TString totName=Form("Tot%06i",i);
950       TString deltaName=Form("Delta%06i",i);
951       fileProf->cd();
952       hSlewingProf->Write(profName);
953       htimetot->Write(timeTotName);
954       hToT->Write(totName);
955       hdeltaTime->Write(deltaName);
956     }
957
958     AliTOFChannelTask * calChannel = new AliTOFChannelTask();
959     calChannel->SetSlewPar(par);
960     arrayCal->AddAt(calChannel,i);
961     delete hToT;
962     hToT=0x0;
963     delete hSlewingProf;
964     hSlewingProf=0x0;
965     delete htimetot;
966     htimetot=0x0;
967     delete hdeltaTime;
968     hdeltaTime=0x0;
969   }
970
971   if(strstr(optionSave,"save") && fileProf){
972     fileProf->Close();
973     delete fileProf;
974     fileProf=0x0;
975   }
976   file->Close();
977   delete file;
978   file=0x0;
979   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
980   arrayCal->Write("array",TObject::kSingleKey);
981   filecalib->Close();
982   delete filecalib;
983   filecalib=0x0;
984   arrayCal->Clear();
985   delete arrayCal;
986   arrayCal=0x0;
987
988   return 0;
989 }
990 //----------------------------------------------------------------------------
991 Int_t AliTOFCalibTask::CalibrateFromProfile(Int_t i, Option_t *optionSave, Option_t *optionFit){
992
993   // computing calibration parameters using the old profiling algo
994   // Returning codes:
995   // 0 -> everything was ok
996   // 1 -> no tree for calibration found
997   // 2 -> not enough statistics to perform calibration
998   // 3 -> problems with arrays
999
1000   TH1::AddDirectory(0);
1001
1002   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
1003   AliInfo(Form("*** Calibrating Histograms From Profile %s", GetName())) ; 
1004   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
1005   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
1006   const char *dirname = "./";
1007   TString filename= "TOFCalib.root";
1008
1009   if((gSystem->FindFile(dirname,filename))==NULL){
1010     AliInfo("No file with tree for calibration found! exiting....");
1011     return 1;
1012   }
1013
1014   TFile *file = new TFile("TOFCalib.root","READ");
1015   TTree *tree = (TTree*)file->Get("T");
1016   Float_t p[MAXCHENTRIESSMALL];
1017   Int_t nentries;
1018   tree->SetBranchAddress("nentries",&nentries);
1019   tree->SetBranchAddress("TOFentries",p);
1020
1021   tree->GetEntry(i);
1022   Int_t ntracks=nentries/3;
1023
1024   if (ntracks < MEANENTRIES) {  
1025     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",ntracks));
1026     return 2;
1027   }
1028
1029   TH1F * hProf = new TH1F();
1030   hProf = Profile(i);
1031   hProf->Fit("pol5",optionFit,"",0,4);
1032   TF1 * calibfunc = (TF1*)hProf->GetFunction("pol5");
1033   Float_t par[6];    
1034   for(Int_t kk=0;kk<6;kk++){
1035     par[kk]=calibfunc->GetParameter(kk);
1036     AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
1037   }
1038
1039   if(strstr(optionSave,"save")){
1040     TFile * fileProf = new TFile("TOFCalibSave.root","recreate");
1041     fileProf->cd(); 
1042     TString profName=Form("Profile%06i",i);
1043     hProf->Write(profName);
1044     fileProf->Close();
1045     delete fileProf;
1046     fileProf=0x0;
1047   }
1048
1049   delete hProf;
1050   hProf=0x0;
1051   file->Close();
1052   delete file;
1053   file=0x0;
1054   AliTOFChannelTask * calChannel = new AliTOFChannelTask();
1055   calChannel->SetSlewPar(par);
1056   arrayCal->AddAt(calChannel,i);
1057   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
1058   arrayCal->Write("array",TObject::kSingleKey);
1059   filecalib->Close();
1060   delete filecalib;
1061   filecalib = 0x0;
1062   delete calChannel;
1063   delete arrayCal;
1064   return 0;
1065 }
1066 //----------------------------------------------------------------------------
1067 Int_t AliTOFCalibTask::Calibrate(Option_t *optionSave, Option_t *optionFit){
1068
1069   // calibrating the whole TOF
1070   // computing calibration parameters
1071   // Returning codes:
1072   // 0 -> everything was ok
1073   // 1 -> no tree for calibration found
1074   // 2 -> not enough statistics to perform calibration
1075   // 3 -> problems with arrays
1076
1077   TH1::AddDirectory(0);
1078
1079   AliInfo(Form("*** Calibrating Histograms %s, all channels", GetName())) ; 
1080   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
1081   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
1082   const char *dirname = "./";
1083   TString filename= "TOFCalib.root";
1084
1085   if((gSystem->FindFile(dirname,filename))==NULL){
1086     AliInfo("No file with tree for calibration found! exiting....");
1087     return 1;
1088   }
1089
1090   TFile * fileProf=0x0;
1091   if(strstr(optionSave,"save")){
1092     fileProf = new TFile("TOFCalibSave.root","recreate");
1093   }
1094   TFile *file = new TFile("TOFCalib.root","READ");
1095   TTree *tree = (TTree*)file->Get("T");
1096   Float_t p[MAXCHENTRIESSMALL];
1097   Int_t nentries;
1098   tree->SetBranchAddress("nentries",&nentries);
1099   tree->SetBranchAddress("TOFentries",p);
1100
1101   Float_t nentriesmean =0;
1102   for (Int_t i=0; i<TOFCHANNELS; i++){
1103     tree->GetEntry(i);
1104     Int_t ntracks=nentries/3;
1105     nentriesmean+=ntracks;
1106   }
1107
1108   nentriesmean/=TOFCHANNELS;
1109   if (nentriesmean < MEANENTRIES) {
1110     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
1111     return 2;
1112   }
1113
1114   //filling ToT and Time arrays
1115
1116   Int_t nbinToT = 100;  // ToT bin width in Profile = 50.0 ps 
1117   Float_t minToT = 0;   // ns
1118   Float_t maxToT = 4.88;// ns
1119   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
1120   arrayCal->SetOwner();
1121   for (Int_t ii=0; ii<TOFCHANNELS; ii++) {
1122     TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
1123     TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
1124     TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
1125     if (ii%1000 == 0) AliDebug(1,Form("Calibrating channel %i ",ii));
1126     Int_t i = 3;
1127     Int_t nusefulbins=0;
1128     Double_t binsProfile[101]; // sized larger than necessary, the correct 
1129                               // dim being set in the booking of the profile
1130     tree->GetEntry(i);
1131     Int_t ntracks=nentries/3;
1132     if (ntracks < MEANENTRIES) {
1133       AliInfo(Form(" Too small mean number of entires in channel %i (number of tracks = %f), not calibrating channel and continuing.....",i,ntracks));
1134       continue;
1135     }
1136     Float_t meantime=0;
1137     for (Int_t j=0;j<ntracks;j++){
1138       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
1139       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
1140       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
1141       Float_t tot = p[idxexToT];
1142       hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
1143       meantime+=p[idxexTime]-p[idxexExTime];
1144       hToT->Fill(tot);
1145     }
1146     nusefulbins = FindBins(hToT,&binsProfile[0]);
1147     meantime/=ntracks;
1148     for (Int_t j=0;j<nusefulbins;j++) {
1149       AliDebug(2,Form(" channel %i, usefulbin = %i, low edge = %f",i,j,binsProfile[j])); 
1150     }
1151     TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G");  // CHECK THE BUILD OPTION, PLEASE!!!!!!
1152     for (Int_t j=0;j<ntracks;j++){
1153       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
1154       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
1155       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
1156       Float_t tot = p[idxexToT];
1157       Float_t time = p[idxexTime]-p[idxexExTime];
1158       AliDebug (2,Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
1159       hSlewingProf->Fill(tot,time-meantime);  // if meantime is not used, the fill may be moved in the loop above
1160       htimetot->Fill(tot,time-meantime);  // if meantime is not used, the fill may be moved in the loop above
1161     }
1162
1163     //    hSlewingProf->Fit("pol5",optionFit,"",1,4);
1164     //TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
1165     Float_t par[6];    
1166     for(Int_t kk=0;kk<6;kk++){
1167       //      par[kk]=calibfunc->GetParameter(kk);
1168       par[kk]=kk;
1169       AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
1170     }
1171
1172     if(strstr(optionSave,"save") && fileProf){
1173       TString profName=Form("Profile%06i",ii);
1174       TString timeTotName=Form("TimeTot%06i",ii);
1175       TString totName=Form("Tot%06i",ii);
1176       TString deltaName=Form("Delta%06i",ii);
1177       fileProf->cd();
1178       hSlewingProf->Write(profName);
1179       htimetot->Write(timeTotName);
1180       hToT->Write(totName);
1181       hdeltaTime->Write(deltaName);
1182     }
1183     AliTOFChannelTask * calChannel = new AliTOFChannelTask();
1184     calChannel->SetSlewPar(par);
1185     arrayCal->AddAt(calChannel,ii);
1186
1187     delete hToT;
1188     hToT=0x0;
1189     delete hSlewingProf;
1190     hSlewingProf=0x0;
1191     delete htimetot;
1192     htimetot=0x0;
1193     delete hdeltaTime;
1194     hdeltaTime=0x0;
1195   }
1196
1197   if(strstr(optionSave,"save")){
1198     fileProf->Close();
1199     delete fileProf;
1200     fileProf=0x0;
1201   }
1202
1203   file->Close();
1204   delete file;
1205   file=0x0;
1206   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
1207   arrayCal->Write("array", TObject::kSingleKey);
1208   filecalib->Close();
1209   delete filecalib;
1210   filecalib=0x0;
1211   arrayCal->Clear();
1212   delete arrayCal;
1213   arrayCal=0x0;
1214
1215   return 0;
1216 }
1217
1218 //_____________________________________________________________________________
1219
1220 Bool_t AliTOFCalibTask::Select(AliESDtrack *t){
1221
1222   // selecting good quality tracks
1223   //track selection: reconstrution to TOF:
1224   if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
1225     return 0;
1226   }
1227   fnESDkTOFout++;
1228   //IsStartedTimeIntegral
1229   if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
1230     return 0;
1231   }
1232   fnESDkTIME++;
1233   if (t->GetStatus() & AliESDtrack::kTRDbackup) { 
1234     Float_t xout= t->GetOuterParam()->GetX();
1235     if (xout<364.25 &&  xout > 300.) return 0;
1236   }
1237   fnESDTRDcut++;
1238   Double_t time=t->GetTOFsignal();      
1239   time*=1.E-3; // tof given in nanoseconds
1240   if(time >= fMinTime){
1241     return 0;
1242   }
1243   fnESDTIMEcut++;
1244   
1245   Double_t mom=t->GetP();
1246   if (!(mom<=UPPERMOMBOUND && mom>=LOWERMOMBOUND)){
1247   //  return 0;  // skipping momentum cut
1248   } 
1249    
1250   UInt_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
1251   if(assignedTOFcluster==0){ // not matched
1252     return 0;
1253   }
1254   fnESDassTOFcl++;
1255   return 1;
1256 }
1257 //_____________________________________________________________________________
1258
1259 Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
1260
1261   //track Combinatorial PID for calibration
1262
1263   //cout << " In comb PID " << endl;
1264   Float_t t0offset=0;
1265   const Int_t kntracksinset=6;
1266   if (kntracksinset != 6) {
1267     AliDebug(1,"Number of tracks in set smaller than expected one, identifying every particle as if it was a pion!");
1268     return 0;
1269   }
1270   Float_t exptof[kntracksinset][3];
1271   Int_t   assparticle[kntracksinset];
1272   Float_t timeofflight[kntracksinset];
1273   Float_t texp[kntracksinset];
1274   Float_t timezero[kntracksinset];
1275   Float_t weightedtimezero[kntracksinset];
1276   Float_t besttimezero[kntracksinset];
1277   Float_t bestchisquare[kntracksinset];
1278   Float_t bestweightedtimezero[kntracksinset];
1279
1280   for (Int_t i=0;i<kntracksinset;i++){
1281     assparticle[i]=3;
1282     timeofflight[i]=0;
1283     texp[i]=0;
1284     timezero[i]=0;
1285     weightedtimezero[i]=0;
1286     besttimezero[i]=0;
1287     bestchisquare[i]=0;
1288     bestweightedtimezero[i]=0;
1289   }
1290
1291   Int_t nset= (Int_t)(size/kntracksinset/NIDX);
1292   for (Int_t i=0; i< nset; i++) {   
1293     for (Int_t j=0; j<kntracksinset; j++) {
1294       Int_t idxtime = ((kntracksinset*i+j)*NIDX)+DELTAIDXTIME; 
1295       Int_t idxextimePi = ((kntracksinset*i+j)*NIDX)+DELTAIDXEXTIMEPI; 
1296       Int_t idxextimeKa = ((kntracksinset*i+j)*NIDX)+DELTAIDXEXTIMEKA; 
1297       Int_t idxextimePr = ((kntracksinset*i+j)*NIDX)+DELTAIDXEXTIMEPR; 
1298       Double_t time=smallarray[idxtime]; // TOF time in ns
1299       timeofflight[j]=time+t0offset;
1300       exptof[j][0]=smallarray[idxextimePi];
1301       exptof[j][1]=smallarray[idxextimeKa];
1302       exptof[j][2]=smallarray[idxextimePr];
1303       AliDebug(2,Form("j = %i, Time = %f, and Exp time in PID: pi = %f, K = %f, p = %f",j,timeofflight[j],exptof[j][0],exptof[j][1],exptof[j][2]));
1304     }
1305     Float_t t0best=999.;
1306     Float_t et0best=999.;
1307     Float_t chisquarebest=999.;
1308     for (Int_t i1=0; i1<3;i1++) {
1309       texp[0]=exptof[0][i1];
1310       for (Int_t i2=0; i2<3;i2++) { 
1311         texp[1]=exptof[1][i2];
1312         for (Int_t i3=0; i3<3;i3++) {
1313           texp[2]=exptof[2][i3];
1314           for (Int_t i4=0; i4<3;i4++) {
1315             texp[3]=exptof[3][i4];
1316             for (Int_t i5=0; i5<3;i5++) {
1317               texp[4]=exptof[4][i5];
1318               for (Int_t i6=0; i6<3;i6++) {
1319                 texp[5]=exptof[5][i6];
1320                                         
1321                 Float_t sumAllweights=0.;
1322                 Float_t meantzero=0.;
1323                 Float_t emeantzero=0.;
1324                 
1325                 for (Int_t itz=0; itz<kntracksinset;itz++) {
1326                   timezero[itz]=texp[itz]-timeofflight[itz];                
1327                   weightedtimezero[itz]=timezero[itz]/TRACKERROR;
1328                   sumAllweights+=1./TRACKERROR;
1329                   meantzero+=weightedtimezero[itz];
1330                   
1331                 } // end loop for (Int_t itz=0; itz<15;itz++)
1332                 
1333                 meantzero=meantzero/sumAllweights; // it is given in [ns]
1334                 emeantzero=sqrt(1./sumAllweights); // it is given in [ns]
1335                 
1336                 // calculate chisquare
1337                 
1338                 Float_t chisquare=0.;           
1339                 for (Int_t icsq=0; icsq<kntracksinset;icsq++) {
1340                   chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/TRACKERROR;
1341                 } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
1342                 
1343                 Int_t npion=0;
1344                 if(i1==0)npion++;
1345                 if(i2==0)npion++;
1346                 if(i3==0)npion++;
1347                 if(i4==0)npion++;
1348                 if(i5==0)npion++;
1349                 if(i6==0)npion++;
1350                 
1351                 if(chisquare<=chisquarebest  && ((Float_t) npion/ ((Float_t) kntracksinset)>0.3)){
1352                   for(Int_t iqsq = 0; iqsq<kntracksinset; iqsq++) {
1353                     besttimezero[iqsq]=timezero[iqsq]; 
1354                     bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
1355                     bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/TRACKERROR; 
1356                   }
1357                   
1358                   assparticle[0]=i1;
1359                   assparticle[1]=i2;
1360                   assparticle[2]=i3;
1361                   assparticle[3]=i4;
1362                   assparticle[4]=i5;
1363                   assparticle[5]=i6;
1364                   
1365                   chisquarebest=chisquare;
1366                   t0best=meantzero;
1367                   et0best=emeantzero;
1368                 } // close if(dummychisquare<=chisquare)
1369               } // end loop on i6
1370             } // end loop on i5
1371           } // end loop on i4
1372         } // end loop on i3
1373       } // end loop on i2
1374     } // end loop on i1
1375                 
1376     Float_t confLevel=999;
1377     if(chisquarebest<999.){
1378       Double_t dblechisquare=(Double_t)chisquarebest;
1379       confLevel=(Float_t)TMath::Prob(dblechisquare,kntracksinset-1); 
1380     }
1381     // assume they are all pions for fake sets
1382     if(confLevel<0.01 || confLevel==999. ){
1383       for (Int_t itrk=0; itrk<kntracksinset; itrk++)assparticle[itrk]=0;
1384     }
1385
1386     AliDebug(2,Form(" Best Assignment, selection %i %i %i %i %i %i ",assparticle[0],assparticle[1],assparticle[2],assparticle[3],assparticle[4],assparticle[5]));
1387
1388     for (Int_t kk=0;kk<kntracksinset;kk++){
1389       Int_t idxextimePi = ((kntracksinset*i+kk)*NIDX)+DELTAIDXEXTIMEPI; 
1390       Int_t idxextimeKa = ((kntracksinset*i+kk)*NIDX)+DELTAIDXEXTIMEKA; 
1391       Int_t idxextimePr = ((kntracksinset*i+kk)*NIDX)+DELTAIDXEXTIMEPR; 
1392       // storing in third slot of smallarray the assigned expected time
1393       fhPID->Fill(assparticle[kk]);
1394       assparticle[kk]=0;
1395       if (assparticle[kk]==0){       //assigned to be a Pi
1396         smallarray[idxextimeKa]=-1;
1397         smallarray[idxextimePr]=-1;
1398       }
1399       else if (assparticle[kk]==1){  //assigned to be a Ka
1400         smallarray[idxextimePi]=smallarray[idxextimeKa];
1401         smallarray[idxextimeKa]=-1;
1402         smallarray[idxextimePr]=-1;
1403       }
1404       else if (assparticle[kk]==2){  //assigned to be a Pr
1405         smallarray[idxextimePi]=smallarray[idxextimePr];
1406         smallarray[idxextimeKa]=-1;
1407         smallarray[idxextimePr]=-1;
1408       }
1409     }
1410   }
1411   return 1;
1412 }
1413 //-----------------------------------------------------------------------
1414 TH1F* AliTOFCalibTask::Profile(Int_t ich)
1415 {
1416   // profiling algo
1417
1418   TFile *file = new TFile("TOFCalib.root","READ");
1419   TTree *tree = (TTree*)file->Get("T");
1420   Float_t p[MAXCHENTRIESSMALL];
1421   Int_t nentries;
1422   tree->SetBranchAddress("nentries",&nentries);
1423   tree->SetBranchAddress("TOFentries",p);
1424
1425   //Prepare histograms for Slewing Correction
1426   const Int_t knbinToT = 100;
1427   Int_t nbinTime = 200;
1428   Float_t minTime = -5.5; //ns
1429   Float_t maxTime = 5.5; //ns
1430   Float_t minToT = 0; //ns
1431   Float_t maxToT = 5.; //ns
1432   Float_t deltaToT = (maxToT-minToT)/knbinToT;
1433   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];
1434   Int_t n[knbinToT+1], nentrx[knbinToT+1];
1435   Double_t sigmaToT[knbinToT+1];
1436   for (Int_t i = 0; i < knbinToT+1 ; i++){
1437     mTime[i]=0;
1438     mToT[i]=0;
1439     n[i]=0;
1440     meanTime[i]=0;
1441     meanTime2[i]=0;
1442     vToT[i]=0;
1443     vToT2[i]=0;
1444     meanToT[i]=0;
1445     meanToT2[i]=0;
1446     vTime[i]=0;
1447     vTime2[i]=0;
1448     xlow[i]=0;
1449     sigmaTime[i]=0;
1450     sigmaToT[i]=0;
1451     n[i]=0;
1452     nentrx[i]=0;
1453   }
1454   TH2F* hSlewing = new TH2F("hSlewing", "hSlewing", knbinToT, minToT, maxToT, nbinTime, minTime, maxTime);
1455   Int_t ntracks = 0;
1456   TH1F *histo = new TH1F("histo", "1D Time vs ToT", knbinToT, minToT, maxToT);
1457   tree->GetEntry(ich);
1458   ntracks=nentries/3;
1459   for (Int_t j=0;j<ntracks;j++){
1460     Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
1461     Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
1462     Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
1463     Float_t tot = p[idxexToT];
1464     Float_t time = p[idxexTime]-p[idxexExTime];
1465     Int_t nx = (Int_t)((tot-minToT)/deltaToT)+1;
1466     if ((tot != 0) && ( time!= 0)){
1467       vTime[nx]+=time;
1468       vTime2[nx]+=time*time;
1469       vToT[nx]+=tot;
1470       vToT2[nx]+=tot*tot;
1471       nentrx[nx]++;
1472       hSlewing->Fill(tot,time);
1473     }
1474   }
1475   Int_t nbinsToT=hSlewing->GetNbinsX();
1476   if (nbinsToT != knbinToT) {
1477     AliError("Profile :: incompatible numbers of bins");
1478     return 0x0;
1479   }
1480   
1481   Int_t usefulBins=0;
1482   for (Int_t i=1;i<=nbinsToT;i++){
1483     if (nentrx[i]!=0){
1484       n[usefulBins]+=nentrx[i];
1485       if (n[usefulBins]==0 && i == nbinsToT) {
1486         break;
1487       }
1488       meanTime[usefulBins]+=vTime[i];
1489       meanTime2[usefulBins]+=vTime2[i];
1490       meanToT[usefulBins]+=vToT[i];
1491       meanToT2[usefulBins]+=vToT2[i];
1492       if (n[usefulBins]<10 && i!=nbinsToT) continue; 
1493       mTime[usefulBins]=meanTime[usefulBins]/n[usefulBins];
1494       mToT[usefulBins]=meanToT[usefulBins]/n[usefulBins];
1495       sigmaTime[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
1496                             *(meanTime2[usefulBins]-meanTime[usefulBins]
1497                             *meanTime[usefulBins]/n[usefulBins]));
1498       if ((1./n[usefulBins]/n[usefulBins]
1499            *(meanToT2[usefulBins]-meanToT[usefulBins]
1500              *meanToT[usefulBins]/n[usefulBins]))< 0) {
1501         AliError(" too small radical" );
1502         sigmaToT[usefulBins]=0;
1503       }
1504       else{       
1505         sigmaToT[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
1506                              *(meanToT2[usefulBins]-meanToT[usefulBins]
1507                              *meanToT[usefulBins]/n[usefulBins]));
1508       }
1509       usefulBins++;
1510     }
1511   }
1512   for (Int_t i=0;i<usefulBins;i++){
1513     Int_t binN = (Int_t)((mToT[i]-minToT)/deltaToT)+1;
1514     histo->Fill(mToT[i],mTime[i]);
1515     histo->SetBinError(binN,sigmaTime[i]);
1516   } 
1517   delete file;
1518   file=0x0;
1519   delete hSlewing;
1520   hSlewing=0x0;
1521
1522   return histo;
1523 }
1524 //----------------------------------------------------------------------------
1525 Int_t AliTOFCalibTask::FindBins(TH1F* h, Double_t *binsProfile) const{
1526
1527   // to determine the bins for ToT histo
1528
1529   Int_t cont = 0;
1530   Int_t startBin = 1;
1531   Int_t nbin = h->GetNbinsX();
1532   Float_t max = h->GetBinLowEdge(nbin);
1533   Int_t nusefulbins=0;
1534   Int_t maxcont=0;
1535   // setting maxvalue of entries per bin
1536   if (nbin <= 60) maxcont = 2;
1537   else  if (nbin <= 100) maxcont = 5;
1538   else  if (nbin <= 500) maxcont = 10;
1539   else  maxcont = 20;
1540   for (Int_t j=1;j<=nbin;j++) {
1541     cont += (Int_t)h->GetBinContent(j);
1542     if (j<nbin){
1543       if (cont>=maxcont){
1544         nusefulbins++;
1545         binsProfile[nusefulbins-1]=h->GetBinLowEdge(startBin);
1546         cont=0;
1547         startBin=j+1;
1548         continue;
1549       }
1550     }
1551     else{
1552       if (cont>=maxcont){
1553         nusefulbins++;
1554         binsProfile[nusefulbins-1]=h->GetBinLowEdge(startBin);
1555         binsProfile[nusefulbins]=max;
1556       }
1557       else {
1558         binsProfile[nusefulbins]=h->GetBinLowEdge(startBin);
1559       }
1560     }
1561   }
1562   return nusefulbins;
1563 }