git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Oversight fixed
[u/mrichter/AliRoot.git] / TOF / AliTOFCalibTask.cxx
CommitLineData
1ed77580 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/*
45af7330 16$Log$
1e76f1ec 17Revision 1.4 2007/10/04 15:36:44 zampolli
18Updates to new TOF offline calibration schema
19
e26353e8 20Revision 1.3 2007/07/31 07:26:16 zampolli
21Bug fixed in the $ field
22
1ed77580 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>
e26353e8 44#include <TDirectory.h>
45#include <TGrid.h>
1ed77580 46
47#include "AliTOFCalibTask.h"
70ddd0ee 48#include "AliESDEvent.h"
1ed77580 49#include "AliESDtrack.h"
50#include "AliLog.h"
e26353e8 51#include "AliCDBManager.h"
52#include "AliCDBMetaData.h"
53#include "AliCDBId.h"
54#include "AliCDBStorage.h"
55
1ed77580 56//______________________________________________________________________________
57AliTOFCalibTask::AliTOFCalibTask(const char *name) :
58 AliAnalysisTask(name,""),
59 fdir(0),
60 fChain(0),
61 fESD(0),
1ed77580 62 fToT(0),
63 fTime(0),
64 fExpTimePi(0),
65 fExpTimeKa(0),
66 fExpTimePr(0),
1ed77580 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),
e26353e8 84 fOutputContainer(0),
85 frun(0)
1ed77580 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();
70ddd0ee 94 fbigarray = new Float_t*[TOFCHANNELS];
95 findexarray = new Int_t[TOFCHANNELS];
96
e26353e8 97 for (Int_t i=0;i<11;i++){
98 fassparticle[i]=-1;
99 }
100
70ddd0ee 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 }
1ed77580 108}
109
110//______________________________________________________________________________
111AliTOFCalibTask::AliTOFCalibTask(const AliTOFCalibTask &calibtask) :
112 AliAnalysisTask("AliTOFCalibTask",""),
113 fdir(0),
114 fChain(0),
115 fESD(0),
1ed77580 116 fToT(0),
117 fTime(0),
118 fExpTimePi(0),
119 fExpTimeKa(0),
120 fExpTimePr(0),
1ed77580 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),
e26353e8 138 fOutputContainer(0),
139 frun(0)
1ed77580 140{
141 // Copy Constructor.
142 fdir=calibtask.fdir;
143 fChain=calibtask.fChain;
144 fESD=calibtask.fESD;
1ed77580 145 fToT=calibtask.fToT;
146 fTime=calibtask.fTime;
147 fExpTimePi=calibtask.fExpTimePi;
148 fExpTimeKa=calibtask.fExpTimeKa;
149 fExpTimePr=calibtask.fExpTimePr;
1ed77580 150 fMinTime=calibtask.fMinTime;
1ed77580 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;
e26353e8 165 frun=calibtask.frun;
1ed77580 166 fOutputContainer=calibtask.fOutputContainer;
70ddd0ee 167
168 fbigarray = new Float_t*[TOFCHANNELS];
169 findexarray = new Int_t[TOFCHANNELS];
170
e26353e8 171 for (Int_t i=0;i<11;i++){
172 fassparticle[i]=calibtask.fassparticle[i];
173 }
174
70ddd0ee 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
1ed77580 183}
184//______________________________________________________________________________
185AliTOFCalibTask:: ~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//______________________________________________________________________________
202AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)
203{
204 //assignment operator
205 this->fdir=calibtask.fdir;
206 this->fChain=calibtask.fChain;
207 this->fESD=calibtask.fESD;
1ed77580 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;
1ed77580 213 this->fMinTime=calibtask.fMinTime;
1ed77580 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;
70ddd0ee 229
230 this->fbigarray = new Float_t*[TOFCHANNELS];
231 this->findexarray = new Int_t[TOFCHANNELS];
e26353e8 232
233 for (Int_t i=0;i<11;i++){
234 this->fassparticle[i]=calibtask.fassparticle[i];
235 }
70ddd0ee 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 }
1ed77580 243 return *this;
244}
245//--------------------------------------------------------------------------
246void 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//----------------------------------------------------------------------------
281void 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//______________________________________________________________________________
348void 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) {
70ddd0ee 365 fESD = (AliESDEvent*)(*address);
1ed77580 366 } else {
70ddd0ee 367 fESD = new AliESDEvent();
1ed77580 368 }
1e76f1ec 369 fESD->ReadFromTree(fChain) ;
1ed77580 370
371 BookHistos();
372
373}
70ddd0ee 374//-----------------------------------------------------------------------
375Bool_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
1e76f1ec 388 char ** address = (char **)GetBranchAddress(0, "ESD");
389 if (address) {
390 fESD = (AliESDEvent*)(*address);
391 } else {
392 fESD = new AliESDEvent();
393 }
e26353e8 394 fESD->ReadFromTree(fChain) ;
70ddd0ee 395
396 return kTRUE;
397}
398
1ed77580 399//________________________________________________________________________
400void AliTOFCalibTask::CreateOutputObjects()
401{
402// Create histograms
403}
404
405//______________________________________________________________________________
406void 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);
e26353e8 440 ich=3;
1ed77580 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();
e26353e8 444 AliDebug(2,Form(" track # %i in channel %i, time = %f \n",ntrk,ich,time));
1ed77580 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//_____________________________________________________________________________
475void AliTOFCalibTask::Terminate(Option_t *)
476{
477 // Processing when the event loop is ended
478
479 // some plots
480
481 TH1::AddDirectory(0);
e26353e8 482 TDirectory *dir = gDirectory;
1ed77580 483 AliInfo("TOF Calib Task: End of events loop");
e26353e8 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
1ed77580 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;
e26353e8 501 if (i==3) AliInfo(Form(" entries %i in channel %i ",findexarray[i],i));
1ed77580 502 AliDebug(2, Form(" entries %i in channel %i ",findexarray[i],i));
e26353e8 503 if (findexarray[i]<=2) {
1ed77580 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 }
70ddd0ee 509
e26353e8 510
1ed77580 511 DrawHistos();
512
513 // saving data in a tree
514
e26353e8 515 AliInfo("Building tree for Calibration");
516 TTree * tree = new TTree("T", "Tree for TOF calibration");
517 Float_t p[CHENTRIESSMALL];
1ed77580 518 Int_t nentries;
e26353e8 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 }
1ed77580 531 }
e26353e8 532 tree->Fill();
1ed77580 533 }
1ed77580 534
e26353e8 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();
1ed77580 544
1ed77580 545}
1ed77580 546//_____________________________________________________________________________
547
548Bool_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
e26353e8 587Int_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
1ed77580 631Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
632
e26353e8 633 // track Combinatorial PID for calibration, only when more than 2 particles
634 // fall in channel
1ed77580 635
1ed77580 636 Float_t t0offset=0;
e26353e8 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!");
1ed77580 640 return 0;
641 }
e26353e8 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 }
1ed77580 652 }
653
e26353e8 654 AliInfo(Form("number of sets = %i", nset));
655 // loop on sets
1ed77580 656 for (Int_t i=0; i< nset; i++) {
e26353e8 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));
1ed77580 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 }
e26353e8 700
701
1ed77580 702 Float_t chisquarebest=999.;
e26353e8 703 AliInfo(Form(" Set = %i with %i tracks ", i,ntrkinset));
704 chisquarebest = LoopCombPID(ntrkinset, ntrkinset,exptof,&texp[0],&timeofflight[0], &index[0],chisquarebest);
705
1ed77580 706 Float_t confLevel=999;
707 if(chisquarebest<999.){
708 Double_t dblechisquare=(Double_t)chisquarebest;
e26353e8 709 confLevel=(Float_t)TMath::Prob(dblechisquare,ntrkinset-1);
1ed77580 710 }
e26353e8 711
1ed77580 712 // assume they are all pions for fake sets
713 if(confLevel<0.01 || confLevel==999. ){
e26353e8 714 for (Int_t itrk=0; itrk<ntrkinset; itrk++)fassparticle[itrk]=0;
1ed77580 715 }
e26353e8 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]));
1ed77580 718
e26353e8 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;
1ed77580 723 // storing in third slot of smallarray the assigned expected time
e26353e8 724 fhPID->Fill(fassparticle[kk]);
725 // fassparticle[kk]=0; //assuming all particles are pions
726 if (fassparticle[kk]==0){ //assigned to be a Pi
1ed77580 727 smallarray[idxextimeKa]=-1;
728 smallarray[idxextimePr]=-1;
729 }
e26353e8 730 else if (fassparticle[kk]==1){ //assigned to be a Ka
1ed77580 731 smallarray[idxextimePi]=smallarray[idxextimeKa];
732 smallarray[idxextimeKa]=-1;
733 smallarray[idxextimePr]=-1;
734 }
e26353e8 735 else if (fassparticle[kk]==2){ //assigned to be a Pr
1ed77580 736 smallarray[idxextimePi]=smallarray[idxextimePr];
737 smallarray[idxextimeKa]=-1;
738 smallarray[idxextimePr]=-1;
739 }
740 }
e26353e8 741
1ed77580 742 }
743 return 1;
744}
e26353e8 745//---------------------------------------------------------------------
746Float_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);
1ed77580 755 }
e26353e8 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++;
1ed77580 786 }
e26353e8 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;
1ed77580 800 }
801 }
802 }
e26353e8 803 return chisquarebest;
1ed77580 804}