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     AliExternalTrackParam trckIn;
 373     track->GetTrackParamIp(trckIn);
 374     if ( (track->GetTrackParamIp(trckIn)) < 0) continue;
 375     AliExternalTrackParam * trackIn = &trckIn;
 376     if (!trackIn) continue;
 377
 378     //const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
 379     AliExternalTrackParam trckOut;
 380     friendTrack->GetTrackParamTPCOut(trckOut);
 381     if ( (friendTrack->GetTrackParamTPCOut(trckOut)) < 0) continue;
 382     AliExternalTrackParam * trackOut = &trckOut;
 383     if (!trackOut) continue;
 384
 385     // calculate necessary track parameters
 386     Double_t meanP = trackIn->GetP();
 387     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
 388     Int_t nclsDeDx = track->GetTPCNcls();
 389
 390     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
 391     if (nclsDeDx < 60) continue;
 392     if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
 393     if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
 394
 395     // Get seeds
 396     TObject *calibObject;
 397     AliTPCseed *seed = 0;
 398     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 399       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 400     }
 401
 402     if (seed) {
 403       Double_t tpcSignal = GetTPCdEdx(seed);
 404       if (nclsDeDx > 100) fHistDeDxTotal->Fill(meanP, tpcSignal);
 405       //
 406       if (fLowMemoryConsumption) {
 407         if (meanP < 20) continue;
 408         meanP = 30; // set all momenta to one in order to save memory
 409       }
 410       //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
 411       Double_t vec[6] = {tpcSignal,static_cast<Double_t>(time),1,meanDrift,meanP,static_cast<Double_t>(runNumber)};
 412       fHistGainTime->Fill(vec);
 413     }
 414
 415   }
 416
 417 }
 418
 419
 420
 421 void AliTPCcalibTimeGain::ProcessBeamEvent(AliVEvent *event) {
 422   //
 423   // Process in case of beam event
 424   //
 425     //Printf("AliTPCcalibTimeGain::ProcessBeamEvent(event)...");
 426
 427   AliVfriendEvent *friendEvent=event->FindFriend();
 428   if (!friendEvent) {
 429    //Printf("ERROR AliTPCcalibTimeGain::ProcessBeamEvent(): ESDfriend not available");
 430    return;
 431   }
 432   //
 433   UInt_t time = event->GetTimeStamp();
 434   if (!time) Printf("ERROR: no time stamp available!");
 435   Int_t nFriendTracks = friendEvent->GetNumberOfTracks();
 436   Int_t runNumber = event->GetRunNumber();
 437   //
 438   // track loop
 439   //
 440   for (Int_t i=0;i<nFriendTracks;++i) { // begin track loop
 441
 442     AliVTrack *track = event->GetVTrack(i);
 443     if (!track) {
 444         //Printf("***ERROR*** : track not available");
 445         continue;}
 446     const AliVfriendTrack *friendTrack = friendEvent->GetTrack(i);
 447     if (!friendTrack) {
 448         //Printf("ERROR ProcessBeamEvent(): friendTrack is not available!");
 449         continue;
 450     }
 451
 452     //const AliExternalTrackParam * trackIn = track->GetInnerParam();
 453     AliExternalTrackParam trckIn;
 454     track->GetTrackParamIp(trckIn);
 455     if ( (track->GetTrackParamIp(trckIn)) < 0) continue;
 456     AliExternalTrackParam * trackIn = &trckIn;
 457     if (!trackIn) continue;
 458
 459     //const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
 460     AliExternalTrackParam trckOut;
 461     friendTrack->GetTrackParamTPCOut(trckOut);
 462     if ( (friendTrack->GetTrackParamTPCOut(trckOut)) < 0) continue;
 463     AliExternalTrackParam * trackOut = &trckOut;
 464     if (!trackOut) continue;
 465
 466     // calculate necessary track parameters
 467     Double_t meanP = trackIn->GetP();
 468     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
 469     Int_t nCrossedRows = track->GetTPCCrossedRows();
 470
 471     //
 472     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
 473     // fCutCrossRows = 80, fCutEtaWindow = 0.8, fCutRequireITSrefit, fCutMaxDcaXY = 3.5, fCutMaxDcaZ = 25
 474     //
 475     if (nCrossedRows < fCutCrossRows) continue;
 476     if (TMath::Abs(trackIn->Eta()) > fCutEtaWindow) continue;
 477     //
 478     UInt_t status = track->GetStatus();
 479     if ((status&AliVTrack::kTPCrefit)==0) continue;
 480     if ((status&AliVTrack::kITSrefit)==0 && fCutRequireITSrefit) continue; // ITS cluster
 481     //
 482     Float_t dca[2], cov[3];
 483     track->GetImpactParameters(dca,cov);
 484     if (TMath::Abs(dca[0]) > fCutMaxDcaXY || TMath::Abs(dca[0]) < 0.0000001) continue;  // cut in xy
 485     if (((status&AliVTrack::kITSrefit) == 1 && TMath::Abs(dca[1]) > 3.) || TMath::Abs(dca[1]) > fCutMaxDcaZ ) continue;
 486     //
 487     Double_t eta = trackIn->Eta();
 488
 489     // Get seeds
 490     TObject *calibObject;
 491     AliTPCseed *seed = 0;
 492     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 493       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 494     }
 495
 496     if (seed) {
 497       Int_t particleCase = 0;
 498       if (meanP < fMaxMomentumMIP && meanP > fMinMomentumMIP)  particleCase = 2; // MIP pions
 499       if (meanP < 0.57 && meanP > 0.56) particleCase = 3; // protons 1
 500       if (meanP < 0.66 && meanP > 0.65) particleCase = 4; // protons 2
 501       //
 502       if (fLowMemoryConsumption && particleCase == 0) continue;
 503       //
 504       Double_t tpcSignal = GetTPCdEdx(seed);
 505       //
 506       // flattens signal in MIP window
 507       //
 508       if (particleCase == 2) {
 509         Float_t corrFactor = AliExternalTrackParam::BetheBlochAleph(meanP/0.13957,
 510                                                                     fAlephParameters[0],
 511                                                                     fAlephParameters[1],
 512                                                                     fAlephParameters[2],
 513                                                                     fAlephParameters[3],
 514                                                                     fAlephParameters[4]);
 515         tpcSignal /= corrFactor;
 516       }
 517       //Printf("Fill DeDx histo..");
 518       fHistDeDxTotal->Fill(meanP, tpcSignal);
 519       //
 520       //dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 protons), momenta, runNumner, eta
 521       Double_t vec[7] = {tpcSignal,static_cast<Double_t>(time),static_cast<Double_t>(particleCase),meanDrift,meanP,static_cast<Double_t>(runNumber), eta};
 522       //Printf("Fill Gain histo in track loop...");
 523       fHistGainTime->Fill(vec);
 524     }
 525
 526   } // end track loop
 527   //
 528   // V0 loop -- in beam events the cosmic part of the histogram is filled with GammaConversions
 529   //
 530   for(Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) {
 531     //AliESDv0 * v0 = event->GetV0(iv0);
 532      AliESDv0 v0dummy;
 533      event->GetV0(v0dummy, iv0);
 534      AliESDv0 *v0 = &v0dummy;
 535
 536      //if (!v0) Printf("ERROR AliTPCcalibTimeGain::ProcessBeamEvent(): ESDv0 not available! ");
 537
 538     if (!v0->GetOnFlyStatus()) continue;
 539     if (v0->GetEffMass(0,0) > 0.02) continue; // select low inv. mass
 540     Double_t xyz[3];
 541     v0->GetXYZ(xyz[0], xyz[1], xyz[2]);
 542     if (TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) < 3 || TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) > 30) continue;
 543     //
 544     // "loop over daughters"
 545     //
 546     for(Int_t idaughter = 0; idaughter < 2; idaughter++) { // daughter loop
 547       Int_t index = idaughter == 0 ? v0->GetPindex() : v0->GetNindex();
 548       AliVTrack * trackP = event->GetVTrack(index);
 549       if (!trackP) continue; //Printf("***ERROR*** trackP not available!");
 550       const AliVfriendTrack *friendTrackP = friendEvent->GetTrack(index);
 551       if (!friendTrackP) continue;
 552       //const AliExternalTrackParam * trackPIn = trackP->GetInnerParam();
 553       AliExternalTrackParam trckPIn;
 554       trackP->GetTrackParamIp(trckPIn);
 555       if ( (trackP->GetTrackParamIp(trckPIn)) < 0) continue;
 556       AliExternalTrackParam * trackPIn = &trckPIn;
 557       if (!trackPIn) continue;
 558
 559       //const AliExternalTrackParam * trackPOut = friendTrackP->GetTPCOut();
 560       AliExternalTrackParam trckPOut;
 561       friendTrackP->GetTrackParamTPCOut(trckPOut);
 562       if ( (friendTrackP->GetTrackParamTPCOut(trckPOut)) < 0) continue;
 563       AliExternalTrackParam * trackPOut = &trckPOut;
 564       if (!trackPOut) continue;
 565
 566       // calculate necessary track parameters
 567       Double_t meanP = trackPIn->GetP();
 568       Double_t meanDrift = 250 - 0.5*TMath::Abs(trackPIn->GetZ() + trackPOut->GetZ());
 569       Int_t nclsDeDx = trackP->GetTPCNcls();
 570       // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
 571       if (nclsDeDx < 60) continue;
 572       if (TMath::Abs(trackPIn->GetTgl()) > 1) continue;
 573       //
 574       TObject *calibObject;
 575       AliTPCseed *seed = 0;
 576       for (Int_t l=0;(calibObject=friendTrackP->GetCalibObject(l));++l) {
 577       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 578       }
 579       if (seed) {
 580         if (fLowMemoryConsumption) {
 581           if (meanP > fMaxMomentumMIP || meanP < fMinMomentumMIP) continue;
 582           meanP = 0.45; // set all momenta to one in order to save memory
 583       }
 584         Double_t tpcSignal = GetTPCdEdx(seed);
 585         //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
 586         Double_t vec[6] = {tpcSignal,static_cast<Double_t>(time),1,meanDrift,meanP,static_cast<Double_t>(runNumber)};
 587     //Printf("Fill Gain histo in v0 loop...");
 588     fHistGainTime->Fill(vec);
 589       }
 590     }
 591
 592   }
 593
 594 }
 595
 596
 597 Float_t AliTPCcalibTimeGain::GetTPCdEdx(AliTPCseed * seed) {
 598   //
 599   // calculate tpc dEdx
 600   //
 601   Double_t signal = 0;
 602   //
 603   if (!fUseCookAnalytical) {
 604     signal = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,fUseMax,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 605   } else {
 606     signal = (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax);
 607   }
 608   //
 609   return signal;
 610 }
 611
 612
 613 void AliTPCcalibTimeGain::AnalyzeRun(Int_t minEntries) {
 614   //
 615   // Analyze results of calibration
 616   //
 617   if (fIsCosmic) {
 618     fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
 619     fHistGainTime->GetAxis(4)->SetRangeUser(20,100);    // only Fermi-Plateau muons
 620   } else {
 621     fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
 622     fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
 623   }
 624   //
 625   fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,minEntries,10);
 626   //
 627   return;
 628 }
 629
 630
 631 TGraphErrors * AliTPCcalibTimeGain::GetGraphGainVsTime(Int_t runNumber, Int_t minEntries) {
 632   //
 633   // Analyze results and get the graph
 634   //
 635   if (runNumber == 0) {
 636     if (!fGainVsTime) {
 637       AnalyzeRun(minEntries);
 638     }
 639   } else {
 640     // 1st check if the current run was cosmic or beam event
 641     fHistGainTime->GetAxis(5)->SetRangeUser(runNumber,runNumber);
 642     AnalyzeRun(minEntries);
 643   }
 644   if (fGainVsTime->GetN() == 0) return 0;
 645   return fGainVsTime;
 646 }
 647
 648 Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
 649   //
 650   // merge component
 651   //
 652   TIterator* iter = li->MakeIterator();
 653   AliTPCcalibTimeGain* cal = 0;
 654
 655   while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
 656     if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
 657       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 658       return -1;
 659     }
 660
 661     // add histograms here...
 662     if (cal->GetHistGainTime() && fHistGainTime )
 663     {
 664       if ((fHistGainTime->GetEntries() + cal->GetHistGainTime()->GetEntries()) < fgMergeEntriesCut)
 665       {
 666         //AliInfo(Form("fHistGainTime has %.0f tracks, going to add %.0f\n",fHistGainTime->GetEntries(),cal->GetHistGainTime()->GetEntries()));
 667         fHistGainTime->Add(cal->GetHistGainTime());
 668       }
 669       else
 670       {
 671         AliInfo(Form("fHistGainTime full (has %.0f merged tracks, max allowed: %.0f)",fHistGainTime->GetEntries(),fgMergeEntriesCut));
 672       }
 673     }
 674
 675     if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
 676
 677   }
 678   delete iter;
 679   return 0;
 680
 681 }
 682
 683
 684 AliSplineFit * AliTPCcalibTimeGain::MakeSplineFit(TGraphErrors * graph) {
 685   //
 686   // make spline fit of gain
 687   //
 688   AliSplineFit *fit = new AliSplineFit();
 689   fit->SetGraph(graph);
 690   fit->SetMinPoints(graph->GetN()+1);
 691   fit->InitKnots(graph,2,0,0.001);
 692   fit->SplineFit(0);
 693   return fit;
 694
 695 }
 696
 697
 698
 699 TGraphErrors * AliTPCcalibTimeGain::GetGraphAttachment(Int_t minEntries, Int_t nmaxBin, Float_t /*fracLow*/, Float_t /*fracUp*/) {
 700   //
 701   // For each time bin the driftlength-dependence of the signal is fitted.
 702   //
 703   TH3D * hist = fHistGainTime->Projection(1, 0, 3);
 704   Double_t *xvec = new Double_t[hist->GetNbinsX()];
 705   Double_t *yvec = new Double_t[hist->GetNbinsX()];
 706   Double_t *xerr = new Double_t[hist->GetNbinsX()];
 707   Double_t *yerr = new Double_t[hist->GetNbinsX()];
 708   Int_t counter  = 0;
 709   TH2D * projectionHist = 0x0;
 710   //
 711   for(Int_t i=1; i < hist->GetNbinsX(); i++) {
 712     Int_t nsum=0;
 713     Int_t imin   =  i;
 714     Int_t imax   =  i;
 715     for (Int_t idelta=0; idelta<nmaxBin; idelta++){
 716       //
 717       imin   =  TMath::Max(i-idelta,1);
 718       imax   =  TMath::Min(i+idelta,hist->GetNbinsX());
 719       nsum = TMath::Nint(hist->Integral(imin,imax,1,hist->GetNbinsY()-1,1,hist->GetNbinsZ()-1));
 720       if (nsum==0) break;
 721       if (nsum>minEntries) break;
 722     }
 723     if (nsum<minEntries) continue;
 724     //
 725     hist->GetXaxis()->SetRange(imin,imax);
 726     projectionHist = (TH2D*)hist->Project3D("yzNUFNOF");
 727     //
 728     TObjArray arr;
 729     projectionHist->FitSlicesY(0,2, projectionHist->GetNbinsX()-2,0,"QNR",&arr);
 730     TH1D * histAttach = (TH1D*)arr.At(1);
 731     TF1 pol1("polynom1","pol1",10,240);
 732     histAttach->Fit(&pol1,"QNR");
 733     xvec[counter] = 0.5*(hist->GetXaxis()->GetBinCenter(imin)+hist->GetXaxis()->GetBinCenter(imax));
 734     yvec[counter] = pol1.GetParameter(1)/pol1.GetParameter(0);
 735     xerr[counter] = 0;
 736     yerr[counter] = pol1.GetParError(1)/pol1.GetParameter(0);
 737     counter++;
 738     //
 739     delete projectionHist;
 740   }
 741
 742   TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
 743   delete [] xvec;
 744   delete [] yvec;
 745   delete [] xerr;
 746   delete [] yerr;
 747   delete hist;
 748   return graphErrors;
 749
 750 }
 751
 752
 753
 754 void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
 755   //
 756   // Method for the correct logarithmic binning of histograms
 757   //
 758   TAxis *axis = h->GetAxis(axisDim);
 759   int bins = axis->GetNbins();
 760
 761   Double_t from = axis->GetXmin();
 762   Double_t to = axis->GetXmax();
 763   Double_t *newBins = new Double_t[bins + 1];
 764
 765   newBins[0] = from;
 766   Double_t factor = pow(to/from, 1./bins);
 767
 768   for (int i = 1; i <= bins; i++) {
 769    newBins[i] = factor * newBins[i-1];
 770   }
 771   axis->Set(bins, newBins);
 772   delete[] newBins;
 773
 774 }
 775
 776
 777 void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
 778   //
 779   // Method for the correct logarithmic binning of histograms
 780   //
 781   TAxis *axis = h->GetXaxis();
 782   int bins = axis->GetNbins();
 783
 784   Double_t from = axis->GetXmin();
 785   Double_t to = axis->GetXmax();
 786   Double_t *newBins = new Double_t[bins + 1];
 787
 788   newBins[0] = from;
 789   Double_t factor = pow(to/from, 1./bins);
 790
 791   for (int i = 1; i <= bins; i++) {
 792    newBins[i] = factor * newBins[i-1];
 793   }
 794   axis->Set(bins, newBins);
 795   delete[] newBins;
 796
 797 }
 798
 799
 800
 801 void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
 802   //
 803   // Fit the bethe bloch params
 804   //
 805   //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
 806   const Double_t sigma = 0.06;
 807
 808   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());
 809   BinLogX(histBG);
 810
 811   TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
 812   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
 813   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
 814   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
 815   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
 816   foElectron->SetParameters(ini);
 817   foMuon->SetParameters(ini);
 818   foPion->SetParameters(ini);
 819   foKaon->SetParameters(ini);
 820   foProton->SetParameters(ini);
 821
 822   TCanvas *canvCheck1 = new TCanvas();
 823   hist->Draw("colz");
 824   foElectron->Draw("same");
 825   foMuon->Draw("same");
 826   foPion->Draw("same");
 827   foKaon->Draw("same");
 828   foProton->Draw("same");
 829
 830   // Loop over all points of the input histogram
 831
 832   for(Int_t i=1; i < hist->GetNbinsX(); i++) {
 833     Double_t x = hist->GetXaxis()->GetBinCenter(i);
 834     for(Int_t j=1; j < hist->GetNbinsY(); j++) {
 835       Long64_t n = hist->GetBin(i, j);
 836       Double_t y = hist->GetYaxis()->GetBinCenter(j);
 837       Double_t entries = hist->GetBinContent(n);
 838       Double_t mass = 0;
 839
 840       // 1. identify protons
 841       mass = 0.938;
 842       if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
 843         for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 844       }
 845
 846       // 2. identify electrons
 847       mass = 0.000511;
 848       if (fIsCosmic) {
 849         if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
 850           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 851         }
 852       } else {
 853         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))) {
 854           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 855         }
 856       }
 857
 858       // 3. identify either muons or pions depending on cosmic or not
 859       if (fIsCosmic) {
 860         mass = 0.1056;
 861         if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
 862           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 863         }
 864       } else {
 865         mass = 0.1396;
 866         if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
 867           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 868         }
 869       }
 870
 871       // 4. for pp also Kaons must be included
 872       if (!fIsCosmic) {
 873         mass = 0.4936;
 874         if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
 875           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 876         }
 877       }
 878     }
 879   }
 880
 881   // Fit Aleph-Parameters to the obtained profile
 882   TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
 883   funcAlephD->SetParameters(ini);
 884
 885   TCanvas *canvCheck2 = new TCanvas();
 886   histBG->Draw();
 887
 888   //FitSlices
 889   TObjArray * arr = new TObjArray();
 890   histBG->FitSlicesY(0,0,-1,0,"QN",arr);
 891   TH1D * fitMean = (TH1D*) arr->At(1);
 892   fitMean->Draw("same");
 893
 894   funcAlephD->SetParLimits(2,1e-3,1e-7);
 895   funcAlephD->SetParLimits(3,0.5,3.5);
 896   funcAlephD->SetParLimits(4,0.5,3.5);
 897   fitMean->Fit(funcAlephD, "QNR");
 898   funcAlephD->Draw("same");
 899
 900   for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
 901
 902   foElectron->SetParameters(ini);
 903   foMuon->SetParameters(ini);
 904   foPion->SetParameters(ini);
 905   foKaon->SetParameters(ini);
 906   foProton->SetParameters(ini);
 907
 908   TCanvas *canvCheck3 = new TCanvas();
 909   hist->Draw("colz");
 910   foElectron->Draw("same");
 911   foMuon->Draw("same");
 912   foPion->Draw("same");
 913   foKaon->Draw("same");
 914   foProton->Draw("same");
 915
 916   canvCheck1->Print();
 917   canvCheck2->Print();
 918   canvCheck3->Print();
 919
 920   return;
 921
 922
 923 }