1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 Revision 1.6 2007/10/05 10:44:11 zampolli
18 Oversight for debugging purpose fixed
20 Revision 1.5 2007/10/05 10:38:10 zampolli
23 Revision 1.4 2007/10/04 15:36:44 zampolli
24 Updates to new TOF offline calibration schema
26 Revision 1.3 2007/07/31 07:26:16 zampolli
27 Bug fixed in the $ field
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
50 #include <TDirectory.h>
53 #include "AliTOFCalibTask.h"
54 #include "AliESDEvent.h"
55 #include "AliESDtrack.h"
57 #include "AliCDBManager.h"
58 #include "AliCDBMetaData.h"
60 #include "AliCDBStorage.h"
62 //______________________________________________________________________________
63 AliTOFCalibTask::AliTOFCalibTask(const char *name) :
64 AliAnalysisTask(name,""),
94 // Input slot #0 works with an Ntuple
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];
103 for (Int_t i=0;i<11;i++){
107 for (Int_t i =0;i<TOFCHANNELS;i++){
108 fbigarray[i] = new Float_t[CHENTRIES];
110 for (Int_t j =0;j<CHENTRIES;j++){
116 //______________________________________________________________________________
117 AliTOFCalibTask::AliTOFCalibTask(const AliTOFCalibTask &calibtask) :
118 AliAnalysisTask("AliTOFCalibTask",""),
149 fChain=calibtask.fChain;
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;
172 fOutputContainer=calibtask.fOutputContainer;
174 fbigarray = new Float_t*[TOFCHANNELS];
175 findexarray = new Int_t[TOFCHANNELS];
177 for (Int_t i=0;i<11;i++){
178 fassparticle[i]=calibtask.fassparticle[i];
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];
190 //______________________________________________________________________________
191 AliTOFCalibTask:: ~AliTOFCalibTask()
195 AliInfo("TOF Calib Task: Deleting");
196 for (Int_t i =0;i<TOFCHANNELS;i++){
197 delete [] fbigarray[i];
200 delete[] findexarray;
201 delete fOutputContainer;
210 //______________________________________________________________________________
211 AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)
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;
239 this->fbigarray = new Float_t*[TOFCHANNELS];
240 this->findexarray = new Int_t[TOFCHANNELS];
242 for (Int_t i=0;i<11;i++){
243 this->fassparticle[i]=calibtask.fassparticle[i];
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];
254 //--------------------------------------------------------------------------
255 void AliTOFCalibTask::BookHistos(){
259 AliInfo(Form("*** Booking Histograms %s", GetName())) ;
262 new TH1F("hToT", " ToT distribution (ns) ", 400, 0, 40);
264 new TH1F("hTime", " Time distribution (ns) ", 400, 0, 40);
266 new TH1F("hExpTimePi", " Exp Time distribution, Pi (ns) ", 400, 0, 40);
268 new TH1F("hExpTimeKa", " Exp Time distribution, Ka (ns) ", 400, 0, 40);
270 new TH1F("hExpTimePr", " Exp Time distribution, Pr (ns) ", 400, 0, 40);
272 new TH1I("hPID", " Combinatorial PID identity ", 3, 0, 3);
274 new TH1D("hch", " TOF channel ", TOFCHANNELS, 0, TOFCHANNELS);
276 //create the putput container
277 fOutputContainer = new TObjArray(7) ;
278 fOutputContainer->SetName(GetName()) ;
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) ;
289 //----------------------------------------------------------------------------
290 void AliTOFCalibTask::DrawHistos(){
292 // drawing output histos
294 AliInfo(Form("*** Drawing Histograms %s", GetName())) ;
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)");
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");
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");
324 canvasExpTime->Print("ExpTime.gif");
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");
333 canvasPID->Print("PID.gif");
335 TCanvas * canvasrndch = new TCanvas("canvasrndch", " TOF channel ",400, 30, 550, 400);
336 fhch->GetXaxis()->SetTitle("TOF ch");
338 Float_t meanTOFch = 0;
339 for (Int_t ibin=0;ibin<TOFCHANNELS;ibin++){
340 meanTOFch+=(Float_t)fhch->GetBinContent(ibin+1);
343 meanTOFch/=TOFCHANNELS;
344 AliDebug(1,Form(" Mean number of tracks/channel = %f ",meanTOFch));
346 canvasrndch->Print("rndch.gif");
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())) ;
356 //______________________________________________________________________________
357 void AliTOFCalibTask::ConnectInputData(const Option_t*)
359 // Initialization of branch container and histograms
361 // AliLog::SetClassDebugLevel("AliTOFCalibTask",1);
362 AliInfo(Form("*** Initialization of %s", GetName())) ;
365 fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
367 AliError(Form("Input 0 for %s not found\n", GetName()));
371 // One should first check if the branch address was taken by some other task
372 char ** address = (char **)GetBranchAddress(0, "ESD");
374 fESD = (AliESDEvent*)(*address);
376 fESD = new AliESDEvent();
378 fESD->ReadFromTree(fChain) ;
383 //-----------------------------------------------------------------------
384 Bool_t AliTOFCalibTask::Notify()
386 // Initialisation of branch container and histograms
388 AliInfo(Form("*** We are in %s::Notify()", GetName())) ;
391 fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
393 AliError(Form("Input 0 for %s not found\n", GetName()));
397 char ** address = (char **)GetBranchAddress(0, "ESD");
399 fESD = (AliESDEvent*)(*address);
401 fESD = new AliESDEvent();
403 fESD->ReadFromTree(fChain) ;
408 //________________________________________________________________________
409 void AliTOFCalibTask::CreateOutputObjects()
414 //______________________________________________________________________________
415 void AliTOFCalibTask::Exec(Option_t * opt)
420 AliInfo(Form("*** Executing %s", GetName())) ;
422 //******* The loop over events -----------------------------------------------
424 // Processing of one event
425 Long64_t entry = fChain->GetReadEntry() ;
427 AliError("fESD is not connected to the input!") ;
431 if ( !((entry-1)%100) )
432 AliDebug(1,Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
434 // ************************ TOF *************************************
436 fMinTime=22E3; //ns; not used
437 Int_t ntrk = fESD->GetNumberOfTracks() ;
443 AliESDtrack * t = fESD->GetTrack(ntrk) ;
444 //selecting only good quality tracks
445 if (!Select(t)) continue;
447 Int_t ich = Int_t(t->GetTOFCalChannel());
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));
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));
478 fnESDselected+=nselected;
480 PostData(0, fOutputContainer);
483 //_____________________________________________________________________________
484 void AliTOFCalibTask::Terminate(Option_t *)
486 // Processing when the event loop is ended
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;
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));
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));
515 if (!CombPID(&fbigarray[i][0], size)) AliError("ERROR!!!!ERROR!!!");
521 // saving data in a tree
523 AliInfo("Building tree for Calibration");
524 TTree * tree = new TTree("T", "Tree for TOF calibration");
525 Float_t p[CHENTRIESSMALL];
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
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];
543 AliInfo("Putting tree for calibration in Reference data");
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");
554 //_____________________________________________________________________________
556 Bool_t AliTOFCalibTask::Select(AliESDtrack *t){
558 // selecting good quality tracks
559 //track selection: reconstrution to TOF:
560 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
564 //IsStartedTimeIntegral
565 if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
569 if (t->GetStatus() & AliESDtrack::kTRDbackup) {
570 Float_t xout= t->GetOuterParam()->GetX();
571 if (xout<364.25 && xout > 300.) return 0;
574 Double_t time=t->GetTOFsignal();
575 time*=1.E-3; // tof given in nanoseconds
576 if(time >= fMinTime){
581 Double_t mom=t->GetP();
582 if (!(mom<=UPPERMOMBOUND && mom>=LOWERMOMBOUND)){
583 // return 0; // skipping momentum cut
586 Int_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
587 if(assignedTOFcluster==-1){ // not matched
593 //_____________________________________________________________________________
595 Int_t AliTOFCalibTask::SelectOnTime(Float_t *charray, Int_t ntracks, Int_t ich){
597 // discarding tracks with time-mintime < MINTIME
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;
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) {
612 AliDebug(1,Form("Deleting %i track from channel %i, time = %f",ndeleted, ich, time));
614 for (Int_t j=itr+1;j<ntracks;j++){
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];
637 //_____________________________________________________________________________
639 Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
641 // track Combinatorial PID for calibration, only when more than 2 particles
647 AliDebug(1,"Number of tracks in channel smaller than 2, identifying every particle as if it was a pion!");
651 Int_t ntracksinchannel = size/NIDX;
653 Int_t nset = ntracksinchannel/6;
657 if (ntracksinchannel < 6) {
658 AliDebug(2,"Number of tracks in set smaller than 6, Combinatorial PID still applied to this set.");
662 AliInfo(Form("number of sets = %i", nset));
664 for (Int_t i=0; i< nset; i++) {
665 if (i<nset-1)ntrkinset=6;
667 if (ntracksinchannel<6){
668 ntrkinset=ntracksinchannel;
671 ntrkinset=6+Int_t((ntracksinchannel)%6);
674 AliInfo(Form("set = %i, number of tracks in set = %i", i,ntrkinset));
676 for (Int_t ii=0;ii<ntrkinset;ii++){
680 Float_t* timeofflight;
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];
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]));
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);
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);
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;
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]));
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;
738 else if (fassparticle[kk]==1){ //assigned to be a Ka
739 smallarray[idxextimePi]=smallarray[idxextimeKa];
740 smallarray[idxextimeKa]=-1;
741 smallarray[idxextimePr]=-1;
743 else if (fassparticle[kk]==2){ //assigned to be a Pr
744 smallarray[idxextimePi]=smallarray[idxextimePr];
745 smallarray[idxextimeKa]=-1;
746 smallarray[idxextimePr]=-1;
750 delete[] timeofflight;
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){
760 Int_t indextr = ntrkinset-itrkinset;
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);
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];
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++)
786 meantzero=meantzero/sumAllweights; // it is given in [ns]
788 // calculate chisquare
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++)
796 for (Int_t j=0;j<ntrkinset;j++){
797 if(index[j]==0)npion++;
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;
807 for (Int_t j=0;j<ntrkinset;j++){
808 fassparticle[j]=index[j];
811 chisquarebest=chisquare;
813 delete[] besttimezero;
814 delete[] bestchisquare;
815 delete[] bestweightedtimezero;
816 delete[] weightedtimezero;
820 return chisquarebest;