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