]> git.uio.no Git - u/mrichter/AliRoot.git/blob - TOF/AliTOFCalibTask.cxx
5d3435fabf9f236d76a84b66cd8c6d370c76f941
[u/mrichter/AliRoot.git] / 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 */
18
19 // task to perform TOF calibration
20 // optimized for pp runs
21 // expect a maximum of 100 entries per channel
22 // different ways to calibrate are implemented
23 // C. Zampolli 
24
25 #include <TChain.h>
26 #include <TObject.h> 
27 #include <TCanvas.h>
28 #include <TROOT.h>
29 #include <TSystem.h>
30 #include <TProfile.h>
31 #include <TF1.h>
32 #include <TTree.h>
33 #include <TFile.h>
34 #include <TH1F.h>
35 #include <TH1I.h>
36 #include <TH1D.h>
37 #include <TH2F.h>
38   
39 #include "AliTOFCalibTask.h" 
40 #include "AliTOFChannelTask.h" 
41 #include "AliESDEvent.h" 
42 #include "AliESDtrack.h" 
43 #include "AliLog.h"
44 //______________________________________________________________________________
45 AliTOFCalibTask::AliTOFCalibTask(const char *name) : 
46   AliAnalysisTask(name,""),  
47   fdir(0),
48   fChain(0),
49   fESD(0),
50   fMinEntries(0),
51   fIch(0),
52   fToT(0),
53   fTime(0),
54   fExpTimePi(0),
55   fExpTimeKa(0),
56   fExpTimePr(0),
57   ftree(0x0),
58   fMinTime(0),
59   fbigarray(0x0),
60   findexarray(0x0),
61   fnESD(0),
62   fnESDselected(0),
63   fnESDkTOFout(0),
64   fnESDkTIME(0),
65   fnESDassTOFcl(0),
66   fnESDTIMEcut(0),
67   fnESDTRDcut(0),
68   fhToT(0),
69   fhTime(0),
70   fhExpTimePi(0),
71   fhExpTimeKa(0),
72   fhExpTimePr(0),
73   fhPID(0),
74   fhch(0),
75   fOutputContainer(0)
76 {
77   // Constructor.
78   // Input slot #0 works with an Ntuple
79
80   DefineInput(0, TChain::Class());
81   // Output slot #0 writes into a TH1 container
82   DefineOutput(0,  TObjArray::Class()) ; 
83   fdir = gDirectory->GetPath();
84   fbigarray = new Float_t*[TOFCHANNELS];
85   findexarray = new Int_t[TOFCHANNELS];
86
87   for (Int_t i =0;i<TOFCHANNELS;i++){
88     fbigarray[i] = new Float_t[CHENTRIES];
89     findexarray[i]=0;
90     for (Int_t j =0;j<CHENTRIES;j++){
91       fbigarray[i][j]=-1;
92     }
93   }
94 }
95
96 //______________________________________________________________________________
97 AliTOFCalibTask::AliTOFCalibTask(const AliTOFCalibTask &calibtask) : 
98   AliAnalysisTask("AliTOFCalibTask",""),  
99   fdir(0),
100   fChain(0),
101   fESD(0),
102   fMinEntries(0),
103   fIch(0),
104   fToT(0),
105   fTime(0),
106   fExpTimePi(0),
107   fExpTimeKa(0),
108   fExpTimePr(0),
109   ftree(0x0),
110   fMinTime(0),
111   fbigarray(0x0),
112   findexarray(0x0),
113   fnESD(0),
114   fnESDselected(0),
115   fnESDkTOFout(0),
116   fnESDkTIME(0),
117   fnESDassTOFcl(0),
118   fnESDTIMEcut(0),
119   fnESDTRDcut(0),
120   fhToT(0),
121   fhTime(0),
122   fhExpTimePi(0),
123   fhExpTimeKa(0),
124   fhExpTimePr(0),
125   fhPID(0),
126   fhch(0),
127   fOutputContainer(0)
128 {
129   // Copy Constructor.
130   fdir=calibtask.fdir;
131   fChain=calibtask.fChain;
132   fESD=calibtask.fESD;
133   fMinEntries=calibtask.fMinEntries;
134   fIch=calibtask.fIch;
135   fToT=calibtask.fToT;
136   fTime=calibtask.fTime;
137   fExpTimePi=calibtask.fExpTimePi;
138   fExpTimeKa=calibtask.fExpTimeKa;
139   fExpTimePr=calibtask.fExpTimePr;
140   ftree=calibtask.ftree;
141   fMinTime=calibtask.fMinTime;
142   fnESD=calibtask.fnESD;
143   fnESDselected=calibtask.fnESDselected;
144   fnESDkTOFout=calibtask.fnESDkTOFout;
145   fnESDkTIME=calibtask.fnESDkTIME;
146   fnESDassTOFcl=calibtask.fnESDassTOFcl;
147   fnESDTIMEcut=calibtask.fnESDTIMEcut;
148   fnESDTRDcut=calibtask.fnESDTRDcut;
149   fhToT=calibtask.fhToT;
150   fhTime=calibtask.fhTime;
151   fhExpTimePi=calibtask.fhExpTimePi;
152   fhExpTimeKa=calibtask.fhExpTimeKa;
153   fhExpTimePr=calibtask.fhExpTimePr;
154   fhPID=calibtask.fhPID;
155   fhch=calibtask.fhch;
156   fOutputContainer=calibtask.fOutputContainer; 
157
158   fbigarray = new Float_t*[TOFCHANNELS];
159   findexarray = new Int_t[TOFCHANNELS];
160
161   for (Int_t i =0;i<TOFCHANNELS;i++){
162     fbigarray[i] = new Float_t[CHENTRIES];
163     findexarray[i]=calibtask.findexarray[i];
164     for (Int_t j =0;j<CHENTRIES;j++){
165       fbigarray[i][j]=calibtask.fbigarray[i][j];
166     }
167   }
168
169 }
170 //______________________________________________________________________________
171 AliTOFCalibTask:: ~AliTOFCalibTask() 
172 {
173   // destructor
174
175   AliInfo("TOF Calib Task: Deleting");
176   delete[] fbigarray;
177   delete findexarray;
178   delete fOutputContainer;
179   delete fhToT;
180   delete fhTime;
181   delete fhExpTimePi; 
182   delete fhExpTimeKa; 
183   delete fhExpTimePr; 
184   delete fhPID; 
185   delete fhch;
186 }
187 //______________________________________________________________________________
188 AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)  
189
190    //assignment operator
191   this->fdir=calibtask.fdir;
192   this->fChain=calibtask.fChain;
193   this->fESD=calibtask.fESD;
194   this->fMinEntries=calibtask.fMinEntries;
195   this->fIch=calibtask.fIch;
196   this->fToT=calibtask.fToT;
197   this->fTime=calibtask.fTime;
198   this->fExpTimePi=calibtask.fExpTimePi;
199   this->fExpTimeKa=calibtask.fExpTimeKa;
200   this->fExpTimePr=calibtask.fExpTimePr;
201   this->ftree=calibtask.ftree;
202   this->fMinTime=calibtask.fMinTime;
203   this->fnESD=calibtask.fnESD;
204   this->fnESDselected=calibtask.fnESDselected;
205   this->fnESDkTOFout=calibtask.fnESDkTOFout;
206   this->fnESDkTIME=calibtask.fnESDkTIME;
207   this->fnESDassTOFcl=calibtask.fnESDassTOFcl;
208   this->fnESDTIMEcut=calibtask.fnESDTIMEcut;
209   this->fnESDTRDcut=calibtask.fnESDTRDcut;
210   this->fhToT=calibtask.fhToT;
211   this->fhTime=calibtask.fhTime;
212   this->fhExpTimePi=calibtask.fhExpTimePi;
213   this->fhExpTimeKa=calibtask.fhExpTimeKa;
214   this->fhExpTimePr=calibtask.fhExpTimePr;
215   this->fOutputContainer=calibtask.fOutputContainer; 
216   this->fhPID=calibtask.fhPID;
217   this->fhch=calibtask.fhch;
218
219   this->fbigarray = new Float_t*[TOFCHANNELS];
220   this->findexarray = new Int_t[TOFCHANNELS];
221   for (Int_t i =0;i<TOFCHANNELS;i++){
222     this->fbigarray[i] = new Float_t[CHENTRIES];
223     this->findexarray[i]=calibtask.findexarray[i];
224     for (Int_t j =0;j<CHENTRIES;j++){
225       this->fbigarray[i][j]=calibtask.fbigarray[i][j];
226     }
227   }
228   return *this;
229 }
230 //--------------------------------------------------------------------------
231 void AliTOFCalibTask::BookHistos(){
232
233   //booking histos
234
235   AliInfo(Form("*** Booking Histograms %s", GetName())) ; 
236
237   fhToT=
238     new TH1F("hToT", " ToT distribution (ns) ", 400, 0, 40);
239   fhTime=
240     new TH1F("hTime", " Time distribution (ns) ", 400, 0, 40);
241   fhExpTimePi=
242     new TH1F("hExpTimePi", " Exp Time distribution, Pi (ns) ", 400, 0, 40);
243   fhExpTimeKa=
244     new TH1F("hExpTimeKa", " Exp Time distribution, Ka (ns) ", 400, 0, 40);
245   fhExpTimePr=
246     new TH1F("hExpTimePr", " Exp Time distribution, Pr (ns) ", 400, 0, 40);
247   fhPID=
248     new TH1I("hPID", " Combinatorial PID identity ", 3, 0, 3);
249   fhch=
250     new TH1D("hch", " TOF channel ", TOFCHANNELS, 0, TOFCHANNELS);
251
252   //create the putput container
253   fOutputContainer = new TObjArray(7) ; 
254   fOutputContainer->SetName(GetName()) ; 
255
256   fOutputContainer->AddAt(fhToT,             0) ; 
257   fOutputContainer->AddAt(fhTime,            1) ; 
258   fOutputContainer->AddAt(fhExpTimePi,       2) ; 
259   fOutputContainer->AddAt(fhExpTimeKa,       3) ; 
260   fOutputContainer->AddAt(fhExpTimePr,       4) ; 
261   fOutputContainer->AddAt(fhPID,             5) ; 
262   fOutputContainer->AddAt(fhch,              6) ; 
263
264 }
265 //----------------------------------------------------------------------------
266 void AliTOFCalibTask::DrawHistos(){
267
268   // drawing output histos
269
270   AliInfo(Form("*** Drawing Histograms %s", GetName())) ; 
271
272   TCanvas * canvasToTTime = new TCanvas("canvasToTTime", " ToT and Time ",400, 30, 550, 630) ;
273   canvasToTTime->Divide(1,2);
274   canvasToTTime->cd(1);
275   fhToT->SetLineColor(4);
276   fhToT->GetXaxis()->SetTitle("ToT (ns)");
277   fhToT->Draw("hist");
278   canvasToTTime->cd(2);
279   fhTime->SetLineColor(4);
280   fhTime->GetXaxis()->SetTitle("Time (ns)");
281   fhTime->Draw("hist");
282   canvasToTTime->Update();
283   canvasToTTime->Print("ToTTime.gif");
284
285   TCanvas * canvasExpTime = new TCanvas("canvasExpTime", " Expected Times ",400, 30, 550, 630) ;
286   canvasExpTime->Divide(1,3);
287   canvasExpTime->cd(1);
288   fhExpTimePi->SetLineColor(4);
289   fhExpTimePi->GetXaxis()->SetTitle("Exp Time (ns), #pi");
290   fhExpTimePi->Draw("hist");
291   canvasExpTime->cd(2);
292   fhExpTimeKa->SetLineColor(4);
293   fhExpTimeKa->GetXaxis()->SetTitle("Exp Time (ns), K");
294   fhExpTimeKa->Draw("hist");
295   canvasExpTime->cd(3);
296   fhExpTimePr->SetLineColor(4);
297   fhExpTimePr->GetXaxis()->SetTitle("Exp Time (ns), p");
298   fhExpTimePr->Draw("hist");
299
300   canvasExpTime->Print("ExpTime.gif");
301
302   TCanvas * canvasPID = new TCanvas("canvasPID", " Combinatorial PID ",400, 30, 550, 400);
303   fhPID->GetXaxis()->SetTitle("Comb PID");
304   fhPID->GetXaxis()->SetBinLabel(1,"#pi");
305   fhPID->GetXaxis()->SetBinLabel(2,"K");
306   fhPID->GetXaxis()->SetBinLabel(3,"p");
307   fhPID->Draw("hist");
308
309   canvasPID->Print("PID.gif");
310
311   TCanvas * canvasrndch = new TCanvas("canvasrndch", " TOF channel ",400, 30, 550, 400);
312   fhch->GetXaxis()->SetTitle("TOF ch");
313   fhch->Draw("hist");
314   Float_t meanTOFch = 0;
315   for (Int_t ibin=0;ibin<TOFCHANNELS;ibin++){
316     meanTOFch+=(Float_t)fhch->GetBinContent(ibin+1);
317   }
318
319   meanTOFch/=TOFCHANNELS;
320   AliDebug(1,Form(" Mean number of tracks/channel = %f ",meanTOFch));
321
322   canvasrndch->Print("rndch.gif");
323
324   char line[1024] ; 
325   sprintf(line, ".!tar -zcvf %s.tar.gz *.gif", GetName()) ; 
326   gROOT->ProcessLine(line);
327   sprintf(line, ".!rm -fR *.gif"); 
328   gROOT->ProcessLine(line);
329   AliInfo(Form("*** TOF Calib Task: plots saved in %s.tar.gz...\n", GetName())) ;
330 }
331
332 //______________________________________________________________________________
333 void AliTOFCalibTask::ConnectInputData(const Option_t*)
334 {
335   // Initialization of branch container and histograms 
336     
337   //  AliLog::SetClassDebugLevel("AliTOFCalibTask",1);
338   AliInfo(Form("*** Initialization of %s", GetName())) ; 
339   
340   // Get input data
341   fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
342   if (!fChain) {
343     AliError(Form("Input 0 for %s not found\n", GetName()));
344     return ;
345   }
346   
347   // One should first check if the branch address was taken by some other task
348   char ** address = (char **)GetBranchAddress(0, "ESD");
349   if (address) {
350     fESD = (AliESDEvent*)(*address);
351   } else {
352     fESD = new AliESDEvent();
353     AliInfo(" qui ok ");
354     //    fESD = (AliESDEvent*)fChain->GetTree()->GetUserInfo()->FindObject("AliESDEvent");
355     fESD->ReadFromTree(fChain) ;  
356   }
357
358   BookHistos();
359
360 }
361 //-----------------------------------------------------------------------
362 Bool_t AliTOFCalibTask::Notify()
363 {
364   // Initialisation of branch container and histograms 
365     
366   AliInfo(Form("*** We are in  %s::Notify()", GetName())) ; 
367   
368   // Get input data
369   fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
370   if (!fChain) {
371     AliError(Form("Input 0 for %s not found\n", GetName()));
372     return kFALSE;
373   }
374   
375   // One should first check if the branch address was taken by some other task
376   char ** address = (char **)GetBranchAddress(0, "ESD") ;
377   if (address) 
378     fESD = (AliESDEvent *)(*address) ; 
379   else {
380     fESD = new AliESDEvent();
381     //    fESD = (AliESDEvent*)fChain->GetTree()->GetUserInfo()->FindObject("AliESDEvent");
382     fESD->ReadFromTree(fChain) ;  
383     //SetBranchAddress(0,"ESD",&fESD);
384   }
385
386   return kTRUE;
387 }
388
389 //________________________________________________________________________
390 void AliTOFCalibTask::CreateOutputObjects()
391 {
392 // Create histograms
393 }
394
395 //______________________________________________________________________________
396 void AliTOFCalibTask::Exec(Option_t * opt) 
397 {
398
399   // main 
400
401   AliInfo(Form("*** Executing %s", GetName())) ; 
402
403 //******* The loop over events -----------------------------------------------
404
405   // Processing of one event
406   Long64_t entry = fChain->GetReadEntry() ;  
407   if (!fESD) {
408     AliError("fESD is not connected to the input!") ; 
409     return ; 
410   }
411
412   if ( !((entry-1)%100) ) 
413     AliDebug(1,Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ; 
414   
415   // ************************  TOF *************************************
416
417   fMinTime=22E3;   //ns; not used
418   Int_t ntrk = fESD->GetNumberOfTracks() ;
419   fnESD+=ntrk;
420   Int_t nselected = 0;
421   Int_t itr = -1;
422   while ( ntrk-- ) {
423     itr++;
424     AliESDtrack * t = fESD->GetTrack(ntrk) ;
425     //selecting only good quality tracks
426     if (!Select(t)) continue;
427     nselected++;
428     Int_t ich = Int_t(t->GetTOFCalChannel()); 
429     fhch->Fill(ich);
430     AliDebug(2,Form(" ESD in channel %i, filling array %i",t->GetTOFCalChannel(),ich));
431     Float_t tot = t->GetTOFsignalToT();
432     Float_t time = t->GetTOFsignalRaw();
433     AliDebug(2,Form(" track # %i in schannel %i, time = %f \n",ntrk,ich,time)); 
434     Double_t expTime[10]; 
435     t->GetIntegratedTimes(expTime);
436     Float_t expTimePi = expTime[2]*1.E-3;
437     Float_t expTimeKa = expTime[3]*1.E-3;
438     Float_t expTimePr = expTime[4]*1.E-3;
439     if (findexarray[ich]==(Int_t)(CHENTRIES/NIDX)) {
440       AliInfo(Form("too many tracks in channel %i, not storing any more...",ich));
441       continue;
442     } 
443     findexarray[ich]++;
444     AliDebug(2,Form("tracks in channel %i = %i, storing... ",ich, findexarray[ich] ));
445     Int_t ientry=(findexarray[ich]-1)*NIDX;
446     fbigarray[ich][ientry+DELTAIDXTOT]=tot;  //in ns
447     fbigarray[ich][ientry+DELTAIDXTIME]=time*1E-3; // in ns
448     fbigarray[ich][ientry+DELTAIDXEXTIMEPI]=expTimePi;
449     fbigarray[ich][ientry+DELTAIDXEXTIMEKA]=expTimeKa;
450     fbigarray[ich][ientry+DELTAIDXEXTIMEPR]=expTimePr;
451     fhToT->Fill(fbigarray[ich][ientry+DELTAIDXTOT]);
452     fhTime->Fill(fbigarray[ich][ientry+DELTAIDXTIME]);
453     fhExpTimePi->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEPI]);
454     fhExpTimeKa->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEKA]);
455     fhExpTimePr->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEPR]);
456     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));
457   }
458   fnESDselected+=nselected;
459
460   PostData(0, fOutputContainer);  
461 }
462
463 //_____________________________________________________________________________
464 void AliTOFCalibTask::Terminate(Option_t *)
465 {
466   // Processing when the event loop is ended
467   
468   // some plots
469
470   TH1::AddDirectory(0);
471   AliInfo("TOF Calib Task: End of events loop");
472   AliInfo(Form(" Number of analyzed ESD tracks: %i\n",fnESD));
473   AliInfo(Form(" Number of selected ESD tracks: %i\n",fnESDselected));
474   AliInfo(Form(" Number of ESD tracks with kTOFout: %i\n",fnESDkTOFout));
475   AliInfo(Form(" Number of ESD tracks with kTIME: %i\n",fnESDkTIME));
476   AliInfo(Form(" Number of ESD tracks with TRDcut: %i\n",fnESDTRDcut));
477   AliInfo(Form(" Number of ESD tracks with TIMEcut: %i\n",fnESDTIMEcut));
478   AliInfo(Form(" Number of ESD tracks with assTOFcl: %i\n",fnESDassTOFcl));
479
480   for (Int_t i = 0;i<TOFCHANNELS;i++){
481     Int_t size=findexarray[i]*NIDX;
482     AliDebug(2, Form(" entries %i in channel %i ",findexarray[i],i));
483     if (findexarray[i]<6) {
484       AliDebug(1, Form(" not enough statistics for combined PID for channel %i, putting all the tracks as if they were pions",i));
485       continue;
486     }
487     if (!CombPID(&fbigarray[i][0], size)) AliError("ERROR!!!!ERROR!!!");
488   }
489   
490   DrawHistos();
491
492   // saving data in a tree
493   
494   TDirectory *dir = gDirectory;
495   TFile *outFile = 0;
496   TTree * tree = 0x0;
497   Bool_t isThere=kFALSE;
498   const char *dirname = "./";
499   TString filename= "TOFCalib.root";
500
501   if((gSystem->FindFile(dirname,filename))!=NULL){
502     isThere=kTRUE;
503   }
504   
505   if(!isThere){
506     AliInfo("File with tree for Calibration not there, creating it");
507     outFile = new TFile( "TOFCalib.root","RECREATE");
508     outFile->cd();
509     tree = new TTree("T", "Tree for TOF calibration");
510     Float_t p[CHENTRIESSMALL];
511     Int_t nentries;
512     tree->Branch("nentries",&nentries,"nentries/I");
513     tree->Branch("TOFentries",p,"TOFentries[nentries]/F");
514     for (Int_t i=0;i<TOFCHANNELS;i++){
515       nentries=findexarray[i]*(NIDXSMALL); // when filling small array, 
516                                            // only first 3 floats taken 
517                                            // into account
518       for (Int_t j=0; j<findexarray[i];j++){
519         for (Int_t k=0; k<NIDXSMALL;k++){
520           Int_t index1= j*NIDXSMALL+k;   // index in small array
521           Int_t index2=j*NIDX+k;         // index in big array
522           p[index1]=fbigarray[i][index2];
523         }
524       }
525       tree->Fill();
526     }
527     tree->Write("",TObject::kOverwrite);
528     outFile->Close();
529     //delete outFile;
530     //outFile=0x0;
531     //delete tree;
532     //tree=0x0;
533   }
534   else{
535     AliInfo("File with tree for Calibration already there, updating it");
536     outFile = new TFile( "TOFCalib.root","READ");
537     outFile->cd();
538     tree=(TTree*)outFile->Get("T");
539     Float_t p[CHENTRIESSMALL];
540     Float_t ptemp[MAXCHENTRIESSMALL];
541     Int_t nentries;
542     Int_t nentries1;
543     Int_t nentriestemp;
544     tree->SetBranchAddress("nentries",&nentries);
545     tree->SetBranchAddress("TOFentries",p);
546     TFile * outFile1 = new TFile( "TOFCalib.root","RECREATE");
547     outFile1->cd();
548     TTree *treeTemp = new TTree("T","Tree for TOF calibration");
549     treeTemp->Branch("nentries",&nentriestemp,"nentries/I");
550     treeTemp->Branch("TOFentries",ptemp,"TOFentries[nentries]/F");
551     for (Int_t i=0;i<TOFCHANNELS;i++){
552       tree->GetEntry(i);
553       nentries1=nentries;
554       if (nentries+findexarray[i]*(NIDXSMALL) > MAXCHENTRIESSMALL){
555         AliDebug(2, Form("findexarray[%i] = %i, nentries = %i",i,findexarray[i],nentries)); 
556         AliInfo(Form("Too many entries in channel %i, stopping at %i  entries, no storing any more...",i,MAXCHENTRIESSMALL/NIDXSMALL));
557         findexarray[i]=(MAXCHENTRIESSMALL-nentries)/NIDXSMALL;
558       }
559       for (Int_t kk = 0;kk<nentries;kk++){
560         ptemp[kk]=p[kk];
561       }
562       nentries+=findexarray[i]*(NIDXSMALL); // when filling small array, 
563                                             // only first 3 floats taken 
564                                             // into account 
565       nentriestemp=nentries; 
566       for (Int_t j=0; j<findexarray[i];j++){
567         for (Int_t k=0; k<NIDXSMALL;k++){
568           Int_t index1= j*NIDXSMALL+k;   // index in small array
569           Int_t index2=j*NIDX+k;   // index in big array
570           ptemp[index1+nentries1]=fbigarray[i][index2];
571         }
572       }
573       treeTemp->Fill();
574     }
575     treeTemp->SetName("T");
576     treeTemp->SetTitle("Tree for TOF calibration");
577     outFile->Close();
578     treeTemp->Write("",TObject::kOverwrite);
579     outFile1->Close();
580   }
581   dir->cd();
582   Int_t calibrationStatus = 0;
583   //  Int_t ch[2]={3,1003};
584   //  calibrationStatus = Calibrate(2,ch,"save");
585   //calibrationStatus = Calibrate(3,4,"save");
586   //calibrationStatus = CalibrateFromProfile(3,"save");
587   //  calibrationStatus = Calibrate(3,"save");
588   calibrationStatus = Calibrate("");
589   AliInfo(Form("Calibration Status = %i, bye.... \n",calibrationStatus));
590   
591 }
592 //----------------------------------------------------------------------------
593 Int_t AliTOFCalibTask::Calibrate(Int_t ichmin, Int_t ichmax, Option_t *optionSave, Option_t *optionFit){
594
595   // calibrating summing more than one channels
596   // computing calibration parameters
597   // Returning codes:
598   // 0 -> everything was ok
599   // 1 -> no tree for calibration found
600   // 2 -> not enough statistics to perform calibration
601   // 3 -> problems with arrays
602   
603   TH1::AddDirectory(0);
604
605   AliInfo(Form("*** Calibrating Histograms %s, summing more channels, from channel %i, to channel %i, storing Calib Pars in channel %i", GetName(),ichmin,ichmax,ichmin)) ; 
606   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
607   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
608   const char *dirname = "./";
609   TString filename= "TOFCalib.root";
610   
611   if((gSystem->FindFile(dirname,filename))==NULL){
612     AliInfo("No file with tree for calibration found! exiting....");
613     return 1;
614   }
615
616   TFile *file = new TFile("TOFCalib.root","READ");
617   TTree *tree = (TTree*)file->Get("T");
618   Float_t p[MAXCHENTRIESSMALL];
619   Int_t nentries;
620   tree->SetBranchAddress("nentries",&nentries);
621   tree->SetBranchAddress("TOFentries",p);
622
623   Float_t nentriesmean =0;
624   for (Int_t i=ichmin; i<ichmax; i++){
625     tree->GetEntry(i);
626     Int_t ntracks=nentries/3;
627     nentriesmean+=ntracks;
628   }
629
630   //  nentriesmean/=(ichmax-ichmin);
631   if (nentriesmean < MEANENTRIES) {
632     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
633     return 2;
634   }
635
636   //filling ToT and Time arrays
637
638   Int_t nbinToT = 100;  // ToT bin width in Profile = 48.8 ps 
639   Float_t minToT = 0;   // ns
640   Float_t maxToT = 4.88;  // ns
641   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
642
643   TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
644   TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
645   Int_t ntracksTot = 0;
646   Int_t ntracks = 0;
647   Double_t binsProfile[101]; // sized larger than necessary, the correct 
648                              // dim being set in the booking of the profile
649   Int_t nusefulbins=0;
650   Float_t meantime=0;
651   for (Int_t i = ichmin;i<ichmax;i++){
652     tree->GetEntry(i);
653     ntracksTot+=nentries/3;
654     ntracks=nentries/3;
655     AliDebug(2,Form("channel %i, nentries = %i, ntracks = %i",i,nentries, ntracks));
656     for (Int_t j=0;j<ntracks;j++){
657       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
658       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
659       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
660       Float_t tot = p[idxexToT];
661       hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
662       meantime+=p[idxexTime]-p[idxexExTime];
663       hToT->Fill(tot);
664     }
665   }
666   nusefulbins = FindBins(hToT,&binsProfile[0]);
667   meantime/=ntracksTot;
668   AliDebug(2, Form("meantime = %f",meantime));
669   
670   for (Int_t j=1;j<=nusefulbins;j++) {
671     AliDebug(2,Form(" summing channels from %i to %i, nusefulbins = %i, bin %i = %f",ichmin,ichmax,nusefulbins,j,binsProfile[j])); 
672   }
673
674   TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G");  // CHECK THE BUILD OPTION, PLEASE!!!!!!
675   TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
676   for (Int_t i=ichmin; i<ichmax; i++){
677     tree->GetEntry(i);
678     ntracks=nentries/3;
679     for (Int_t j=0;j<ntracks;j++){
680       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
681       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
682       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
683       Float_t tot = p[idxexToT];
684       Float_t time = p[idxexTime]-p[idxexExTime];
685       AliDebug (2, Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
686       hSlewingProf->Fill(tot,time);  // if meantime is not used, the fill may be moved in the loop above
687       htimetot->Fill(tot,time-meantime);  // if meantime is not used, the fill may be moved in the loop above
688     }
689   }
690   hSlewingProf->Fit("pol5",optionFit,"",0,4);
691   TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
692   Float_t par[6];    
693   for(Int_t kk=0;kk<6;kk++){
694     par[kk]=calibfunc->GetParameter(kk);
695     AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
696   }
697
698   if(strstr(optionSave,"save")){
699     TFile * fileProf = new TFile("TOFCalibSave.root","recreate");
700     fileProf->cd(); 
701     TString profName=Form("Profile%06i_%06i",ichmin,ichmax);
702     TString timeTotName=Form("TimeTot%06i_%06i",ichmin,ichmax);
703     TString totName=Form("Tot%06i_%06i",ichmin,ichmax);
704     TString deltaName=Form("Delta%06i_%06i",ichmin,ichmax);
705     hSlewingProf->Write(profName);
706     htimetot->Write(timeTotName);
707     hToT->Write(totName);
708     hdeltaTime->Write(deltaName);
709     fileProf->Close();
710     delete fileProf;
711     fileProf=0x0;
712   }
713
714   delete hToT;
715   hToT=0x0;
716   delete hSlewingProf;
717   hSlewingProf=0x0;
718   delete htimetot;
719   htimetot=0x0;
720   delete hdeltaTime;
721   hdeltaTime=0x0;
722   file->Close();
723   delete file;
724   file=0x0;
725
726   AliTOFChannelTask * calChannel = new AliTOFChannelTask();
727   calChannel->SetSlewPar(par);
728   // saving parameters in chmin
729   arrayCal->AddAt(calChannel,ichmin);
730   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
731   arrayCal->Write("array",TObject::kSingleKey);
732   filecalib->Close();
733   delete filecalib;
734   filecalib = 0x0;
735   delete arrayCal;
736   arrayCal = 0x0;
737   delete calChannel;
738   calChannel =0x0;
739   return 0;
740 }
741 //----------------------------------------------------------------------------
742 Int_t AliTOFCalibTask::Calibrate(Int_t i, Option_t *optionSave, Option_t *optionFit){
743
744   // computing calibration parameters for channel i
745   // Returning codes:
746   // 0 -> everything was ok
747   // 1 -> no tree for calibration found
748   // 2 -> not enough statistics to perform calibration
749   // 3 -> problems with arrays
750
751   TH1::AddDirectory(0);
752   
753   AliInfo(Form("*** Calibrating Histograms (one channel) %s", GetName())) ; 
754   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
755   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
756   const char *dirname = "./";
757   TString filename= "TOFCalib.root";
758   
759   if((gSystem->FindFile(dirname,filename))==NULL){
760     AliInfo("No file with tree for calibration found! exiting....");
761     return 1;
762   }
763
764   TFile *file = new TFile("TOFCalib.root","READ");
765   TTree *tree = (TTree*)file->Get("T");
766   Float_t p[MAXCHENTRIESSMALL];
767   Int_t nentries;
768   tree->SetBranchAddress("nentries",&nentries);
769   tree->SetBranchAddress("TOFentries",p);
770
771   tree->GetEntry(i);
772   Int_t nentriesmean=nentries/3;
773
774   if (nentriesmean < MEANENTRIES) {  
775     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
776     return 2;
777   }
778
779   //filling ToT and Time arrays
780
781   Int_t nbinToT = 100;  // ToT bin width in Profile = 48.8 ps 
782   Float_t minToT = 0;   // ns
783   Float_t maxToT = 4.88;  // ns
784   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
785
786   TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
787   TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
788   Int_t ntracks = 0;
789   Double_t binsProfile[101]; // sized larger than necessary, the correct 
790                              // dim being set in the booking of the profile
791   Int_t nusefulbins=0;
792   Float_t meantime=0;
793   ntracks=nentries/3;
794   AliDebug(2,Form("channel %i, nentries = %i, ntracks = %i",i ,nentries, ntracks));
795   for (Int_t j=0;j<ntracks;j++){
796     Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
797     Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
798     Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
799     Float_t tot = p[idxexToT];
800     meantime+=p[idxexTime]-p[idxexExTime];
801     hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
802     hToT->Fill(tot);
803   }
804   
805   nusefulbins = FindBins(hToT,&binsProfile[0]);
806   meantime/=ntracks;
807   AliDebug(2,Form("meantime = %f",meantime));
808   
809   for (Int_t j=1;j<=nusefulbins;j++) {
810     AliDebug(2,Form(" channel %i, nusefulbins = %i, bin %i = %f",i,nusefulbins,j,binsProfile[j])); 
811   }
812
813   TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G");  // CHECK THE BUILD OPTION, PLEASE!!!!!!
814   TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
815   tree->GetEntry(i);
816   ntracks=nentries/3;
817   for (Int_t j=0;j<ntracks;j++){
818     Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
819     Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
820     Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
821     Float_t tot = p[idxexToT];
822     Float_t time = p[idxexTime]-p[idxexExTime];
823     AliDebug (2,Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
824     hSlewingProf->Fill(tot,time);  // if meantime is not used, the fill may be moved in the loop above
825     htimetot->Fill(tot,time-meantime);  // if meantime is not used, the fill may be moved in the loop above
826   }
827   hSlewingProf->Fit("pol5",optionFit,"",0,4);
828   TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
829   Float_t par[6];    
830   for(Int_t kk=0;kk<6;kk++){
831     par[kk]=calibfunc->GetParameter(kk);
832     AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
833   }
834
835   file->Close();
836   delete file;
837   file=0x0;
838
839   if(strstr(optionSave,"save")){
840     TFile * fileProf = new TFile("TOFCalibSave.root","recreate");
841     fileProf->cd();   
842     TString profName=Form("Profile%06i",i);
843     TString timeTotName=Form("TimeTot%06i",i);
844     TString totName=Form("Tot%06i",i);
845     TString deltaName=Form("Delta%06i",i);
846     hSlewingProf->Write(profName);
847     htimetot->Write(timeTotName);
848     hToT->Write(totName);
849     hdeltaTime->Write(deltaName);
850     fileProf->Close();
851     delete fileProf;
852     fileProf=0x0;
853   }
854
855   delete hToT;
856   hToT=0x0; 
857   delete hSlewingProf;
858   hSlewingProf=0x0;
859   delete htimetot;
860   htimetot=0x0;
861   delete hdeltaTime;
862   hdeltaTime=0x0;
863
864   AliTOFChannelTask * calChannel = new AliTOFChannelTask();
865   calChannel->SetSlewPar(par);
866   arrayCal->AddAt(calChannel,i);
867   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
868   arrayCal->Write("array",TObject::kSingleKey);
869   filecalib->Close();
870   delete filecalib;
871   filecalib=0x0;
872   delete arrayCal;
873   arrayCal = 0x0;
874   delete calChannel;
875   calChannel =0x0;
876   
877   return 0;
878 }
879 //----------------------------------------------------------------------------
880 Int_t AliTOFCalibTask::Calibrate(Int_t nch, Int_t *ch, Option_t *optionSave, Option_t *optionFit){
881
882   // calibrating an array of channels
883   // computing calibration parameters
884   // Returning codes:
885   // 0 -> everything was ok
886   // 1 -> no tree for calibration found
887   // 2 -> not enough statistics to perform calibration
888   // 3 -> problems with arrays
889   
890   TH1::AddDirectory(0);
891
892   AliInfo(Form("*** Calibrating Histograms %s, number of channels = %i", GetName(),nch)) ; 
893   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
894   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
895   for (Int_t ich=0; ich<nch; ich++){
896     Int_t i = ch[ich];
897     AliInfo(Form("Calibrating channel = %i",i )) ; 
898   }
899   const char *dirname = "./";
900   TString filename= "TOFCalib.root";
901   
902   if((gSystem->FindFile(dirname,filename))==NULL){
903     AliInfo("No file with tree for calibration found! exiting....");
904     return 1;
905   }
906
907   TFile *file = new TFile("TOFCalib.root","READ");
908   TTree *tree = (TTree*)file->Get("T");
909   Float_t p[MAXCHENTRIESSMALL];
910   Int_t nentries;
911   tree->SetBranchAddress("nentries",&nentries);
912   tree->SetBranchAddress("TOFentries",p);
913
914   Float_t nentriesmean =0;
915   for (Int_t ich=0; ich<nch; ich++){
916     Int_t i = ch[ich];
917     tree->GetEntry(i);
918     Int_t ntracks=nentries/3;
919     nentriesmean+=ntracks;
920   }
921
922   nentriesmean/=nch;
923   if (nentriesmean < MEANENTRIES) { 
924     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
925     return 2;
926   }
927
928   //filling ToT and Time arrays
929
930   Int_t nbinToT = 100;  // ToT bin width in Profile = 48.8 ps 
931   Float_t minToT = 0;   // ns
932   Float_t maxToT = 4.88;  // ns
933   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
934   arrayCal->SetOwner();
935   TFile * fileProf=0x0;
936   if(strstr(optionSave,"save")){
937     fileProf = new TFile("TOFCalibSave.root","recreate");
938   }
939   for (Int_t ich=0; ich<nch; ich++) {
940     TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
941     TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
942     Double_t binsProfile[101]; // sized larger than necessary, the correct 
943                               // dim being set in the booking of the profile
944     TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
945     Int_t ntracksTot = 0;
946     Int_t ntracks = 0;
947     Int_t nusefulbins=0;
948     Float_t meantime=0;
949     Int_t i = ch[ich];
950     AliDebug(2,Form("Calibrating channel %i",i));
951     tree->GetEntry(i);
952     ntracksTot+=nentries/3;
953     ntracks=nentries/3;
954     if (ntracks < MEANENTRIES) {
955       AliInfo(Form(" Too small mean number of entires in channel %i (number of tracks = %f), not calibrating channel and continuing.....",i,ntracks));
956       continue;
957     }
958     for (Int_t j=0;j<ntracks;j++){
959       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
960       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
961       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
962       Float_t tot = p[idxexToT];
963       hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
964       meantime+=p[idxexTime]-p[idxexExTime];
965       hToT->Fill(tot);
966     }
967
968     nusefulbins = FindBins(hToT,&binsProfile[0]);
969     meantime/=ntracksTot;
970     for (Int_t j=1;j<=nusefulbins;j++) {
971       AliDebug(2,Form(" channel %i, nusefulbins = %i, bin %i = %f",i,nusefulbins,j,binsProfile[j])); 
972     }
973
974     TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G");  // CHECK THE BUILD OPTION, PLEASE!!!!!!
975     tree->GetEntry(i);
976     ntracks=nentries/3;
977     for (Int_t j=0;j<ntracks;j++){
978       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
979       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
980       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
981       Float_t tot = p[idxexToT];
982       Float_t time = p[idxexTime]-p[idxexExTime];
983       AliDebug(2,Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
984       hSlewingProf->Fill(tot,time);  // if meantime is not used, the fill may be moved in the loop above
985       htimetot->Fill(tot,time-meantime);  // if meantime is not used, the fill may be moved in the loop above
986     }
987     
988     hSlewingProf->Fit("pol5",optionFit,"",1,4);
989     TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
990     Float_t par[6];    
991     for(Int_t kk=0;kk<6;kk++){
992       par[kk]=calibfunc->GetParameter(kk);
993       AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
994     }
995     
996     if(strstr(optionSave,"save") && fileProf){
997       TString profName=Form("Profile%06i",i);
998       TString timeTotName=Form("TimeTot%06i",i);
999       TString totName=Form("Tot%06i",i);
1000       TString deltaName=Form("Delta%06i",i);
1001       fileProf->cd();
1002       hSlewingProf->Write(profName);
1003       htimetot->Write(timeTotName);
1004       hToT->Write(totName);
1005       hdeltaTime->Write(deltaName);
1006     }
1007
1008     AliTOFChannelTask * calChannel = new AliTOFChannelTask();
1009     calChannel->SetSlewPar(par);
1010     arrayCal->AddAt(calChannel,i);
1011     delete hToT;
1012     hToT=0x0;
1013     delete hSlewingProf;
1014     hSlewingProf=0x0;
1015     delete htimetot;
1016     htimetot=0x0;
1017     delete hdeltaTime;
1018     hdeltaTime=0x0;
1019   }
1020
1021   if(strstr(optionSave,"save") && fileProf){
1022     fileProf->Close();
1023     delete fileProf;
1024     fileProf=0x0;
1025   }
1026   file->Close();
1027   delete file;
1028   file=0x0;
1029   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
1030   arrayCal->Write("array",TObject::kSingleKey);
1031   filecalib->Close();
1032   delete filecalib;
1033   filecalib=0x0;
1034   arrayCal->Clear();
1035   delete arrayCal;
1036   arrayCal=0x0;
1037
1038   return 0;
1039 }
1040 //----------------------------------------------------------------------------
1041 Int_t AliTOFCalibTask::CalibrateFromProfile(Int_t i, Option_t *optionSave, Option_t *optionFit){
1042
1043   // computing calibration parameters using the old profiling algo
1044   // Returning codes:
1045   // 0 -> everything was ok
1046   // 1 -> no tree for calibration found
1047   // 2 -> not enough statistics to perform calibration
1048   // 3 -> problems with arrays
1049
1050   TH1::AddDirectory(0);
1051
1052   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
1053   AliInfo(Form("*** Calibrating Histograms From Profile %s", GetName())) ; 
1054   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
1055   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
1056   const char *dirname = "./";
1057   TString filename= "TOFCalib.root";
1058
1059   if((gSystem->FindFile(dirname,filename))==NULL){
1060     AliInfo("No file with tree for calibration found! exiting....");
1061     return 1;
1062   }
1063
1064   TFile *file = new TFile("TOFCalib.root","READ");
1065   TTree *tree = (TTree*)file->Get("T");
1066   Float_t p[MAXCHENTRIESSMALL];
1067   Int_t nentries;
1068   tree->SetBranchAddress("nentries",&nentries);
1069   tree->SetBranchAddress("TOFentries",p);
1070
1071   tree->GetEntry(i);
1072   Int_t ntracks=nentries/3;
1073
1074   if (ntracks < MEANENTRIES) {  
1075     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",ntracks));
1076     return 2;
1077   }
1078
1079   TH1F * hProf = new TH1F();
1080   hProf = Profile(i);
1081   hProf->Fit("pol5",optionFit,"",0,4);
1082   TF1 * calibfunc = (TF1*)hProf->GetFunction("pol5");
1083   Float_t par[6];    
1084   for(Int_t kk=0;kk<6;kk++){
1085     par[kk]=calibfunc->GetParameter(kk);
1086     AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
1087   }
1088
1089   if(strstr(optionSave,"save")){
1090     TFile * fileProf = new TFile("TOFCalibSave.root","recreate");
1091     fileProf->cd(); 
1092     TString profName=Form("Profile%06i",i);
1093     hProf->Write(profName);
1094     fileProf->Close();
1095     delete fileProf;
1096     fileProf=0x0;
1097   }
1098
1099   delete hProf;
1100   hProf=0x0;
1101   file->Close();
1102   delete file;
1103   file=0x0;
1104   AliTOFChannelTask * calChannel = new AliTOFChannelTask();
1105   calChannel->SetSlewPar(par);
1106   arrayCal->AddAt(calChannel,i);
1107   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
1108   arrayCal->Write("array",TObject::kSingleKey);
1109   filecalib->Close();
1110   delete filecalib;
1111   filecalib = 0x0;
1112   delete calChannel;
1113   delete arrayCal;
1114   return 0;
1115 }
1116 //----------------------------------------------------------------------------
1117 Int_t AliTOFCalibTask::Calibrate(Option_t *optionSave, Option_t *optionFit){
1118
1119   // calibrating the whole TOF
1120   // computing calibration parameters
1121   // Returning codes:
1122   // 0 -> everything was ok
1123   // 1 -> no tree for calibration found
1124   // 2 -> not enough statistics to perform calibration
1125   // 3 -> problems with arrays
1126
1127   TH1::AddDirectory(0);
1128
1129   AliInfo(Form("*** Calibrating Histograms %s, all channels", GetName())) ; 
1130   AliInfo(Form("Option for Saving histos = %s",optionSave )) ; 
1131   AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ; 
1132   const char *dirname = "./";
1133   TString filename= "TOFCalib.root";
1134
1135   if((gSystem->FindFile(dirname,filename))==NULL){
1136     AliInfo("No file with tree for calibration found! exiting....");
1137     return 1;
1138   }
1139
1140   TFile * fileProf=0x0;
1141   if(strstr(optionSave,"save")){
1142     fileProf = new TFile("TOFCalibSave.root","recreate");
1143   }
1144   TFile *file = new TFile("TOFCalib.root","READ");
1145   TTree *tree = (TTree*)file->Get("T");
1146   Float_t p[MAXCHENTRIESSMALL];
1147   Int_t nentries;
1148   tree->SetBranchAddress("nentries",&nentries);
1149   tree->SetBranchAddress("TOFentries",p);
1150
1151   Float_t nentriesmean =0;
1152   for (Int_t i=0; i<TOFCHANNELS; i++){
1153     tree->GetEntry(i);
1154     Int_t ntracks=nentries/3;
1155     nentriesmean+=ntracks;
1156   }
1157
1158   nentriesmean/=TOFCHANNELS;
1159   if (nentriesmean < MEANENTRIES) {
1160     AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
1161     return 2;
1162   }
1163
1164   //filling ToT and Time arrays
1165
1166   Int_t nbinToT = 100;  // ToT bin width in Profile = 50.0 ps 
1167   Float_t minToT = 0;   // ns
1168   Float_t maxToT = 4.88;// ns
1169   TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
1170   arrayCal->SetOwner();
1171   for (Int_t ii=0; ii<TOFCHANNELS; ii++) {
1172     TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
1173     TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
1174     TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
1175     if (ii%1000 == 0) AliDebug(1,Form("Calibrating channel %i ",ii));
1176     Int_t i = 3;
1177     Int_t nusefulbins=0;
1178     Double_t binsProfile[101]; // sized larger than necessary, the correct 
1179                               // dim being set in the booking of the profile
1180     tree->GetEntry(i);
1181     Int_t ntracks=nentries/3;
1182     if (ntracks < MEANENTRIES) {
1183       AliInfo(Form(" Too small mean number of entires in channel %i (number of tracks = %f), not calibrating channel and continuing.....",i,ntracks));
1184       continue;
1185     }
1186     Float_t meantime=0;
1187     for (Int_t j=0;j<ntracks;j++){
1188       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
1189       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
1190       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
1191       Float_t tot = p[idxexToT];
1192       hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
1193       meantime+=p[idxexTime]-p[idxexExTime];
1194       hToT->Fill(tot);
1195     }
1196     nusefulbins = FindBins(hToT,&binsProfile[0]);
1197     meantime/=ntracks;
1198     for (Int_t j=0;j<nusefulbins;j++) {
1199       AliDebug(2,Form(" channel %i, usefulbin = %i, low edge = %f",i,j,binsProfile[j])); 
1200     }
1201     TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G");  // CHECK THE BUILD OPTION, PLEASE!!!!!!
1202     for (Int_t j=0;j<ntracks;j++){
1203       Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
1204       Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
1205       Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
1206       Float_t tot = p[idxexToT];
1207       Float_t time = p[idxexTime]-p[idxexExTime];
1208       AliDebug (2,Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
1209       hSlewingProf->Fill(tot,time);  // if meantime is not used, the fill may be moved in the loop above
1210       htimetot->Fill(tot,time-meantime);  // if meantime is not used, the fill may be moved in the loop above
1211     }
1212
1213     //    hSlewingProf->Fit("pol5",optionFit,"",1,4);
1214     //TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
1215     Float_t par[6];    
1216     for(Int_t kk=0;kk<6;kk++){
1217       //      par[kk]=calibfunc->GetParameter(kk);
1218       par[kk]=kk;
1219       AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
1220     }
1221
1222     if(strstr(optionSave,"save") && fileProf){
1223       TString profName=Form("Profile%06i",ii);
1224       TString timeTotName=Form("TimeTot%06i",ii);
1225       TString totName=Form("Tot%06i",ii);
1226       TString deltaName=Form("Delta%06i",ii);
1227       fileProf->cd();
1228       hSlewingProf->Write(profName);
1229       htimetot->Write(timeTotName);
1230       hToT->Write(totName);
1231       hdeltaTime->Write(deltaName);
1232     }
1233     AliTOFChannelTask * calChannel = new AliTOFChannelTask();
1234     calChannel->SetSlewPar(par);
1235     arrayCal->AddAt(calChannel,ii);
1236
1237     delete hToT;
1238     hToT=0x0;
1239     delete hSlewingProf;
1240     hSlewingProf=0x0;
1241     delete htimetot;
1242     htimetot=0x0;
1243     delete hdeltaTime;
1244     hdeltaTime=0x0;
1245   }
1246
1247   if(strstr(optionSave,"save")){
1248     fileProf->Close();
1249     delete fileProf;
1250     fileProf=0x0;
1251   }
1252
1253   file->Close();
1254   delete file;
1255   file=0x0;
1256   TFile *filecalib = new TFile("outCalArray.root","RECREATE");
1257   arrayCal->Write("array", TObject::kSingleKey);
1258   filecalib->Close();
1259   delete filecalib;
1260   filecalib=0x0;
1261   arrayCal->Clear();
1262   delete arrayCal;
1263   arrayCal=0x0;
1264
1265   return 0;
1266 }
1267
1268 //_____________________________________________________________________________
1269
1270 Bool_t AliTOFCalibTask::Select(AliESDtrack *t){
1271
1272   // selecting good quality tracks
1273   //track selection: reconstrution to TOF:
1274   if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
1275     return 0;
1276   }
1277   fnESDkTOFout++;
1278   //IsStartedTimeIntegral
1279   if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
1280     return 0;
1281   }
1282   fnESDkTIME++;
1283   if (t->GetStatus() & AliESDtrack::kTRDbackup) { 
1284     Float_t xout= t->GetOuterParam()->GetX();
1285     if (xout<364.25 &&  xout > 300.) return 0;
1286   }
1287   fnESDTRDcut++;
1288   Double_t time=t->GetTOFsignal();      
1289   time*=1.E-3; // tof given in nanoseconds
1290   if(time >= fMinTime){
1291     return 0;
1292   }
1293   fnESDTIMEcut++;
1294   
1295   Double_t mom=t->GetP();
1296   if (!(mom<=UPPERMOMBOUND && mom>=LOWERMOMBOUND)){
1297   //  return 0;  // skipping momentum cut
1298   } 
1299    
1300   UInt_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
1301   if(assignedTOFcluster==0){ // not matched
1302     return 0;
1303   }
1304   fnESDassTOFcl++;
1305   return 1;
1306 }
1307 //_____________________________________________________________________________
1308
1309 Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
1310
1311   //track Combinatorial PID for calibration
1312
1313   //cout << " In comb PID " << endl;
1314   Float_t t0offset=0;
1315   const Int_t kntracksinset=6;
1316   if (kntracksinset != 6) {
1317     AliDebug(1,"Number of tracks in set smaller than expected one, identifying every particle as if it was a pion!");
1318     return 0;
1319   }
1320   Float_t exptof[kntracksinset][3];
1321   Int_t   assparticle[kntracksinset];
1322   Float_t timeofflight[kntracksinset];
1323   Float_t texp[kntracksinset];
1324   Float_t timezero[kntracksinset];
1325   Float_t weightedtimezero[kntracksinset];
1326   Float_t besttimezero[kntracksinset];
1327   Float_t bestchisquare[kntracksinset];
1328   Float_t bestweightedtimezero[kntracksinset];
1329
1330   for (Int_t i=0;i<kntracksinset;i++){
1331     assparticle[i]=3;
1332     timeofflight[i]=0;
1333     texp[i]=0;
1334     timezero[i]=0;
1335     weightedtimezero[i]=0;
1336     besttimezero[i]=0;
1337     bestchisquare[i]=0;
1338     bestweightedtimezero[i]=0;
1339   }
1340
1341   Int_t nset= (Int_t)(size/kntracksinset/NIDX);
1342   for (Int_t i=0; i< nset; i++) {   
1343     for (Int_t j=0; j<kntracksinset; j++) {
1344       Int_t idxtime = ((kntracksinset*i+j)*NIDX)+DELTAIDXTIME; 
1345       Int_t idxextimePi = ((kntracksinset*i+j)*NIDX)+DELTAIDXEXTIMEPI; 
1346       Int_t idxextimeKa = ((kntracksinset*i+j)*NIDX)+DELTAIDXEXTIMEKA; 
1347       Int_t idxextimePr = ((kntracksinset*i+j)*NIDX)+DELTAIDXEXTIMEPR; 
1348       Double_t time=smallarray[idxtime]; // TOF time in ns
1349       timeofflight[j]=time+t0offset;
1350       exptof[j][0]=smallarray[idxextimePi];
1351       exptof[j][1]=smallarray[idxextimeKa];
1352       exptof[j][2]=smallarray[idxextimePr];
1353       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]));
1354     }
1355     Float_t t0best=999.;
1356     Float_t et0best=999.;
1357     Float_t chisquarebest=999.;
1358     for (Int_t i1=0; i1<3;i1++) {
1359       texp[0]=exptof[0][i1];
1360       for (Int_t i2=0; i2<3;i2++) { 
1361         texp[1]=exptof[1][i2];
1362         for (Int_t i3=0; i3<3;i3++) {
1363           texp[2]=exptof[2][i3];
1364           for (Int_t i4=0; i4<3;i4++) {
1365             texp[3]=exptof[3][i4];
1366             for (Int_t i5=0; i5<3;i5++) {
1367               texp[4]=exptof[4][i5];
1368               for (Int_t i6=0; i6<3;i6++) {
1369                 texp[5]=exptof[5][i6];
1370                                         
1371                 Float_t sumAllweights=0.;
1372                 Float_t meantzero=0.;
1373                 Float_t emeantzero=0.;
1374                 
1375                 for (Int_t itz=0; itz<kntracksinset;itz++) {
1376                   timezero[itz]=texp[itz]-timeofflight[itz];                
1377                   weightedtimezero[itz]=timezero[itz]/TRACKERROR;
1378                   sumAllweights+=1./TRACKERROR;
1379                   meantzero+=weightedtimezero[itz];
1380                   
1381                 } // end loop for (Int_t itz=0; itz<15;itz++)
1382                 
1383                 meantzero=meantzero/sumAllweights; // it is given in [ns]
1384                 emeantzero=sqrt(1./sumAllweights); // it is given in [ns]
1385                 
1386                 // calculate chisquare
1387                 
1388                 Float_t chisquare=0.;           
1389                 for (Int_t icsq=0; icsq<kntracksinset;icsq++) {
1390                   chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/TRACKERROR;
1391                 } // end loop for (Int_t icsq=0; icsq<15;icsq++) 
1392                 
1393                 Int_t npion=0;
1394                 if(i1==0)npion++;
1395                 if(i2==0)npion++;
1396                 if(i3==0)npion++;
1397                 if(i4==0)npion++;
1398                 if(i5==0)npion++;
1399                 if(i6==0)npion++;
1400                 
1401                 if(chisquare<=chisquarebest  && ((Float_t) npion/ ((Float_t) kntracksinset)>0.3)){
1402                   for(Int_t iqsq = 0; iqsq<kntracksinset; iqsq++) {
1403                     besttimezero[iqsq]=timezero[iqsq]; 
1404                     bestweightedtimezero[iqsq]=weightedtimezero[iqsq]; 
1405                     bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/TRACKERROR; 
1406                   }
1407                   
1408                   assparticle[0]=i1;
1409                   assparticle[1]=i2;
1410                   assparticle[2]=i3;
1411                   assparticle[3]=i4;
1412                   assparticle[4]=i5;
1413                   assparticle[5]=i6;
1414                   
1415                   chisquarebest=chisquare;
1416                   t0best=meantzero;
1417                   et0best=emeantzero;
1418                 } // close if(dummychisquare<=chisquare)
1419               } // end loop on i6
1420             } // end loop on i5
1421           } // end loop on i4
1422         } // end loop on i3
1423       } // end loop on i2
1424     } // end loop on i1
1425                 
1426     Float_t confLevel=999;
1427     if(chisquarebest<999.){
1428       Double_t dblechisquare=(Double_t)chisquarebest;
1429       confLevel=(Float_t)TMath::Prob(dblechisquare,kntracksinset-1); 
1430     }
1431     // assume they are all pions for fake sets
1432     if(confLevel<0.01 || confLevel==999. ){
1433       for (Int_t itrk=0; itrk<kntracksinset; itrk++)assparticle[itrk]=0;
1434     }
1435
1436     AliDebug(2,Form(" Best Assignment, selection %i %i %i %i %i %i ",assparticle[0],assparticle[1],assparticle[2],assparticle[3],assparticle[4],assparticle[5]));
1437
1438     for (Int_t kk=0;kk<kntracksinset;kk++){
1439       Int_t idxextimePi = ((kntracksinset*i+kk)*NIDX)+DELTAIDXEXTIMEPI; 
1440       Int_t idxextimeKa = ((kntracksinset*i+kk)*NIDX)+DELTAIDXEXTIMEKA; 
1441       Int_t idxextimePr = ((kntracksinset*i+kk)*NIDX)+DELTAIDXEXTIMEPR; 
1442       // storing in third slot of smallarray the assigned expected time
1443       fhPID->Fill(assparticle[kk]);
1444       //      assparticle[kk]=0;  //assuming all particles are pions
1445       if (assparticle[kk]==0){       //assigned to be a Pi
1446         smallarray[idxextimeKa]=-1;
1447         smallarray[idxextimePr]=-1;
1448       }
1449       else if (assparticle[kk]==1){  //assigned to be a Ka
1450         smallarray[idxextimePi]=smallarray[idxextimeKa];
1451         smallarray[idxextimeKa]=-1;
1452         smallarray[idxextimePr]=-1;
1453       }
1454       else if (assparticle[kk]==2){  //assigned to be a Pr
1455         smallarray[idxextimePi]=smallarray[idxextimePr];
1456         smallarray[idxextimeKa]=-1;
1457         smallarray[idxextimePr]=-1;
1458       }
1459     }
1460   }
1461   return 1;
1462 }
1463 //-----------------------------------------------------------------------
1464 TH1F* AliTOFCalibTask::Profile(Int_t ich)
1465 {
1466   // profiling algo
1467
1468   TFile *file = new TFile("TOFCalib.root","READ");
1469   TTree *tree = (TTree*)file->Get("T");
1470   Float_t p[MAXCHENTRIESSMALL];
1471   Int_t nentries;
1472   tree->SetBranchAddress("nentries",&nentries);
1473   tree->SetBranchAddress("TOFentries",p);
1474
1475   //Prepare histograms for Slewing Correction
1476   const Int_t knbinToT = 100;
1477   Int_t nbinTime = 200;
1478   Float_t minTime = -5.5; //ns
1479   Float_t maxTime = 5.5; //ns
1480   Float_t minToT = 0; //ns
1481   Float_t maxToT = 5.; //ns
1482   Float_t deltaToT = (maxToT-minToT)/knbinToT;
1483   Double_t mTime[knbinToT+1],mToT[knbinToT+1],meanTime[knbinToT+1], meanTime2[knbinToT+1],vToT[knbinToT+1], vToT2[knbinToT+1],meanToT[knbinToT+1],meanToT2[knbinToT+1],vTime[knbinToT+1],vTime2[knbinToT+1],xlow[knbinToT+1],sigmaTime[knbinToT+1];
1484   Int_t n[knbinToT+1], nentrx[knbinToT+1];
1485   Double_t sigmaToT[knbinToT+1];
1486   for (Int_t i = 0; i < knbinToT+1 ; i++){
1487     mTime[i]=0;
1488     mToT[i]=0;
1489     n[i]=0;
1490     meanTime[i]=0;
1491     meanTime2[i]=0;
1492     vToT[i]=0;
1493     vToT2[i]=0;
1494     meanToT[i]=0;
1495     meanToT2[i]=0;
1496     vTime[i]=0;
1497     vTime2[i]=0;
1498     xlow[i]=0;
1499     sigmaTime[i]=0;
1500     sigmaToT[i]=0;
1501     n[i]=0;
1502     nentrx[i]=0;
1503   }
1504   TH2F* hSlewing = new TH2F("hSlewing", "hSlewing", knbinToT, minToT, maxToT, nbinTime, minTime, maxTime);
1505   Int_t ntracks = 0;
1506   TH1F *histo = new TH1F("histo", "1D Time vs ToT", knbinToT, minToT, maxToT);
1507   tree->GetEntry(ich);
1508   ntracks=nentries/3;
1509   for (Int_t j=0;j<ntracks;j++){
1510     Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT; 
1511     Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME; 
1512     Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID; 
1513     Float_t tot = p[idxexToT];
1514     Float_t time = p[idxexTime]-p[idxexExTime];
1515     Int_t nx = (Int_t)((tot-minToT)/deltaToT)+1;
1516     if ((tot != 0) && ( time!= 0)){
1517       vTime[nx]+=time;
1518       vTime2[nx]+=time*time;
1519       vToT[nx]+=tot;
1520       vToT2[nx]+=tot*tot;
1521       nentrx[nx]++;
1522       hSlewing->Fill(tot,time);
1523     }
1524   }
1525   Int_t nbinsToT=hSlewing->GetNbinsX();
1526   if (nbinsToT != knbinToT) {
1527     AliError("Profile :: incompatible numbers of bins");
1528     return 0x0;
1529   }
1530   
1531   Int_t usefulBins=0;
1532   for (Int_t i=1;i<=nbinsToT;i++){
1533     if (nentrx[i]!=0){
1534       n[usefulBins]+=nentrx[i];
1535       if (n[usefulBins]==0 && i == nbinsToT) {
1536         break;
1537       }
1538       meanTime[usefulBins]+=vTime[i];
1539       meanTime2[usefulBins]+=vTime2[i];
1540       meanToT[usefulBins]+=vToT[i];
1541       meanToT2[usefulBins]+=vToT2[i];
1542       if (n[usefulBins]<10 && i!=nbinsToT) continue; 
1543       mTime[usefulBins]=meanTime[usefulBins]/n[usefulBins];
1544       mToT[usefulBins]=meanToT[usefulBins]/n[usefulBins];
1545       sigmaTime[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
1546                             *(meanTime2[usefulBins]-meanTime[usefulBins]
1547                             *meanTime[usefulBins]/n[usefulBins]));
1548       if ((1./n[usefulBins]/n[usefulBins]
1549            *(meanToT2[usefulBins]-meanToT[usefulBins]
1550              *meanToT[usefulBins]/n[usefulBins]))< 0) {
1551         AliError(" too small radical" );
1552         sigmaToT[usefulBins]=0;
1553       }
1554       else{       
1555         sigmaToT[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
1556                              *(meanToT2[usefulBins]-meanToT[usefulBins]
1557                              *meanToT[usefulBins]/n[usefulBins]));
1558       }
1559       usefulBins++;
1560     }
1561   }
1562   for (Int_t i=0;i<usefulBins;i++){
1563     Int_t binN = (Int_t)((mToT[i]-minToT)/deltaToT)+1;
1564     histo->Fill(mToT[i],mTime[i]);
1565     histo->SetBinError(binN,sigmaTime[i]);
1566   } 
1567   delete file;
1568   file=0x0;
1569   delete hSlewing;
1570   hSlewing=0x0;
1571
1572   return histo;
1573 }
1574 //----------------------------------------------------------------------------
1575 Int_t AliTOFCalibTask::FindBins(TH1F* h, Double_t *binsProfile) const{
1576
1577   // to determine the bins for ToT histo
1578
1579   Int_t cont = 0;
1580   Int_t startBin = 1;
1581   Int_t nbin = h->GetNbinsX();
1582   Int_t nentries = (Int_t)h->GetEntries();
1583   Float_t max = h->GetBinLowEdge(nbin);
1584   Int_t nusefulbins=0;
1585   Int_t maxcont=0;
1586   // setting maxvalue of entries per bin
1587   if (nentries <= 60) maxcont = 2;
1588   else  if (nentries <= 100) maxcont = 5;
1589   else  if (nentries <= 500) maxcont = 10;
1590   else  maxcont = 20;
1591   for (Int_t j=1;j<=nbin;j++) {
1592     cont += (Int_t)h->GetBinContent(j);
1593     if (j<nbin){
1594       if (cont>=maxcont){
1595         nusefulbins++;
1596         binsProfile[nusefulbins-1]=h->GetBinLowEdge(startBin);
1597         cont=0;
1598         startBin=j+1;
1599         continue;
1600       }
1601     }
1602     else{
1603       if (cont>=maxcont){
1604         nusefulbins++;
1605         binsProfile[nusefulbins-1]=h->GetBinLowEdge(startBin);
1606         binsProfile[nusefulbins]=max;
1607       }
1608       else {
1609         binsProfile[nusefulbins]=h->GetBinLowEdge(startBin);
1610       }
1611     }
1612   }
1613   return nusefulbins;
1614 }