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