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