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 }