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 }