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