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 }