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