Changes to compile with Root6 on macosx64
[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    fMinMomentumMIP(0),
198    fMaxMomentumMIP(0),
199    fAlephParameters(),
200    fUseMax(0),
201    fLowerTrunc(0),
202    fUpperTrunc(0),
203    fUseShapeNorm(0),
204    fUsePosNorm(0),
205    fUsePadNorm(0),
206    fUseCookAnalytical(0),
207    fIsCosmic(0),
208    fLowMemoryConsumption(0)
209 {  
210   //
211   // Default constructor
212   //
213   AliDebug(5,"Default Constructor");  
214 }
215
216
217 AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain)
218   :AliTPCcalibBase(),
219    fHistGainTime(0),
220    fGainVsTime(0),
221    fHistDeDxTotal(0),
222    fIntegrationTimeDeDx(0),
223    fMIP(0),
224    fCutCrossRows(0),
225    fCutEtaWindow(0),
226    fCutRequireITSrefit(0),
227    fCutMaxDcaXY(0),
228    fCutMaxDcaZ(0),
229    fMinMomentumMIP(0),
230    fMaxMomentumMIP(0),
231    fAlephParameters(),
232    fUseMax(0),
233    fLowerTrunc(0),
234    fUpperTrunc(0),
235    fUseShapeNorm(0),
236    fUsePosNorm(0),
237    fUsePadNorm(0),
238    fUseCookAnalytical(0),
239    fIsCosmic(0),
240    fLowMemoryConsumption(0)
241 {
242   //
243   // No default constructor
244   //
245   SetName(name);
246   SetTitle(title);
247   
248   AliDebug(5,"Non Default Constructor");
249   
250   fIntegrationTimeDeDx = deltaIntegrationTimeGain;
251   
252   Double_t deltaTime = EndTime - StartTime;
253   
254   
255   // 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
256   Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
257   Int_t binsGainTime[7]    = {150,  static_cast<Int_t>(timeBins),    4,  25, 200, 10000000, 20};
258   Double_t xminGainTime[7] = {0.5, static_cast<Double_t>(StartTime),  0.5,   0, 0.1,    -0.5,  -1};
259   Double_t xmaxGainTime[7] = {  8,   static_cast<Double_t>(EndTime),  4.5, 250,  50, 9999999.5, 1};
260   fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number, eta;dEdx",7,binsGainTime,xminGainTime,xmaxGainTime);
261   BinLogX(fHistGainTime, 4);
262   //
263   fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,8);
264   BinLogX(fHistDeDxTotal);
265   
266   
267   // default track selection cuts
268   fCutCrossRows = 80;
269   fCutEtaWindow = 0.8;
270   fCutRequireITSrefit = kFALSE;
271   fCutMaxDcaXY = 3.5;
272   fCutMaxDcaZ  = 25.;
273
274   // default values for MIP window selection
275   fMinMomentumMIP = 0.4;
276   fMaxMomentumMIP = 0.6;
277   fAlephParameters[0] = 0.07657; // the following parameters work for most of the periods and are therefore default
278   fAlephParameters[1] = 10.6654; // but they can be overwritten in the train setup of cpass0
279   fAlephParameters[2] = 2.51466e-14;
280   fAlephParameters[3] = 2.05379;
281   fAlephParameters[4] = 1.84288;
282
283   // default values for dE/dx
284   fMIP = 50.;
285   fUseMax = kTRUE;
286   fLowerTrunc = 0.02;
287   fUpperTrunc = 0.6;
288   fUseShapeNorm = kTRUE;
289   fUsePosNorm = kFALSE;
290   fUsePadNorm = kFALSE;
291   fUseCookAnalytical = kFALSE;
292   //
293   fIsCosmic = kTRUE;
294   fLowMemoryConsumption = kTRUE;
295   //
296   
297 }
298
299
300
301 AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
302   //
303   // Destructor
304   //
305   delete fHistGainTime;
306   delete fGainVsTime;
307   delete fHistDeDxTotal;
308 }
309
310
311 void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
312   //
313   // main track loop
314   //
315   if (!event) {
316     Printf("ERROR: ESD not available");
317     return;
318   }
319   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
320   if (!ESDfriend) {
321    return;
322   }
323   if (ESDfriend->TestSkipBit()) return;
324
325   if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
326     ProcessCosmicEvent(event);
327   } else {
328     ProcessBeamEvent(event);
329   }
330   
331
332   
333   
334 }
335
336
337 void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
338   //
339   // Process in case of cosmic event
340   //
341   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
342   if (!esdFriend) {
343    Printf("ERROR: ESDfriend not available");
344    return;
345   }
346   //
347   UInt_t time = event->GetTimeStamp();
348   Int_t nFriendTracks = esdFriend->GetNumberOfTracks();
349   Int_t runNumber = event->GetRunNumber();
350   //
351   // track loop
352   //
353   for (Int_t i=0;i<nFriendTracks;++i) {
354
355     AliESDtrack *track = event->GetTrack(i);
356     if (!track) continue;
357     AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
358     if (!friendTrack) continue;        
359     const AliExternalTrackParam * trackIn = track->GetInnerParam();
360     const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
361     if (!trackIn) continue;
362     if (!trackOut) continue;
363
364     // calculate necessary track parameters
365     Double_t meanP = trackIn->GetP();
366     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
367     Int_t nclsDeDx = track->GetTPCNcls();
368
369     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
370     if (nclsDeDx < 60) continue;     
371     if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
372     if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
373     
374     // Get seeds
375     TObject *calibObject;
376     AliTPCseed *seed = 0;
377     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
378       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
379     }    
380
381     if (seed) { 
382       Double_t tpcSignal = GetTPCdEdx(seed);
383       if (nclsDeDx > 100) fHistDeDxTotal->Fill(meanP, tpcSignal);
384       //
385       if (fLowMemoryConsumption) {
386         if (meanP < 20) continue;
387         meanP = 30; // set all momenta to one in order to save memory
388       }
389       //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
390       Double_t vec[6] = {tpcSignal,static_cast<Double_t>(time),1,meanDrift,meanP,static_cast<Double_t>(runNumber)};
391       fHistGainTime->Fill(vec);
392     }
393     
394   }
395
396 }
397
398
399
400 void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
401   //
402   // Process in case of beam event
403   //
404   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
405   if (!esdFriend) {
406    Printf("ERROR: ESDfriend not available");
407    return;
408   }
409   //
410   UInt_t time = event->GetTimeStamp();
411   Int_t nFriendTracks = esdFriend->GetNumberOfTracks();
412   Int_t runNumber = event->GetRunNumber();
413   //
414   // track loop
415   //
416   for (Int_t i=0;i<nFriendTracks;++i) { // begin track loop
417
418     AliESDtrack *track = event->GetTrack(i);
419     if (!track) continue;
420     AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
421     if (!friendTrack) continue;
422         
423     const AliExternalTrackParam * trackIn = track->GetInnerParam();
424     const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
425     if (!trackIn) continue;
426     if (!trackOut) continue;
427
428     // calculate necessary track parameters
429     Double_t meanP = trackIn->GetP();
430     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
431     Int_t nCrossedRows = track->GetTPCCrossedRows();
432
433     //
434     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
435     // fCutCrossRows = 80, fCutEtaWindow = 0.8, fCutRequireITSrefit, fCutMaxDcaXY = 3.5, fCutMaxDcaZ = 25
436     //
437     if (nCrossedRows < fCutCrossRows) continue;     
438     if (TMath::Abs(trackIn->Eta()) > fCutEtaWindow) continue;
439     //
440     UInt_t status = track->GetStatus();
441     if ((status&AliESDtrack::kTPCrefit)==0) continue;
442     if ((status&AliESDtrack::kITSrefit)==0 && fCutRequireITSrefit) continue; // ITS cluster
443     //
444     Float_t dca[2], cov[3];
445     track->GetImpactParameters(dca,cov);
446     if (TMath::Abs(dca[0]) > fCutMaxDcaXY || TMath::Abs(dca[0]) < 0.0000001) continue;  // cut in xy
447     if (((status&AliESDtrack::kITSrefit) == 1 && TMath::Abs(dca[1]) > 3.) || TMath::Abs(dca[1]) > fCutMaxDcaZ ) continue;
448     //
449     Double_t eta = trackIn->Eta();
450     
451     // Get seeds
452     TObject *calibObject;
453     AliTPCseed *seed = 0;
454     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
455       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
456     }    
457
458     if (seed) {
459       Int_t particleCase = 0;
460       if (meanP < fMaxMomentumMIP && meanP > fMinMomentumMIP)  particleCase = 2; // MIP pions
461       if (meanP < 0.57 && meanP > 0.56) particleCase = 3; // protons 1
462       if (meanP < 0.66 && meanP > 0.65) particleCase = 4; // protons 2
463       //
464       if (fLowMemoryConsumption && particleCase == 0) continue;
465       //
466       Double_t tpcSignal = GetTPCdEdx(seed);
467       //
468       // flattens signal in MIP window
469       //
470       if (particleCase == 2) {
471         Float_t corrFactor = AliExternalTrackParam::BetheBlochAleph(meanP/0.13957, 
472                                                                     fAlephParameters[0], 
473                                                                     fAlephParameters[1], 
474                                                                     fAlephParameters[2], 
475                                                                     fAlephParameters[3],
476                                                                     fAlephParameters[4]);
477         tpcSignal /= corrFactor; 
478       } 
479       fHistDeDxTotal->Fill(meanP, tpcSignal);
480       //
481       //dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 protons), momenta, runNumner, eta
482       Double_t vec[7] = {tpcSignal,static_cast<Double_t>(time),static_cast<Double_t>(particleCase),meanDrift,meanP,static_cast<Double_t>(runNumber), eta};
483       fHistGainTime->Fill(vec);
484     }
485     
486   } // end track loop
487   //
488   // V0 loop -- in beam events the cosmic part of the histogram is filled with GammaConversions
489   //
490   for(Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) {
491     AliESDv0 * v0 = event->GetV0(iv0);
492     if (!v0->GetOnFlyStatus()) continue;
493     if (v0->GetEffMass(0,0) > 0.02) continue; // select low inv. mass
494     Double_t xyz[3];
495     v0->GetXYZ(xyz[0], xyz[1], xyz[2]);
496     if (TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) < 3 || TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) > 30) continue;
497     //
498     // "loop over daughters" 
499     //
500     for(Int_t idaughter = 0; idaughter < 2; idaughter++) { // daughter loop
501       Int_t index = idaughter == 0 ? v0->GetPindex() : v0->GetNindex();
502       AliESDtrack * trackP = event->GetTrack(index);
503       AliESDfriendTrack *friendTrackP = esdFriend->GetTrack(index);
504       if (!friendTrackP) continue;
505       const AliExternalTrackParam * trackPIn = trackP->GetInnerParam();
506       const AliExternalTrackParam * trackPOut = friendTrackP->GetTPCOut();
507       if (!trackPIn) continue;
508       if (!trackPOut) continue;
509       // calculate necessary track parameters
510       Double_t meanP = trackPIn->GetP();
511       Double_t meanDrift = 250 - 0.5*TMath::Abs(trackPIn->GetZ() + trackPOut->GetZ());
512       Int_t nclsDeDx = trackP->GetTPCNcls();
513       // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
514       if (nclsDeDx < 60) continue;     
515       if (TMath::Abs(trackPIn->GetTgl()) > 1) continue;
516       //
517       TObject *calibObject;
518       AliTPCseed *seed = 0;
519       for (Int_t l=0;(calibObject=friendTrackP->GetCalibObject(l));++l) {
520       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
521       }    
522       if (seed) { 
523         if (fLowMemoryConsumption) {
524           if (meanP > fMaxMomentumMIP || meanP < fMinMomentumMIP) continue;
525           meanP = 0.45; // set all momenta to one in order to save memory
526       }
527         Double_t tpcSignal = GetTPCdEdx(seed);
528         //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
529         Double_t vec[6] = {tpcSignal,static_cast<Double_t>(time),1,meanDrift,meanP,static_cast<Double_t>(runNumber)};
530         fHistGainTime->Fill(vec);
531       }
532     }
533     
534   }
535
536 }
537
538
539 Float_t AliTPCcalibTimeGain::GetTPCdEdx(AliTPCseed * seed) {
540   //
541   // calculate tpc dEdx
542   //
543   Double_t signal = 0;
544   //
545   if (!fUseCookAnalytical) {
546     signal = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,fUseMax,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
547   } else {
548     signal = (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax);
549   }
550   //
551   return signal;
552 }
553
554
555 void AliTPCcalibTimeGain::AnalyzeRun(Int_t minEntries) {
556   //
557   // Analyze results of calibration
558   //
559   if (fIsCosmic) {
560     fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
561     fHistGainTime->GetAxis(4)->SetRangeUser(20,100);    // only Fermi-Plateau muons
562   } else {
563     fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
564     fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
565   }
566   //
567   fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,minEntries,10);
568   //
569   return;
570 }
571
572
573 TGraphErrors * AliTPCcalibTimeGain::GetGraphGainVsTime(Int_t runNumber, Int_t minEntries) {
574   //
575   // Analyze results and get the graph 
576   //
577   if (runNumber == 0) {
578     if (!fGainVsTime) {
579       AnalyzeRun(minEntries);
580     }
581   } else {
582     // 1st check if the current run was cosmic or beam event
583     fHistGainTime->GetAxis(5)->SetRangeUser(runNumber,runNumber);
584     AnalyzeRun(minEntries);
585   }
586   if (fGainVsTime->GetN() == 0) return 0;
587   return fGainVsTime;
588 }
589
590 Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
591   //
592   // merge component
593   //
594   TIterator* iter = li->MakeIterator();
595   AliTPCcalibTimeGain* cal = 0;
596
597   while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
598     if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
599       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
600       return -1;
601     }
602
603     // add histograms here...
604     if (cal->GetHistGainTime() && fHistGainTime ) 
605     {
606       if ((fHistGainTime->GetEntries() + cal->GetHistGainTime()->GetEntries()) < fgMergeEntriesCut) 
607       {
608         //AliInfo(Form("fHistGainTime has %.0f tracks, going to add %.0f\n",fHistGainTime->GetEntries(),cal->GetHistGainTime()->GetEntries()));
609         fHistGainTime->Add(cal->GetHistGainTime());
610       }
611       else
612       {
613         AliInfo(Form("fHistGainTime full (has %.0f merged tracks, max allowed: %.0f)",fHistGainTime->GetEntries(),fgMergeEntriesCut));
614       }
615     }
616
617     if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
618
619   }
620   delete iter;
621   return 0;
622   
623 }
624
625
626 AliSplineFit * AliTPCcalibTimeGain::MakeSplineFit(TGraphErrors * graph) {
627   //
628   // make spline fit of gain
629   //
630   AliSplineFit *fit = new AliSplineFit();
631   fit->SetGraph(graph);
632   fit->SetMinPoints(graph->GetN()+1);
633   fit->InitKnots(graph,2,0,0.001);
634   fit->SplineFit(0);
635   return fit;
636   
637 }
638
639
640
641 TGraphErrors * AliTPCcalibTimeGain::GetGraphAttachment(Int_t minEntries, Int_t nmaxBin, Float_t /*fracLow*/, Float_t /*fracUp*/) {
642   //
643   // For each time bin the driftlength-dependence of the signal is fitted.
644   //
645   TH3D * hist = fHistGainTime->Projection(1, 0, 3);
646   Double_t *xvec = new Double_t[hist->GetNbinsX()];
647   Double_t *yvec = new Double_t[hist->GetNbinsX()];
648   Double_t *xerr = new Double_t[hist->GetNbinsX()];
649   Double_t *yerr = new Double_t[hist->GetNbinsX()];
650   Int_t counter  = 0;
651   TH2D * projectionHist = 0x0;
652   //
653   for(Int_t i=1; i < hist->GetNbinsX(); i++) {
654     Int_t nsum=0;
655     Int_t imin   =  i;
656     Int_t imax   =  i;    
657     for (Int_t idelta=0; idelta<nmaxBin; idelta++){
658       //
659       imin   =  TMath::Max(i-idelta,1);
660       imax   =  TMath::Min(i+idelta,hist->GetNbinsX());
661       nsum = TMath::Nint(hist->Integral(imin,imax,1,hist->GetNbinsY()-1,1,hist->GetNbinsZ()-1));
662       if (nsum==0) break;
663       if (nsum>minEntries) break;
664     }
665     if (nsum<minEntries) continue;
666     //
667     hist->GetXaxis()->SetRange(imin,imax);
668     projectionHist = (TH2D*)hist->Project3D("yzNUFNOF");
669     //
670     TObjArray arr;
671     projectionHist->FitSlicesY(0,2, projectionHist->GetNbinsX()-2,0,"QNR",&arr);
672     TH1D * histAttach = (TH1D*)arr.At(1);
673     TF1 pol1("polynom1","pol1",10,240);
674     histAttach->Fit(&pol1,"QNR");
675     xvec[counter] = 0.5*(hist->GetXaxis()->GetBinCenter(imin)+hist->GetXaxis()->GetBinCenter(imax));
676     yvec[counter] = pol1.GetParameter(1)/pol1.GetParameter(0);
677     xerr[counter] = 0;
678     yerr[counter] = pol1.GetParError(1)/pol1.GetParameter(0);
679     counter++;
680     //
681     delete projectionHist;
682   }
683   
684   TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
685   delete [] xvec;
686   delete [] yvec;
687   delete [] xerr;
688   delete [] yerr;
689   delete hist;
690   return graphErrors;
691   
692 }
693
694
695
696 void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
697   //
698   // Method for the correct logarithmic binning of histograms
699   //
700   TAxis *axis = h->GetAxis(axisDim);
701   int bins = axis->GetNbins();
702
703   Double_t from = axis->GetXmin();
704   Double_t to = axis->GetXmax();
705   Double_t *newBins = new Double_t[bins + 1];
706    
707   newBins[0] = from;
708   Double_t factor = pow(to/from, 1./bins);
709   
710   for (int i = 1; i <= bins; i++) {
711    newBins[i] = factor * newBins[i-1];
712   }
713   axis->Set(bins, newBins);
714   delete[] newBins;
715   
716 }
717
718
719 void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
720   //
721   // Method for the correct logarithmic binning of histograms
722   //
723   TAxis *axis = h->GetXaxis();
724   int bins = axis->GetNbins();
725
726   Double_t from = axis->GetXmin();
727   Double_t to = axis->GetXmax();
728   Double_t *newBins = new Double_t[bins + 1];
729    
730   newBins[0] = from;
731   Double_t factor = pow(to/from, 1./bins);
732   
733   for (int i = 1; i <= bins; i++) {
734    newBins[i] = factor * newBins[i-1];
735   }
736   axis->Set(bins, newBins);
737   delete[] newBins;
738   
739 }
740
741
742
743 void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
744   //
745   // Fit the bethe bloch params
746   //
747   //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
748   const Double_t sigma = 0.06;
749
750   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());
751   BinLogX(histBG);
752
753   TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
754   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
755   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
756   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
757   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
758   foElectron->SetParameters(ini);
759   foMuon->SetParameters(ini);
760   foPion->SetParameters(ini);
761   foKaon->SetParameters(ini);
762   foProton->SetParameters(ini);
763   
764   TCanvas *canvCheck1 = new TCanvas();
765   hist->Draw("colz");
766   foElectron->Draw("same");
767   foMuon->Draw("same");
768   foPion->Draw("same");
769   foKaon->Draw("same");  
770   foProton->Draw("same");
771  
772   // Loop over all points of the input histogram
773   
774   for(Int_t i=1; i < hist->GetNbinsX(); i++) {
775     Double_t x = hist->GetXaxis()->GetBinCenter(i);   
776     for(Int_t j=1; j < hist->GetNbinsY(); j++) {
777       Long64_t n = hist->GetBin(i, j);
778       Double_t y = hist->GetYaxis()->GetBinCenter(j);
779       Double_t entries = hist->GetBinContent(n);
780       Double_t mass = 0;
781
782       // 1. identify protons
783       mass = 0.938;
784       if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
785         for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
786       }
787
788       // 2. identify electrons
789       mass = 0.000511;
790       if (fIsCosmic) {
791         if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
792           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
793         }
794       } else {
795         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))) {
796           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
797         }
798       }
799       
800       // 3. identify either muons or pions depending on cosmic or not
801       if (fIsCosmic) {
802         mass = 0.1056;
803         if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
804           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
805         }
806       } else {
807         mass = 0.1396;
808         if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
809           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
810         }
811       }
812       
813       // 4. for pp also Kaons must be included
814       if (!fIsCosmic) {
815         mass = 0.4936;
816         if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
817           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
818         }
819       }
820     }
821   }
822
823   // Fit Aleph-Parameters to the obtained profile
824   TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
825   funcAlephD->SetParameters(ini);
826
827   TCanvas *canvCheck2 = new TCanvas();
828   histBG->Draw();
829   
830   //FitSlices
831   TObjArray * arr = new TObjArray();
832   histBG->FitSlicesY(0,0,-1,0,"QN",arr);
833   TH1D * fitMean = (TH1D*) arr->At(1);
834   fitMean->Draw("same");
835
836   funcAlephD->SetParLimits(2,1e-3,1e-7);
837   funcAlephD->SetParLimits(3,0.5,3.5);
838   funcAlephD->SetParLimits(4,0.5,3.5);
839   fitMean->Fit(funcAlephD, "QNR");
840   funcAlephD->Draw("same");
841
842   for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
843
844   foElectron->SetParameters(ini);
845   foMuon->SetParameters(ini);
846   foPion->SetParameters(ini);
847   foKaon->SetParameters(ini);
848   foProton->SetParameters(ini);
849   
850   TCanvas *canvCheck3 = new TCanvas();
851   hist->Draw("colz");
852   foElectron->Draw("same");
853   foMuon->Draw("same");
854   foPion->Draw("same");
855   foKaon->Draw("same");  
856   foProton->Draw("same");
857   
858   canvCheck1->Print();
859   canvCheck2->Print();
860   canvCheck3->Print();
861
862   return;
863
864
865 }