e8bcc6bb197e222216f46d6e72f4cf4995e30ba6
[u/mrichter/AliRoot.git] / TPC / Calib / AliTPCcalibTimeGain.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 /*
17
18 This class provides the calibration of the time dependence of the TPC gain due to pressure and temperature changes etc.
19
20
21 //0.  Libraries to load
22 gSystem->Load("libANALYSIS");
23 gSystem->Load("libSTAT");
24 gSystem->Load("libTPCcalib");
25
26
27 //1. Do calibration ...
28 //
29 // compare reference
30
31 //
32 //2. Visualize results
33 //
34 TFile fcalib("CalibObjects.root");
35 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
36 AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
37 TGraphErrors * gr = gain->GetGraphGainVsTime(0,1000)
38
39 // gain->GetHistGainTime()->GetAxis(1)->SetRangeUser(1213.8e6,1214.3e6)
40 TH2D * GainTime = gain->GetHistGainTime()->Projection(0,1)
41 GainTime->GetXaxis()->SetTimeDisplay(kTRUE)
42 GainTime->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}")
43 GainTime->Draw("colz")
44
45 //
46 // MakeSlineFit example
47 //
48 AliSplineFit *fit = AliTPCcalibTimeGain::MakeSplineFit(gr)
49
50 TGraph * grfit = fit.MakeGraph(gr->GetX()[0],gr->GetX()[gr->GetN()-1],50000,0);
51
52 gr->SetMarkerStyle(25);
53 gr->Draw("lp");
54 grfit->SetLineColor(2);
55 grfit->Draw("lu");
56
57 //
58 // QA - dE/dx resoultion as a function of time
59 //TCa
60
61 TGraph * grSigma = AliTPCcalibBase::FitSlicesSigma(gain->GetHistGainTime(),0,1,1800,5)
62
63 TCanvas *c1 = new TCanvas("c1","transparent pad",200,10,700,500);
64    TPad *pad1 = new TPad("pad1","",0,0,1,1);
65    TPad *pad2 = new TPad("pad2","",0,0,1,1);
66    pad2->SetFillStyle(4000); //will be transparent
67    pad1->Draw();
68    pad1->cd();
69
70 GainTime->Draw("colz")
71 gr->Draw("lp")
72
73
74
75   c1->cd();
76  Double_t ymin = -0.04;
77  Double_t ymax = 0.12;
78 Double_t dy = (ymax-ymin)/0.8;
79 Double_t xmin = GainTime->GetXaxis()->GetXmin()
80 Double_t xmax = GainTime->GetXaxis()->GetXmax()
81 Double_t dx = (xmax-xmin)/0.8; //10 per cent margins left and right
82 pad2->Range(xmin-0.1*dx,ymin-0.1*dy,xmax+0.1*dx,ymax+0.1*dy);
83    pad2->Draw();
84    pad2->cd();
85 grSigma->SetLineColor(2);
86 grSigma->SetLineWidth(2);
87 grSigma->Draw("lp")
88 TGaxis *axis = new TGaxis(xmax,ymin,xmax,ymax,ymin,ymax,50510,"+L");
89    axis->SetLabelColor(kRed);
90    axis->SetTitle("dE/dx resolution #sigma_{dE/dx}");
91    axis->Draw();
92
93 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
94
95  ----> make Attachment study
96
97 TFile fcalib("CalibObjects40366b.root");
98 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
99 AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
100 TGraphErrors * grAttach = gain->GetGraphAttachment(2000,4)
101
102 TCanvas *c1 = new TCanvas("c1","transparent pad",200,10,700,500);
103    TPad *pad1 = new TPad("pad1","",0,0,1,1);
104    TPad *pad2 = new TPad("pad2","",0,0,1,1);
105    pad2->SetFillStyle(4000); //will be transparent
106    pad1->Draw();
107    pad1->cd();
108
109 gain->GetHistGainTime()->GetAxis(1)->SetRangeUser(1213.8e6,1214.3e6)
110 TH2D * GainTime = gain->GetHistGainTime()->Projection(0,1)
111 GainTime->GetXaxis()->SetTimeDisplay(kTRUE)
112 GainTime->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}")
113 GainTime->Draw("colz")
114 //gr->Draw("lp")
115
116   c1->cd();
117  Double_t ymin = -0.001;
118  Double_t ymax = 0.001;
119 Double_t dy = (ymax-ymin)/0.8;
120 Double_t xmin = GainTime->GetXaxis()->GetXmin()
121 Double_t xmax = GainTime->GetXaxis()->GetXmax()
122 Double_t dx = (xmax-xmin)/0.8; //10 per cent margins left and right
123 pad2->Range(xmin-0.1*dx,ymin-0.1*dy,xmax+0.1*dx,ymax+0.1*dy);
124    pad2->Draw();
125    pad2->cd();
126 grAttach->SetLineColor(2);
127 grAttach->SetLineWidth(2);
128 grAttach->Draw("lp")
129 TGaxis *axis = new TGaxis(xmax,ymin,xmax,ymax,ymin,ymax,50510,"+L");
130    axis->SetLabelColor(kRed);
131    axis->SetTitle("attachment coefficient b");
132    axis->Draw();
133
134
135 */
136
137
138 #include "Riostream.h"
139 #include "TChain.h"
140 #include "TTree.h"
141 #include "TH1F.h"
142 #include "TH2F.h"
143 #include "TH3F.h"
144 #include "THnSparse.h"
145 #include "TGraphErrors.h"
146 #include "TList.h"
147 #include "TMath.h"
148 #include "TCanvas.h"
149 #include "TFile.h"
150 #include "TF1.h"
151 #include "TVectorD.h"
152 #include "TProfile.h"
153 #include "TGraphErrors.h"
154 #include "TCanvas.h"
155
156 #include "AliTPCclusterMI.h"
157 #include "AliTPCseed.h"
158 #include "AliESDVertex.h"
159 #include "AliESDEvent.h"
160 #include "AliESDfriend.h"
161 #include "AliESDInputHandler.h"
162 #include "AliAnalysisManager.h"
163
164 #include "AliTracker.h"
165 #include "AliMagF.h"
166 #include "AliTPCCalROC.h"
167 #include "AliESDv0.h"
168
169 #include "AliLog.h"
170
171 #include "AliTPCcalibTimeGain.h"
172
173 #include "TTreeStream.h"
174 #include "AliTPCTracklet.h"
175 #include "TTimeStamp.h"
176 #include "AliTPCcalibDB.h"
177 #include "AliTPCcalibLaser.h"
178 #include "AliDCSSensorArray.h"
179 #include "AliDCSSensor.h"
180
181 ClassImp(AliTPCcalibTimeGain)
182
183 Double_t AliTPCcalibTimeGain::fgMergeEntriesCut=10000000.;
184
185 AliTPCcalibTimeGain::AliTPCcalibTimeGain() 
186   :AliTPCcalibBase(), 
187    fHistGainTime(0),
188    fGainVsTime(0),
189    fHistDeDxTotal(0),
190    fIntegrationTimeDeDx(0),
191    fMIP(0),
192    fCutCrossRows(0),
193    fCutEtaWindow(0),
194    fCutRequireITSrefit(0),
195    fCutMaxDcaXY(0),
196    fCutMaxDcaZ(0),
197    fUseMax(0),
198    fLowerTrunc(0),
199    fUpperTrunc(0),
200    fUseShapeNorm(0),
201    fUsePosNorm(0),
202    fUsePadNorm(0),
203    fUseCookAnalytical(0),
204    fIsCosmic(0),
205    fLowMemoryConsumption(0)
206 {  
207   //
208   // Default constructor
209   //
210   AliDebug(5,"Default Constructor");  
211 }
212
213
214 AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain)
215   :AliTPCcalibBase(),
216    fHistGainTime(0),
217    fGainVsTime(0),
218    fHistDeDxTotal(0),
219    fIntegrationTimeDeDx(0),
220    fMIP(0),
221    fCutCrossRows(0),
222    fCutEtaWindow(0),
223    fCutRequireITSrefit(0),
224    fCutMaxDcaXY(0),
225    fCutMaxDcaZ(0),
226    fUseMax(0),
227    fLowerTrunc(0),
228    fUpperTrunc(0),
229    fUseShapeNorm(0),
230    fUsePosNorm(0),
231    fUsePadNorm(0),
232    fUseCookAnalytical(0),
233    fIsCosmic(0),
234    fLowMemoryConsumption(0)
235 {
236   //
237   // No default constructor
238   //
239   SetName(name);
240   SetTitle(title);
241   
242   AliDebug(5,"Non Default Constructor");
243   
244   fIntegrationTimeDeDx = deltaIntegrationTimeGain;
245   
246   Double_t deltaTime = EndTime - StartTime;
247   
248   
249   // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 - proton points at higher dE/dx), meanDriftlength, momenta (only filled if enough space is available), run number, eta
250   Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
251   Int_t binsGainTime[7]    = {150,  timeBins,    4,  25, 200, 10000000, 20};
252   Double_t xminGainTime[7] = {0.5, StartTime,  0.5,   0, 0.1,    -0.5,  -1};
253   Double_t xmaxGainTime[7] = {  8,   EndTime,  4.5, 250,  50, 9999999.5, 1};
254   fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number, eta;dEdx",7,binsGainTime,xminGainTime,xmaxGainTime);
255   BinLogX(fHistGainTime, 4);
256   //
257   fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,8);
258   BinLogX(fHistDeDxTotal);
259   
260   
261   // default track selection cuts
262   fCutCrossRows = 80;
263   fCutEtaWindow = 0.8;
264   fCutRequireITSrefit = kFALSE;
265   fCutMaxDcaXY = 3.5;
266   fCutMaxDcaZ  = 25.;
267
268   // default values for dE/dx
269   fMIP = 50.;
270   fUseMax = kTRUE;
271   fLowerTrunc = 0.02;
272   fUpperTrunc = 0.6;
273   fUseShapeNorm = kTRUE;
274   fUsePosNorm = kFALSE;
275   fUsePadNorm = kFALSE;
276   fUseCookAnalytical = kFALSE;
277   //
278   fIsCosmic = kTRUE;
279   fLowMemoryConsumption = kTRUE;
280   //
281   
282 }
283
284
285
286 AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
287   //
288   // Destructor
289   //
290   delete fHistGainTime;
291   delete fGainVsTime;
292   delete fHistDeDxTotal;
293 }
294
295
296 void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
297   //
298   // main track loop
299   //
300   if (!event) {
301     Printf("ERROR: ESD not available");
302     return;
303   }
304   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
305   if (!ESDfriend) {
306    return;
307   }
308   if (ESDfriend->TestSkipBit()) return;
309
310   if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
311     ProcessCosmicEvent(event);
312   } else {
313     ProcessBeamEvent(event);
314   }
315   
316
317   
318   
319 }
320
321
322 void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
323   //
324   // Process in case of cosmic event
325   //
326   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
327   if (!esdFriend) {
328    Printf("ERROR: ESDfriend not available");
329    return;
330   }
331   //
332   UInt_t time = event->GetTimeStamp();
333   Int_t nFriendTracks = esdFriend->GetNumberOfTracks();
334   Int_t runNumber = event->GetRunNumber();
335   //
336   // track loop
337   //
338   for (Int_t i=0;i<nFriendTracks;++i) {
339
340     AliESDtrack *track = event->GetTrack(i);
341     if (!track) continue;
342     AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
343     if (!friendTrack) continue;        
344     const AliExternalTrackParam * trackIn = track->GetInnerParam();
345     const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
346     if (!trackIn) continue;
347     if (!trackOut) continue;
348
349     // calculate necessary track parameters
350     Double_t meanP = trackIn->GetP();
351     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
352     Int_t nclsDeDx = track->GetTPCNcls();
353
354     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
355     if (nclsDeDx < 60) continue;     
356     if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
357     if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
358     
359     // Get seeds
360     TObject *calibObject;
361     AliTPCseed *seed = 0;
362     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
363       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
364     }    
365
366     if (seed) { 
367       Double_t tpcSignal = GetTPCdEdx(seed);
368       if (nclsDeDx > 100) fHistDeDxTotal->Fill(meanP, tpcSignal);
369       //
370       if (fLowMemoryConsumption) {
371         if (meanP < 20) continue;
372         meanP = 30; // set all momenta to one in order to save memory
373       }
374       //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
375       Double_t vec[6] = {tpcSignal,time,1,meanDrift,meanP,runNumber};
376       fHistGainTime->Fill(vec);
377     }
378     
379   }
380
381 }
382
383
384
385 void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
386   //
387   // Process in case of beam event
388   //
389   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
390   if (!esdFriend) {
391    Printf("ERROR: ESDfriend not available");
392    return;
393   }
394   //
395   UInt_t time = event->GetTimeStamp();
396   Int_t nFriendTracks = esdFriend->GetNumberOfTracks();
397   Int_t runNumber = event->GetRunNumber();
398   //
399   // track loop
400   //
401   for (Int_t i=0;i<nFriendTracks;++i) { // begin track loop
402
403     AliESDtrack *track = event->GetTrack(i);
404     if (!track) continue;
405     AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
406     if (!friendTrack) continue;
407         
408     const AliExternalTrackParam * trackIn = track->GetInnerParam();
409     const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
410     if (!trackIn) continue;
411     if (!trackOut) continue;
412
413     // calculate necessary track parameters
414     Double_t meanP = trackIn->GetP();
415     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
416     Int_t nCrossedRows = track->GetTPCCrossedRows();
417
418     //
419     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
420     // fCutCrossRows = 80, fCutEtaWindow = 0.8, fCutRequireITSrefit, fCutMaxDcaXY = 3.5, fCutMaxDcaZ = 25
421     //
422     if (nCrossedRows < fCutCrossRows) continue;     
423     if (TMath::Abs(trackIn->Eta()) > fCutEtaWindow) continue;
424     //
425     UInt_t status = track->GetStatus();
426     if ((status&AliESDtrack::kTPCrefit)==0) continue;
427     if ((status&AliESDtrack::kITSrefit)==0 && fCutRequireITSrefit) continue; // ITS cluster
428     //
429     Float_t dca[2], cov[3];
430     track->GetImpactParameters(dca,cov);
431     if (TMath::Abs(dca[0]) > fCutMaxDcaXY || TMath::Abs(dca[0]) < 0.0000001) continue;  // cut in xy
432     if (((status&AliESDtrack::kITSrefit) == 1 && TMath::Abs(dca[1]) > 3.) || TMath::Abs(dca[1]) > fCutMaxDcaZ ) continue;
433     //
434     Double_t eta = trackIn->Eta();
435     
436     // Get seeds
437     TObject *calibObject;
438     AliTPCseed *seed = 0;
439     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
440       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
441     }    
442
443     if (seed) {
444       Int_t particleCase = 0;
445       if (meanP < 0.55 && meanP > 0.5)  particleCase = 2; // MIP pions
446       if (meanP < 0.57 && meanP > 0.56) particleCase = 3; // protons 1
447       if (meanP < 0.66 && meanP > 0.65) particleCase = 4; // protons 2
448       //
449       if (fLowMemoryConsumption && particleCase == 0) continue;
450       //
451       Double_t tpcSignal = GetTPCdEdx(seed);
452       /*
453       if (particleCalse == 0) {
454         Float_t corrFactor = AliExternalTrackParam::BetheBlochAleph(x/0.13957);
455           tpcSignal /= corrFactor; 
456       }
457       */
458       fHistDeDxTotal->Fill(meanP, tpcSignal);
459       //
460       //dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 protons), momenta, runNumner, eta
461       Double_t vec[7] = {tpcSignal,time,particleCase,meanDrift,meanP,runNumber, eta};
462       fHistGainTime->Fill(vec);
463     }
464     
465   } // end track loop
466   //
467   // V0 loop -- in beam events the cosmic part of the histogram is filled with GammaConversions
468   //
469   for(Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) {
470     AliESDv0 * v0 = event->GetV0(iv0);
471     if (!v0->GetOnFlyStatus()) continue;
472     if (v0->GetEffMass(0,0) > 0.02) continue; // select low inv. mass
473     Double_t xyz[3];
474     v0->GetXYZ(xyz[0], xyz[1], xyz[2]);
475     if (TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) < 3 || TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) > 30) continue;
476     //
477     // "loop over daughters" 
478     //
479     for(Int_t idaughter = 0; idaughter < 2; idaughter++) { // daughter loop
480       Int_t index = idaughter == 0 ? v0->GetPindex() : v0->GetNindex();
481       AliESDtrack * trackP = event->GetTrack(index);
482       AliESDfriendTrack *friendTrackP = esdFriend->GetTrack(index);
483       if (!friendTrackP) continue;
484       const AliExternalTrackParam * trackPIn = trackP->GetInnerParam();
485       const AliExternalTrackParam * trackPOut = friendTrackP->GetTPCOut();
486       if (!trackPIn) continue;
487       if (!trackPOut) continue;
488       // calculate necessary track parameters
489       Double_t meanP = trackPIn->GetP();
490       Double_t meanDrift = 250 - 0.5*TMath::Abs(trackPIn->GetZ() + trackPOut->GetZ());
491       Int_t nclsDeDx = trackP->GetTPCNcls();
492       // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
493       if (nclsDeDx < 60) continue;     
494       if (TMath::Abs(trackPIn->GetTgl()) > 1) continue;
495       //
496       TObject *calibObject;
497       AliTPCseed *seed = 0;
498       for (Int_t l=0;(calibObject=friendTrackP->GetCalibObject(l));++l) {
499       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
500       }    
501       if (seed) { 
502         if (fLowMemoryConsumption) {
503           if (meanP > 0.5 || meanP < 0.4) continue;
504           meanP = 0.45; // set all momenta to one in order to save memory
505       }
506         Double_t tpcSignal = GetTPCdEdx(seed);
507         //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
508         Double_t vec[6] = {tpcSignal,time,1,meanDrift,meanP,runNumber};
509         fHistGainTime->Fill(vec);
510       }
511     }
512     
513   }
514
515 }
516
517
518 Float_t AliTPCcalibTimeGain::GetTPCdEdx(AliTPCseed * seed) {
519   //
520   // calculate tpc dEdx
521   //
522   Double_t signal = 0;
523   //
524   if (!fUseCookAnalytical) {
525     signal = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,fUseMax,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
526   } else {
527     signal = (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax);
528   }
529   //
530   return signal;
531 }
532
533
534 void AliTPCcalibTimeGain::AnalyzeRun(Int_t minEntries) {
535   //
536   // Analyze results of calibration
537   //
538   if (fIsCosmic) {
539     fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
540     fHistGainTime->GetAxis(4)->SetRangeUser(20,100);    // only Fermi-Plateau muons
541   } else {
542     fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
543     fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
544   }
545   //
546   fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,minEntries,10);
547   //
548   return;
549 }
550
551
552 TGraphErrors * AliTPCcalibTimeGain::GetGraphGainVsTime(Int_t runNumber, Int_t minEntries) {
553   //
554   // Analyze results and get the graph 
555   //
556   if (runNumber == 0) {
557     if (!fGainVsTime) {
558       AnalyzeRun(minEntries);
559     }
560   } else {
561     // 1st check if the current run was cosmic or beam event
562     fHistGainTime->GetAxis(5)->SetRangeUser(runNumber,runNumber);
563     AnalyzeRun(minEntries);
564   }
565   if (fGainVsTime->GetN() == 0) return 0;
566   return fGainVsTime;
567 }
568
569 Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
570   //
571   // merge component
572   //
573   TIterator* iter = li->MakeIterator();
574   AliTPCcalibTimeGain* cal = 0;
575
576   while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
577     if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
578       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
579       return -1;
580     }
581
582     // add histograms here...
583     if (cal->GetHistGainTime() && fHistGainTime ) 
584     {
585       if ((fHistGainTime->GetEntries() + cal->GetHistGainTime()->GetEntries()) < fgMergeEntriesCut) 
586       {
587         //AliInfo(Form("fHistGainTime has %.0f tracks, going to add %.0f\n",fHistGainTime->GetEntries(),cal->GetHistGainTime()->GetEntries()));
588         fHistGainTime->Add(cal->GetHistGainTime());
589       }
590       else
591       {
592         AliInfo(Form("fHistGainTime full (has %.0f merged tracks, max allowed: %.0f)",fHistGainTime->GetEntries(),fgMergeEntriesCut));
593       }
594     }
595
596     if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
597
598   }
599   delete iter;
600   return 0;
601   
602 }
603
604
605 AliSplineFit * AliTPCcalibTimeGain::MakeSplineFit(TGraphErrors * graph) {
606   //
607   // make spline fit of gain
608   //
609   AliSplineFit *fit = new AliSplineFit();
610   fit->SetGraph(graph);
611   fit->SetMinPoints(graph->GetN()+1);
612   fit->InitKnots(graph,2,0,0.001);
613   fit->SplineFit(0);
614   return fit;
615   
616 }
617
618
619
620 TGraphErrors * AliTPCcalibTimeGain::GetGraphAttachment(Int_t minEntries, Int_t nmaxBin, Float_t /*fracLow*/, Float_t /*fracUp*/) {
621   //
622   // For each time bin the driftlength-dependence of the signal is fitted.
623   //
624   TH3D * hist = fHistGainTime->Projection(1, 0, 3);
625   Double_t *xvec = new Double_t[hist->GetNbinsX()];
626   Double_t *yvec = new Double_t[hist->GetNbinsX()];
627   Double_t *xerr = new Double_t[hist->GetNbinsX()];
628   Double_t *yerr = new Double_t[hist->GetNbinsX()];
629   Int_t counter  = 0;
630   TH2D * projectionHist = 0x0;
631   //
632   for(Int_t i=1; i < hist->GetNbinsX(); i++) {
633     Int_t nsum=0;
634     Int_t imin   =  i;
635     Int_t imax   =  i;    
636     for (Int_t idelta=0; idelta<nmaxBin; idelta++){
637       //
638       imin   =  TMath::Max(i-idelta,1);
639       imax   =  TMath::Min(i+idelta,hist->GetNbinsX());
640       nsum = TMath::Nint(hist->Integral(imin,imax,1,hist->GetNbinsY()-1,1,hist->GetNbinsZ()-1));
641       if (nsum==0) break;
642       if (nsum>minEntries) break;
643     }
644     if (nsum<minEntries) continue;
645     //
646     hist->GetXaxis()->SetRange(imin,imax);
647     projectionHist = (TH2D*)hist->Project3D("yzNUFNOF");
648     //
649     TObjArray arr;
650     projectionHist->FitSlicesY(0,2, projectionHist->GetNbinsX()-2,0,"QNR",&arr);
651     TH1D * histAttach = (TH1D*)arr.At(1);
652     TF1 pol1("polynom1","pol1",10,240);
653     histAttach->Fit(&pol1,"QNR");
654     xvec[counter] = 0.5*(hist->GetXaxis()->GetBinCenter(imin)+hist->GetXaxis()->GetBinCenter(imax));
655     yvec[counter] = pol1.GetParameter(1)/pol1.GetParameter(0);
656     xerr[counter] = 0;
657     yerr[counter] = pol1.GetParError(1)/pol1.GetParameter(0);
658     counter++;
659     //
660     delete projectionHist;
661   }
662   
663   TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
664   delete [] xvec;
665   delete [] yvec;
666   delete [] xerr;
667   delete [] yerr;
668   delete hist;
669   return graphErrors;
670   
671 }
672
673
674
675 void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
676   //
677   // Method for the correct logarithmic binning of histograms
678   //
679   TAxis *axis = h->GetAxis(axisDim);
680   int bins = axis->GetNbins();
681
682   Double_t from = axis->GetXmin();
683   Double_t to = axis->GetXmax();
684   Double_t *newBins = new Double_t[bins + 1];
685    
686   newBins[0] = from;
687   Double_t factor = pow(to/from, 1./bins);
688   
689   for (int i = 1; i <= bins; i++) {
690    newBins[i] = factor * newBins[i-1];
691   }
692   axis->Set(bins, newBins);
693   delete[] newBins;
694   
695 }
696
697
698 void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
699   //
700   // Method for the correct logarithmic binning of histograms
701   //
702   TAxis *axis = h->GetXaxis();
703   int bins = axis->GetNbins();
704
705   Double_t from = axis->GetXmin();
706   Double_t to = axis->GetXmax();
707   Double_t *newBins = new Double_t[bins + 1];
708    
709   newBins[0] = from;
710   Double_t factor = pow(to/from, 1./bins);
711   
712   for (int i = 1; i <= bins; i++) {
713    newBins[i] = factor * newBins[i-1];
714   }
715   axis->Set(bins, newBins);
716   delete[] newBins;
717   
718 }
719
720
721
722 void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
723   //
724   // Fit the bethe bloch params
725   //
726   //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
727   const Double_t sigma = 0.06;
728
729   TH2F *histBG = new TH2F("histBG","dEdxBg; #beta #gamma; TPC signal (a.u.)",TMath::Nint(0.5*hist->GetNbinsX()),0.1,5000.,hist->GetNbinsY(),hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
730   BinLogX(histBG);
731
732   TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
733   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
734   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
735   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
736   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
737   foElectron->SetParameters(ini);
738   foMuon->SetParameters(ini);
739   foPion->SetParameters(ini);
740   foKaon->SetParameters(ini);
741   foProton->SetParameters(ini);
742   
743   TCanvas *canvCheck1 = new TCanvas();
744   hist->Draw("colz");
745   foElectron->Draw("same");
746   foMuon->Draw("same");
747   foPion->Draw("same");
748   foKaon->Draw("same");  
749   foProton->Draw("same");
750  
751   // Loop over all points of the input histogram
752   
753   for(Int_t i=1; i < hist->GetNbinsX(); i++) {
754     Double_t x = hist->GetXaxis()->GetBinCenter(i);   
755     for(Int_t j=1; j < hist->GetNbinsY(); j++) {
756       Long64_t n = hist->GetBin(i, j);
757       Double_t y = hist->GetYaxis()->GetBinCenter(j);
758       Double_t entries = hist->GetBinContent(n);
759       Double_t mass = 0;
760
761       // 1. identify protons
762       mass = 0.938;
763       if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
764         for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
765       }
766
767       // 2. identify electrons
768       mass = 0.000511;
769       if (fIsCosmic) {
770         if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
771           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
772         }
773       } else {
774         if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 2*sigma && ((x>0.25&&x<0.35) || (x>1.5&&x<1.8) || (x>0.65&&x<0.7))) {
775           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
776         }
777       }
778       
779       // 3. identify either muons or pions depending on cosmic or not
780       if (fIsCosmic) {
781         mass = 0.1056;
782         if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
783           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
784         }
785       } else {
786         mass = 0.1396;
787         if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
788           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
789         }
790       }
791       
792       // 4. for pp also Kaons must be included
793       if (!fIsCosmic) {
794         mass = 0.4936;
795         if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
796           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
797         }
798       }
799     }
800   }
801
802   // Fit Aleph-Parameters to the obtained profile
803   TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
804   funcAlephD->SetParameters(ini);
805
806   TCanvas *canvCheck2 = new TCanvas();
807   histBG->Draw();
808   
809   //FitSlices
810   TObjArray * arr = new TObjArray();
811   histBG->FitSlicesY(0,0,-1,0,"QN",arr);
812   TH1D * fitMean = (TH1D*) arr->At(1);
813   fitMean->Draw("same");
814
815   funcAlephD->SetParLimits(2,1e-3,1e-7);
816   funcAlephD->SetParLimits(3,0.5,3.5);
817   funcAlephD->SetParLimits(4,0.5,3.5);
818   fitMean->Fit(funcAlephD, "QNR");
819   funcAlephD->Draw("same");
820
821   for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
822
823   foElectron->SetParameters(ini);
824   foMuon->SetParameters(ini);
825   foPion->SetParameters(ini);
826   foKaon->SetParameters(ini);
827   foProton->SetParameters(ini);
828   
829   TCanvas *canvCheck3 = new TCanvas();
830   hist->Draw("colz");
831   foElectron->Draw("same");
832   foMuon->Draw("same");
833   foPion->Draw("same");
834   foKaon->Draw("same");  
835   foProton->Draw("same");
836   
837   canvCheck1->Print();
838   canvCheck2->Print();
839   canvCheck3->Print();
840
841   return;
842
843
844 }