corrected TRD/TOF MV position
[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 Revision 1.6  2007/10/05 10:44:11  zampolli
18 Oversight for debugging purpose fixed
19
20 Revision 1.5  2007/10/05 10:38:10  zampolli
21 Oversight fixed
22
23 Revision 1.4  2007/10/04 15:36:44  zampolli
24 Updates to new TOF offline calibration schema
25
26 Revision 1.3  2007/07/31 07:26:16  zampolli
27 Bug fixed in the $ field
28
29 */
30
31 // task to perform TOF calibration
32 // optimized for pp runs
33 // expect a maximum of 100 entries per channel
34 // different ways to calibrate are implemented
35 // C. Zampolli 
36
37 #include <TObject.h> 
38 #include <TCanvas.h>
39 #include <TROOT.h>
40 #include <TSystem.h>
41 #include <TProfile.h>
42 #include <TF1.h>
43 #include <TTree.h>
44 #include <TFile.h>
45 #include <TH1F.h>
46 #include <TH1I.h>
47 #include <TH1D.h>
48 #include <TH2F.h>
49 #include <TGrid.h>
50 #include <TList.h>
51 #include <TArrayF.h>
52 #include <TBenchmark.h>
53   
54 #include "AliTOFCalibTask.h" 
55 #include "AliESDEvent.h" 
56 #include "AliESDtrack.h" 
57 #include "AliLog.h"
58 #include "AliTOFArray.h"
59
60 //_________________________________________________________________________
61 AliTOFCalibTask::AliTOFCalibTask(const char *name) : 
62   AliAnalysisTaskSE(name),  
63   fESD(0),
64   fToT(0),
65   fTime(0),
66   fExpTimePi(0),
67   fExpTimeKa(0),
68   fExpTimePr(0),
69   fMinTime(0),
70   fTOFArray(0x0),
71   fnESD(0),
72   fnESDselected(0),
73   fnESDkTOFout(0),
74   fnESDkTIME(0),
75   fnESDassTOFcl(0),
76   fnESDTIMEcut(0),
77   fnESDTRDcut(0),
78   fhToT(0),
79   fhTime(0),
80   fhExpTimePi(0),
81   fhExpTimeKa(0),
82   fhExpTimePr(0),
83   fhPID(0),
84   fhch(0),
85   fhESD(0),
86   fhESDselected(0),
87   fhESDkTOFout(0),
88   fhESDkTIME(0),
89   fhESDassTOFcl(0),
90   fhESDTIMEcut(0),
91   fhESDTRDcut(0),
92   fListOfHistos(0x0),
93   fListArray(0x0)
94 {
95   // Constructor.
96
97   DefineOutput(1,TList::Class()); 
98   DefineOutput(2,TList::Class()); 
99
100   for (Int_t i=0;i<11;i++){
101     fassparticle[i]=-1;
102   } 
103
104 }
105
106 //______________________________________________________________________________
107 AliTOFCalibTask::AliTOFCalibTask(const AliTOFCalibTask &calibtask) : 
108   AliAnalysisTaskSE("AliTOFCalibTask"),  
109   fESD(calibtask.fESD),
110   fToT(calibtask.fToT),
111   fTime(calibtask.fTime),
112   fExpTimePi(calibtask.fExpTimePi),
113   fExpTimeKa(calibtask.fExpTimeKa),
114   fExpTimePr(calibtask.fExpTimePr),
115   fMinTime(calibtask.fMinTime),
116   fTOFArray(calibtask.fTOFArray),
117   fnESD(calibtask.fnESD),
118   fnESDselected(calibtask.fnESDselected),
119   fnESDkTOFout(calibtask.fnESDkTOFout),
120   fnESDkTIME(calibtask.fnESDkTIME),
121   fnESDassTOFcl(calibtask.fnESDassTOFcl),
122   fnESDTIMEcutcalibtask.fnESDTIMEcut(),
123   fnESDTRDcutcalibtask.fnESDTRDcut(),
124   fhToT(calibtask.fhToT),
125   fhTime(calibtask.fhTime),
126   fhExpTimePi(calibtask.fhExpTimePi),
127   fhExpTimeKa(calibtask.fhExpTimeKa),
128   fhExpTimePr(calibtask.fhExpTimePr),
129   fhPID(calibtask.fhPID),
130   fhch(calibtask.fhch),
131   fhESD(calibtask.fhESD),
132   fhESDselected(calibtask.fhESDselected),
133   fhESDkTOFout(calibtask.fhESDkTOFout),
134   fhESDkTIME(calibtask.fhESDkTIME),
135   fhESDassTOFcl(calibtask.fhESDassTOFcl),
136   fhESDTIMEcut(calibtask.fhESDTIMEcut),
137   fhESDTRDcut(calibtask.fhESDTRDcut),
138   fListOfHistos(calibtask.fListOfHistos),
139   fListArray(calibtask.fListArray)
140
141 {
142   // Copy Constructor.
143
144   for (Int_t i=0;i<11;i++){
145     fassparticle[i]=calibtask.fassparticle[i];
146   } 
147
148 }
149 //______________________________________________________________________________
150 AliTOFCalibTask:: ~AliTOFCalibTask() 
151 {
152   // destructor
153
154   AliInfo("TOF Calib Task: Deleting");
155   
156   delete fTOFArray;  
157   delete fhToT;
158   delete fhTime;
159   delete fhExpTimePi; 
160   delete fhExpTimeKa; 
161   delete fhExpTimePr; 
162   delete fhPID; 
163   delete fhch;
164   delete fhESD;
165   delete fhESDselected;
166   delete fhESDkTOFout;
167   delete fhESDkTIME;
168   delete fhESDTRDcut;
169   delete fhESDTIMEcut;
170   delete fhESDassTOFcl;
171   //  delete fListOfHistos;
172 }
173 //______________________________________________________________________________
174 AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)  
175
176    //assignment operator
177   fESD=calibtask.fESD;
178   fToT=calibtask.fToT;
179   fTime=calibtask.fTime;
180   fExpTimePi=calibtask.fExpTimePi;
181   fExpTimeKa=calibtask.fExpTimeKa;
182   fExpTimePr=calibtask.fExpTimePr;
183   fMinTime=calibtask.fMinTime;
184   fTOFArray=calibtask.fTOFArray;
185   fnESD=calibtask.fnESD;
186   fnESDselected=calibtask.fnESDselected;
187   fnESDkTOFout=calibtask.fnESDkTOFout;
188   fnESDkTIME=calibtask.fnESDkTIME;
189   fnESDassTOFcl=calibtask.fnESDassTOFcl;
190   fnESDTIMEcut=calibtask.fnESDTIMEcut;
191   fnESDTRDcut=calibtask.fnESDTRDcut;
192   fhToT=calibtask.fhToT;
193   fhTime=calibtask.fhTime;
194   fhExpTimePi=calibtask.fhExpTimePi;
195   fhExpTimeKa=calibtask.fhExpTimeKa;
196   fhExpTimePr=calibtask.fhExpTimePr;
197   fhPID=calibtask.fhPID;
198   fhch=calibtask.fhch;
199   fhESD=calibtask.fhESD;
200   fhESDselected=calibtask.fhESDselected;
201   fhESDkTOFout=calibtask.fhESDkTOFout;
202   fhESDkTIME=calibtask.fhESDkTIME;
203   fhESDassTOFcl=calibtask.fhESDassTOFcl;
204   fhESDTIMEcut=calibtask.fhESDTIMEcut;
205   fhESDTRDcut=calibtask.fhESDTRDcut;
206   fListOfHistos=calibtask.fListOfHistos;
207   fListArray=calibtask.fListArray;
208
209   for (Int_t i=0;i<11;i++){
210     fassparticle[i]=calibtask.fassparticle[i];
211   } 
212
213   return *this;
214 }
215 //--------------------------------------------------------------------------
216 void AliTOFCalibTask::BookHistos(){
217
218   //booking histos
219
220   AliInfo(Form("*** Booking Histograms %s", GetName())) ; 
221
222   fhToT=
223     new TH1F("hToT", " ToT distribution (ns) ", 400, 0, 40);
224   fhTime=
225     new TH1F("hTime", " Time distribution (ns) ", 400, 0, 40);
226   fhExpTimePi=
227     new TH1F("hExpTimePi", " Exp Time distribution, Pi (ns) ", 400, 0, 40);
228   fhExpTimeKa=
229     new TH1F("hExpTimeKa", " Exp Time distribution, Ka (ns) ", 400, 0, 40);
230   fhExpTimePr=
231     new TH1F("hExpTimePr", " Exp Time distribution, Pr (ns) ", 400, 0, 40);
232   fhPID=
233     new TH1I("hPID", " Combinatorial PID identity ", 3, 0, 3);
234   fhch=
235     new TH1D("hch", " TOF channel ", TOFCHANNELS, 0, TOFCHANNELS);
236
237
238   // create the output list of histos
239   if (!fListOfHistos) fListOfHistos = new TList();
240   fListOfHistos->SetOwner();
241
242   fListOfHistos->AddAt(fhToT,             0) ; 
243   fListOfHistos->AddAt(fhTime,            1) ; 
244   fListOfHistos->AddAt(fhExpTimePi,       2) ; 
245   fListOfHistos->AddAt(fhExpTimeKa,       3) ; 
246   fListOfHistos->AddAt(fhExpTimePr,       4) ; 
247   fListOfHistos->AddAt(fhPID,             5) ; 
248   fListOfHistos->AddAt(fhch,              6) ; 
249
250   fhESD=
251     new TH1I("hESD","Number of analyzed ESDs",1,0,1);
252   fhESDselected=
253     new TH1I("hESDselected","Number of selected ESDs",1,0,1);
254   fhESDkTOFout=
255     new TH1I("hESDkTOFout","Number of ESDs with kTOFout",1,0,1);
256   fhESDkTIME=
257     new TH1I("hESDkTIME","Number of ESDs with kTime",1,0,1);
258   fhESDTRDcut=
259     new TH1I("hESDTRDcut","Number of ESDs with TRDcut",1,0,1);
260   fhESDTIMEcut=
261     new TH1I("hESDTIMEcut","Number of ESDs with TIMEcut",1,0,1);
262   fhESDassTOFcl=
263     new TH1I("hESDassTOFcl","Number of ESDs with assTOFcl",1,0,1);
264
265   fListOfHistos->AddAt(fhESD,             7) ; 
266   fListOfHistos->AddAt(fhESDselected,     8) ; 
267   fListOfHistos->AddAt(fhESDkTOFout,      9) ; 
268   fListOfHistos->AddAt(fhESDkTIME,        10) ; 
269   fListOfHistos->AddAt(fhESDTRDcut,       11) ; 
270   fListOfHistos->AddAt(fhESDTIMEcut,      12) ; 
271   fListOfHistos->AddAt(fhESDassTOFcl,     13) ; 
272
273   return;
274   
275 }
276 //----------------------------------------------------------------------------
277 void AliTOFCalibTask::DrawHistos(){
278
279   // drawing output histos
280
281   AliInfo(Form("*** Drawing Histograms %s", GetName())) ; 
282
283   if (!gROOT->IsBatch()){
284
285           fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
286           
287           fhToT = (TH1F*)fListOfHistos->At(0);
288           fhTime = (TH1F*)fListOfHistos->At(1);
289           fhExpTimePi = (TH1F*)fListOfHistos->At(2);
290           fhExpTimeKa = (TH1F*)fListOfHistos->At(3);
291           fhExpTimePr = (TH1F*)fListOfHistos->At(4);
292           fhPID = (TH1I*)fListOfHistos->At(5);
293           fhch = (TH1D*)fListOfHistos->At(6);
294           
295           fhESD = (TH1I*)fListOfHistos->At(7);
296           fhESDselected = (TH1I*)fListOfHistos->At(8);
297           fhESDkTOFout = (TH1I*)fListOfHistos->At(9);
298           fhESDkTIME = (TH1I*)fListOfHistos->At(10);
299           fhESDassTOFcl = (TH1I*)fListOfHistos->At(11);
300           fhESDTIMEcut = (TH1I*)fListOfHistos->At(12);
301           fhESDTRDcut = (TH1I*)fListOfHistos->At(13);
302           
303           TCanvas * canvasToTTime = new TCanvas("canvasToTTime", " ToT and Time ",400, 30, 550, 630) ;
304           canvasToTTime->Divide(1,2);
305           canvasToTTime->cd(1);
306           fhToT->SetLineColor(4);
307           fhToT->GetXaxis()->SetTitle("ToT (ns)");
308           fhToT->DrawCopy("hist");
309           canvasToTTime->cd(2);
310           fhTime->SetLineColor(4);
311           fhTime->GetXaxis()->SetTitle("Time (ns)");
312           fhTime->DrawCopy("hist");
313           canvasToTTime->Update();
314           //canvasToTTime->Print("ToTTime.gif");
315           
316           TCanvas * canvasExpTime = new TCanvas("canvasExpTime", " Expected Times ",400, 30, 550, 630) ;
317           canvasExpTime->Divide(1,3);
318           canvasExpTime->cd(1);
319           fhExpTimePi->SetLineColor(4);
320           fhExpTimePi->GetXaxis()->SetTitle("Exp Time (ns), #pi");
321           fhExpTimePi->DrawCopy("hist");
322           canvasExpTime->cd(2);
323           fhExpTimeKa->SetLineColor(4);
324           fhExpTimeKa->GetXaxis()->SetTitle("Exp Time (ns), K");
325           fhExpTimeKa->DrawCopy("hist");
326           canvasExpTime->cd(3);
327           fhExpTimePr->SetLineColor(4);
328           fhExpTimePr->GetXaxis()->SetTitle("Exp Time (ns), p");
329           fhExpTimePr->DrawCopy("hist");
330           
331           //canvasExpTime->Print("ExpTime.gif");
332           
333           TCanvas * canvasPID = new TCanvas("canvasPID", " Combinatorial PID ",400, 30, 550, 400);
334           canvasPID->cd();
335           fhPID->GetXaxis()->SetTitle("Comb PID");
336           fhPID->GetXaxis()->SetBinLabel(1,"#pi");
337           fhPID->GetXaxis()->SetBinLabel(2,"K");
338           fhPID->GetXaxis()->SetBinLabel(3,"p");
339           fhPID->DrawCopy("hist");
340           
341           //canvasPID->Print("PID.gif");
342           
343           TCanvas * canvasrndch = new TCanvas("canvasrndch", " TOF channel ",400, 30, 550, 400);
344           canvasrndch->cd();
345           fhch->GetXaxis()->SetTitle("TOF ch");
346           fhch->Draw("hist");
347           Float_t meanTOFch = 0;
348           for (Int_t ibin=0;ibin<TOFCHANNELS;ibin++){
349                   meanTOFch+=(Float_t)fhch->GetBinContent(ibin+1);
350           }
351           
352           meanTOFch/=TOFCHANNELS;
353           AliDebug(1,Form(" Mean number of tracks/channel = %f ",meanTOFch));
354           
355           //canvasrndch->Print("rndch.gif");
356
357           /*
358           char line[1024] ; 
359           sprintf(line, ".!tar -zcvf %s.tar.gz *.gif", GetName()) ; 
360           gROOT->ProcessLine(line);
361           sprintf(line, ".!rm -fR *.gif"); 
362           gROOT->ProcessLine(line);
363           AliInfo(Form("*** TOF Calib Task: plots saved in %s.tar.gz...\n", GetName())) ;
364           */
365           
366           AliInfo(Form(" Number of analyzed ESD tracks: %i",(Int_t)fhESD->GetEntries()));
367           AliInfo(Form(" Number of selected ESD tracks: %i",(Int_t)fhESDselected->GetEntries()));
368           AliInfo(Form(" Number of ESD tracks with kTOFout: %i",(Int_t)fhESDkTOFout->GetEntries()));
369           AliInfo(Form(" Number of ESD tracks with kTIME: %i",(Int_t)fhESDkTIME->GetEntries()));
370           AliInfo(Form(" Number of ESD tracks with TRDcut: %i",(Int_t)fhESDTRDcut->GetEntries()));
371           AliInfo(Form(" Number of ESD tracks with TIMEcut: %i",(Int_t)fhESDTIMEcut->GetEntries()));
372           AliInfo(Form(" Number of ESD tracks with assTOFcl: %i",(Int_t)fhESDassTOFcl->GetEntries()));
373   }
374   return;
375 }
376
377 //________________________________________________________________________
378 void AliTOFCalibTask::UserCreateOutputObjects()
379 {
380 // Create histograms and AliTOFArray
381
382         AliInfo("UserCreateObjects");
383         OpenFile(1);
384         BookHistos();
385
386         if (!fTOFArray) fTOFArray = new AliTOFArray(TOFCHANNELS);
387         for (Int_t i =0;i<TOFCHANNELS;i++){
388                 fTOFArray->SetArray(i);
389         }
390         if (!fListArray) fListArray = new TList();
391         fListArray->SetOwner();
392         fListArray->Add(fTOFArray);
393         return;
394 }
395
396 //______________________________________________________________________________
397 void AliTOFCalibTask::UserExec(Option_t * /*opt*/) 
398 {
399
400   // main 
401
402
403 //******* The loop over events -----------------------------------------------
404
405   AliVEvent*    fESD = fInputEvent ;
406   if (!fESD) {
407     Error("UserExec","NO EVENT FOUND!");
408     return;
409   }
410
411
412   fMinTime=22E3;   //ns; not used
413   Int_t ntrk = fESD->GetNumberOfTracks() ;
414   fnESD+=ntrk;
415   Int_t nselected = 0;
416   Int_t itr = -1;
417   while ( ntrk-- ) {
418     fhESD->Fill(0);
419     itr++;
420     AliESDtrack * t = (AliESDtrack*)fESD->GetTrack(ntrk) ;
421     //selecting only good quality tracks
422     if (!Select(t)) continue;
423     nselected++;
424     fhESDselected->Fill(0);
425     Int_t ich = Int_t(t->GetTOFCalChannel()); 
426     fhch->Fill(ich);
427     //    ich=3; //only for debug purpose
428     AliDebug(2,Form(" ESD in channel %i, filling array %i",t->GetTOFCalChannel(),ich));
429     Float_t tot = t->GetTOFsignalToT();
430     Float_t time = t->GetTOFsignalRaw();
431     AliDebug(2,Form(" track # %i in channel %i, time = %f \n",ntrk,ich,time)); 
432     Double_t expTime[AliPID::kSPECIESC]; 
433     t->GetIntegratedTimes(expTime,AliPID::kSPECIESC);
434     Float_t expTimePi = expTime[2]*1.E-3;
435     Float_t expTimeKa = expTime[3]*1.E-3;
436     Float_t expTimePr = expTime[4]*1.E-3;
437     Int_t currentSize = fTOFArray->GetArraySize(ich);
438     if (currentSize==CHENTRIES){ 
439             AliDebug(2,Form("too many tracks in channel %i, not storing any more...",ich));
440       continue;
441     } 
442     AliDebug(2,Form("tracks in channel %i = %i, storing... ",ich, currentSize/NIDX ));
443
444     fTOFArray->ReSetArraySize(ich,currentSize+5);
445     fTOFArray->SetAt(ich,currentSize+DELTAIDXTOT,tot);
446     fTOFArray->SetAt(ich,currentSize+DELTAIDXTIME,time*1E-3);
447     fTOFArray->SetAt(ich,currentSize+DELTAIDXEXTIMEPI,expTimePi);
448     fTOFArray->SetAt(ich,currentSize+DELTAIDXEXTIMEKA,expTimeKa);
449     fTOFArray->SetAt(ich,currentSize+DELTAIDXEXTIMEPR,expTimePr);
450     fhToT->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXTOT));
451     fhTime->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXTIME));
452     fhExpTimePi->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXEXTIMEPI));
453     fhExpTimeKa->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXEXTIMEKA));
454     fhExpTimePr->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXEXTIMEPR));
455     AliDebug(2,Form("track = %i, tot = %f, time = %f, and Exp time in TOF: pi = %f, K = %f, p = %f",itr, fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXTOT), fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXTIME), expTimePi,expTimeKa,expTimePr));
456
457
458   }
459   fnESDselected+=nselected;
460   PostData(1, fListOfHistos);  
461   PostData(2, fListArray);  
462
463 }
464
465 //_____________________________________________________________________________
466 void AliTOFCalibTask::Terminate(Option_t *)
467 {
468   // Processing when the event loop is ended
469   
470   // some plots
471
472   TH1::AddDirectory(0);
473   AliInfo("TOF Calib Task: End of events loop");
474   DrawHistos();
475   fListArray = dynamic_cast<TList*>(GetOutputData(2));
476   fTOFArray = dynamic_cast<AliTOFArray*>(fListArray->At(0));
477
478   for (Int_t ich = 0;ich<TOFCHANNELS;ich++){
479           Int_t entriesChannel = fTOFArray->GetArraySize(ich)/NIDX;
480           if (entriesChannel>0){
481                   Int_t ncuttime = SelectOnTime(fTOFArray->GetArray(ich),entriesChannel,ich);
482                   fnESDselected-=ncuttime;
483           }
484   }
485
486   TBenchmark bench;
487   bench.Start("CombPID");
488   for (Int_t i = 0;i<TOFCHANNELS;i++){
489     Int_t size=fTOFArray->GetArraySize(i);
490     AliDebug(2, Form(" entries %i in channel %i ",size/NIDX,i));
491     if (size/NIDX<=2) {
492       AliDebug(1, Form(" not enough statistics for combined PID for channel %i, putting all the tracks as if they were pions",i));
493       continue;
494     }
495     if (i%1000 == 0) AliInfo(Form("At channel %d",i));
496     if (!CombPID(fTOFArray->GetArray(i), size)) AliError("ERROR!!!!ERROR!!!");
497   }
498   bench.Stop("CombPID");
499   bench.Print("CombPID");
500
501   // saving data in a tree --> obsolete code; keeping for backup with new structure
502   // using AliTOFArray
503   /*  
504   AliInfo("Building tree for Calibration");
505   TTree * tree = new TTree("T", "Tree for TOF calibration");
506   AliTOFArray* tempArray = new AliTOFArray(TOFCHANNELS);
507   tree->Branch("TOFentries","AliTOFArray",&tempArray,32000,0);
508   for (Int_t i = 0;i<TOFCHANNELS;i++){
509           Int_t sizeChannel = fTOFArray->GetArraySize(i)/NIDX;
510           if (i==3) AliDebug(2,Form("Entries in channel %d = %d ",i,sizeChannel));  // just for debug 
511           if (sizeChannel!=0){
512                   tempArray->SetArray(i,NIDXSMALL*sizeChannel);
513           }
514           for (Int_t j =0; j<sizeChannel;j++){
515                   for (Int_t k=0; k<NIDXSMALL;k++){
516                           tempArray->SetAt(i,j*NIDXSMALL+k,fTOFArray->GetArrayAt(i,j*NIDX+k));
517                   }
518           }
519   }
520   tree->Fill();
521   */      
522   
523   AliInfo("Putting object for calibration in Reference data");
524
525   TFile *fileout = TFile::Open("outputTOFCalibration.root","RECREATE");
526   tempArray->Write();
527   fileout->Close();
528   delete tempArray;
529 }
530 //_____________________________________________________________________________
531
532 Bool_t AliTOFCalibTask::Select(AliESDtrack *t){
533
534   // selecting good quality tracks
535   //track selection: reconstrution to TOF:
536   if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
537     return 0;
538   }
539   fnESDkTOFout++;
540   fhESDkTOFout->Fill(0);
541   //IsStartedTimeIntegral
542   if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
543     return 0;
544   }
545   fnESDkTIME++;
546   fhESDkTIME->Fill(0);
547   if (t->GetStatus() & AliESDtrack::kTRDbackup) { 
548     Float_t xout= t->GetOuterParam()->GetX();
549     if (xout<364.25 &&  xout > 300.) return 0;
550   }
551   fnESDTRDcut++;
552   fhESDTRDcut->Fill(0);
553   Double_t time=t->GetTOFsignal();      
554   time*=1.E-3; // tof given in nanoseconds
555   if(time >= fMinTime){
556     return 0;
557   }
558   fnESDTIMEcut++;
559   fhESDTIMEcut->Fill(0);
560   
561   Double_t mom=t->GetP();
562   if (!(mom<=UPPERMOMBOUND && mom>=LOWERMOMBOUND)){
563   //  return 0;  // skipping momentum cut
564   } 
565    
566   Int_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
567   if(assignedTOFcluster==-1){ // not matched
568     return 0;
569   }
570   fnESDassTOFcl++;
571   fhESDassTOFcl->Fill(0);
572   AliDebug(2,"selecting the track");
573   return 1;
574 }
575 //_____________________________________________________________________________
576
577 Int_t AliTOFCalibTask::SelectOnTime(Float_t *charray, Int_t ntracks, Int_t ich){
578
579         // to be re-implemented with new object
580
581   // discarding tracks with time-mintime < MINTIME
582   
583   Int_t ndeleted=0;
584   Float_t mintime = 1E6;
585   for (Int_t itr = 0;itr<ntracks;itr++){
586     Int_t ientry=itr*NIDX;
587     Float_t time = charray[ientry+DELTAIDXTIME];// in ns
588       if (time<mintime) mintime = time;
589   }
590   AliDebug(1,Form("Mintime for channel %i = %f",ich, mintime));
591   /*
592   Int_t nTracksInChannel = fTOFArray->GetArray(ich)->GetSize();
593   for (Int_t itr = 0;itr<ntracks;itr++){
594     Int_t ientry=itr*NIDX;
595     Float_t time = charray[ientry+DELTAIDXTIME];// in ns
596     if ((time-mintime)>MINTIME) {
597       ndeleted++;
598       AliDebug(1,Form("Deleting %i track from channel %i, time = %f",ndeleted, ich, time));
599       nTracksInChannel;
600       for (Int_t j=itr+1;j<ntracks;j++){
601         Int_t ientry=j*NIDX;
602         Int_t idxtot = ientry+DELTAIDXTOT; 
603         Int_t idxtime = ientry+DELTAIDXTIME; 
604         Int_t idxextimePi = ientry+DELTAIDXEXTIMEPI; 
605         Int_t idxextimeKa = ientry+DELTAIDXEXTIMEKA; 
606         Int_t idxextimePr = ientry+DELTAIDXEXTIMEPR;
607         Int_t ientrydel=(j-1)*NIDX;
608         Int_t idxtotdel = ientrydel+DELTAIDXTOT; 
609         Int_t idxtimedel = ientrydel+DELTAIDXTIME; 
610         Int_t idxextimePidel = ientrydel+DELTAIDXEXTIMEPI; 
611         Int_t idxextimeKadel = ientrydel+DELTAIDXEXTIMEKA; 
612         Int_t idxextimePrdel = ientrydel+DELTAIDXEXTIMEPR;
613         charray[idxtotdel]=charray[idxtot]; 
614         charray[idxtimedel]=charray[idxtime]; 
615         charray[idxextimePidel]=charray[idxextimePi]; 
616         charray[idxextimeKadel]=charray[idxextimeKa]; 
617         charray[idxextimePrdel]=charray[idxextimePr];
618       }
619     } 
620   }
621   */
622   return ndeleted;
623 }
624 //_____________________________________________________________________________
625
626 Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
627
628   // track Combinatorial PID for calibration, only when more than 2 particles
629   // fall in channel
630
631   Float_t t0offset=0;
632
633   if (size/NIDX<=2){
634     AliDebug(1,"Number of tracks in channel smaller than 2, identifying every particle as if it was a pion!");
635     return 0;
636   }
637
638   Int_t ntracksinchannel = size/NIDX;
639   Int_t ntrkinset=-1;
640   Int_t nset = ntracksinchannel/6;
641
642   if (nset ==0) {
643     nset=1;
644     if (ntracksinchannel < 6) {
645       AliDebug(2,"Number of tracks in set smaller than 6, Combinatorial PID still applied to this set.");
646     }
647   }
648
649   AliDebug(2,Form("number of sets = %i", nset));
650   // loop on sets
651   for (Int_t i=0; i< nset; i++) {   
652     if (i<nset-1)ntrkinset=6;
653     else {
654       if (ntracksinchannel<6){
655         ntrkinset=ntracksinchannel;
656       }
657       else{
658         ntrkinset=6+Int_t((ntracksinchannel)%6);
659       }
660     }
661     AliDebug(2,Form("set = %i, number of tracks in set = %i", i,ntrkinset));
662       
663     for (Int_t ii=0;ii<ntrkinset;ii++){
664       fassparticle[ii]=-1;
665     }
666     Float_t **exptof;
667     Float_t* timeofflight;
668     Float_t* texp;
669     Int_t* index;
670
671     exptof = new Float_t*[ntrkinset]; 
672     timeofflight = new Float_t[ntrkinset];
673     texp = new Float_t[ntrkinset];
674     index = new Int_t[ntrkinset];
675     for (Int_t ii=0;ii<ntrkinset;ii++){
676       exptof[ii] = new Float_t[3];
677       timeofflight[ii]=0;
678       texp[ii]=0;
679       index[ii]=-1;
680     }
681
682     for (Int_t j=0; j<ntrkinset; j++) {
683       Int_t idxtime = ((6*i+j)*NIDX)+DELTAIDXTIME; 
684       Int_t idxextimePi = ((6*i+j)*NIDX)+DELTAIDXEXTIMEPI; 
685       Int_t idxextimeKa = ((6*i+j)*NIDX)+DELTAIDXEXTIMEKA; 
686       Int_t idxextimePr = ((6*i+j)*NIDX)+DELTAIDXEXTIMEPR; 
687       AliDebug(2,Form("idxtime = %i, idxextimePi = %i, idxextimeKa = %i, idxextimePr = %i", idxtime, idxextimePi, idxextimeKa, idxextimePr));
688       Double_t time=smallarray[idxtime]; // TOF time in ns
689       timeofflight[j]=time+t0offset;
690       exptof[j][0]=smallarray[idxextimePi];
691       exptof[j][1]=smallarray[idxextimeKa];
692       exptof[j][2]=smallarray[idxextimePr];
693       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]));
694     }
695     
696     
697     Float_t chisquarebest=999.;
698     AliDebug(2,Form(" Set = %i with %i tracks ", i,ntrkinset));
699     chisquarebest = LoopCombPID(ntrkinset, ntrkinset,exptof,&texp[0],&timeofflight[0], &index[0],chisquarebest); 
700     
701     Float_t confLevel=999;
702     if(chisquarebest<999.){
703       Double_t dblechisquare=(Double_t)chisquarebest;
704       confLevel=(Float_t)TMath::Prob(dblechisquare,ntrkinset-1); 
705     }
706
707     // assume they are all pions for fake sets
708     if(confLevel<0.01 || confLevel==999. ){
709       for (Int_t itrk=0; itrk<ntrkinset; itrk++)fassparticle[itrk]=0;
710     }
711     
712     AliDebug(2,Form(" Best Assignment, selection %i %i %i %i %i %i %i %i %i %i",fassparticle[0],fassparticle[1],fassparticle[2],fassparticle[3],fassparticle[4],fassparticle[5],fassparticle[6],fassparticle[7],fassparticle[8],fassparticle[9]));
713
714     for (Int_t kk=0;kk<ntrkinset;kk++){
715       Int_t idxextimePi = ((6*i+kk)*NIDX)+DELTAIDXEXTIMEPI; 
716       Int_t idxextimeKa = ((6*i+kk)*NIDX)+DELTAIDXEXTIMEKA; 
717       Int_t idxextimePr = ((6*i+kk)*NIDX)+DELTAIDXEXTIMEPR; 
718       // storing in third slot of smallarray the assigned expected time
719       fhPID->Fill(fassparticle[kk]);
720       //      fassparticle[kk]=0;  //assuming all particles are pions
721       if (fassparticle[kk]==0){       //assigned to be a Pi
722         smallarray[idxextimeKa]=-1;
723         smallarray[idxextimePr]=-1;
724       }
725       else if (fassparticle[kk]==1){  //assigned to be a Ka
726         smallarray[idxextimePi]=smallarray[idxextimeKa];
727         smallarray[idxextimeKa]=-1;
728         smallarray[idxextimePr]=-1;
729       }
730       else if (fassparticle[kk]==2){  //assigned to be a Pr
731         smallarray[idxextimePi]=smallarray[idxextimePr];
732         smallarray[idxextimeKa]=-1;
733         smallarray[idxextimePr]=-1;
734       }
735     }
736     delete[] exptof; 
737     delete[] timeofflight;
738     delete[] texp;
739     delete[] index;
740     
741   }
742   return 1;
743 }
744 //---------------------------------------------------------------------
745 Float_t  AliTOFCalibTask::LoopCombPID(Int_t itrkinset, Int_t ntrkinset, Float_t **exptof, Float_t *texp, Float_t *timeofflight, Int_t *index, Float_t chisquarebest){
746
747         // performing combinatorial PID in recursive way
748
749   Int_t indextr = ntrkinset-itrkinset;
750  
751   for (index[indextr]=0;index[indextr]<3;index[indextr]++){
752    Int_t ii = index[indextr];
753    texp[indextr]=exptof[indextr][ii];
754     if (indextr<ntrkinset-1){
755       chisquarebest = LoopCombPID(ntrkinset-indextr-1,ntrkinset,exptof,&texp[0],&timeofflight[0],&index[0], chisquarebest);
756     }
757     
758     else {
759       Float_t *besttimezero = new Float_t[ntrkinset];
760       Float_t *bestchisquare = new Float_t[ntrkinset];
761       Float_t *bestweightedtimezero = new Float_t[ntrkinset];
762       Float_t sumAllweights=0.;
763       Float_t meantzero=0.;
764       Float_t *weightedtimezero = new Float_t[ntrkinset];
765       Float_t *timezero = new Float_t[ntrkinset];
766       
767       AliDebug(2,Form("Configuration = %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, so far chisquarebest = %f ",index[0],index[1],index[2],index[3],index[4],index[5],index[6],index[7],index[8],index[9],chisquarebest)); 
768       for (Int_t itz=0; itz<ntrkinset;itz++) {
769         timezero[itz]=texp[itz]-timeofflight[itz];                  
770         weightedtimezero[itz]=timezero[itz]/TRACKERROR;
771         sumAllweights+=1./TRACKERROR;
772         meantzero+=weightedtimezero[itz];
773       } // end loop for (Int_t itz=0; itz<ntrkinset;itz++)
774       
775       meantzero=meantzero/sumAllweights; // it is given in [ns]
776       
777       // calculate chisquare
778       
779       Float_t chisquare=0.;             
780       for (Int_t icsq=0; icsq<ntrkinset;icsq++) {
781         chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/TRACKERROR;
782       } // end loop for (Int_t icsq=0; icsq<ntrkinset;icsq++) 
783       
784       Int_t npion=0;
785       for (Int_t j=0;j<ntrkinset;j++){
786         if(index[j]==0)npion++;
787       }
788       
789       if(chisquare<=chisquarebest  && ((Float_t) npion/ ((Float_t) ntrkinset)>0.3)){
790         for(Int_t iqsq = 0; iqsq<ntrkinset; iqsq++) {
791           besttimezero[iqsq]=timezero[iqsq]; 
792           bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
793           bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/TRACKERROR; 
794         }
795         
796         for (Int_t j=0;j<ntrkinset;j++){
797           fassparticle[j]=index[j];
798         }                 
799         
800         chisquarebest=chisquare;
801       }
802       delete[] besttimezero;
803       delete[] bestchisquare;
804       delete[] bestweightedtimezero;
805       delete[] weightedtimezero;
806       delete[] timezero;
807     }
808   }
809   return chisquarebest;
810 }