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 }