[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTimeGain.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*

//0.  Libraries to lod
gSystem->Load("libANALYSIS");
gSystem->Load("libSTAT");
gSystem->Load("libTPCcalib");


//1. Do calibration ...
//
// compare reference

//
//2. Visualize results
//
TFile fcalib("CalibObjects.root");
TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
TGraphErrors * gr = AliTPCcalibBase::FitSlices(gain->GetHistGainTime(),0,1,2000,10,0.2,0.8);
gain->GetHistGainTime()->Projection(0,1)->Draw("colz")
gr->SetMarkerStyle(25);
gr->Draw("lp")


//
// MakeSlineFit example
//
AliSplineFit fit;
fit.SetGraph(gr)
fit->SetMinPoints(gr->GetN()+1);
fit->InitKnots(gr,2,0,0.001)
fit.SplineFit(0)
fit.MakeDiffHisto(gr)->Draw();
TGraph * grfit = fit.MakeGraph(gr->GetX()[0],gr->GetX()[gr->GetN()-1],50000,0);

gr->SetMarkerStyle(25);
gr->Draw("alp");
grfit->SetLineColor(2);
grfit->Draw("lu");

*/


#include "Riostream.h"
#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "THnSparse.h"
#include "TGraphErrors.h"
#include "TList.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TF1.h"
#include "TVectorD.h"
#include "TProfile.h"
#include "TGraphErrors.h"
#include "TCanvas.h"

#include "AliTPCclusterMI.h"
#include "AliTPCseed.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliESDfriend.h"
#include "AliESDInputHandler.h"
#include "AliAnalysisManager.h"

#include "AliTracker.h"
#include "AliMagF.h"
#include "AliTPCCalROC.h"

#include "AliLog.h"

#include "AliTPCcalibTimeGain.h"

#include "TTreeStream.h"
#include "AliTPCTracklet.h"
#include "TTimeStamp.h"
#include "AliTPCcalibDB.h"
#include "AliTPCcalibLaser.h"
#include "AliDCSSensorArray.h"
#include "AliDCSSensor.h"

ClassImp(AliTPCcalibTimeGain)


AliTPCcalibTimeGain::AliTPCcalibTimeGain() 
  :AliTPCcalibBase(), 
   fHistGainTime(0),
   fGainVsTime(0),
   fHistDeDxTotal(0),
   fIntegrationTimeDeDx(0),
   fMIP(0),
   fLowerTrunc(0),
   fUpperTrunc(0),
   fUseShapeNorm(0),
   fUsePosNorm(0),
   fUsePadNorm(0),
   fIsCosmic(0),
   fLowMemoryConsumption(0)
{  
  AliInfo("Default Constructor");  
}


AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain)
  :AliTPCcalibBase(),
   fHistGainTime(0),
   fGainVsTime(0),
   fHistDeDxTotal(0),
   fIntegrationTimeDeDx(0),
   fMIP(0),
   fLowerTrunc(0),
   fUpperTrunc(0),
   fUseShapeNorm(0),
   fUsePosNorm(0),
   fUsePadNorm(0),
   fIsCosmic(0),
   fLowMemoryConsumption(0)
 {
  
  SetName(name);
  SetTitle(title);

  AliInfo("Non Default Constructor");

  fIntegrationTimeDeDx = deltaIntegrationTimeGain;
 
  Double_t deltaTime = EndTime - StartTime;
  

  // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available)
  Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
  Int_t binsGainTime[5]    = {100,  timeBins,    2,  25, 200};
  Double_t xminGainTime[5] = {0.5, StartTime,  0.5,   0, 0.1};
  Double_t xmaxGainTime[5] = {  4,   EndTime,  2.5, 250, 50};
  fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx;dEdx,time,type,driftlength,momenta",5,binsGainTime,xminGainTime,xmaxGainTime);
  BinLogX(fHistGainTime, 4);
  //
  fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,1000);
  BinLogX(fHistDeDxTotal);
  
  // default values for dE/dx
  fMIP = 50.;
  fLowerTrunc = 0.0;
  fUpperTrunc = 0.7;
  fUseShapeNorm = kTRUE;
  fUsePosNorm = kFALSE;
  fUsePadNorm = kFALSE;
  //
  fIsCosmic = kTRUE;
  fLowMemoryConsumption = kTRUE;
  //
  
 }


AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
  //
  //
  //
}


void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
  //
  // main track loop
  //
  if (!event) {
    Printf("ERROR: ESD not available");
    return;
  }
  
  if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
    ProcessCosmicEvent(event);
  } else {
    ProcessBeamEvent(event);
  }
  

}


void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {

  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
  if (!ESDfriend) {
   Printf("ERROR: ESDfriend not available");
   return;
  }
  //
  UInt_t time = event->GetTimeStamp();
  Int_t ntracks = event->GetNumberOfTracks();
  //
  // track loop
  //
  for (Int_t i=0;i<ntracks;++i) {

    AliESDtrack *track = event->GetTrack(i);
    if (!track) continue;
        
    const AliExternalTrackParam * trackIn = track->GetInnerParam();
    const AliExternalTrackParam * trackOut = track->GetOuterParam();
    if (!trackIn) continue;
    if (!trackOut) continue;

    // calculate necessary track parameters
    Double_t meanP = trackIn->GetP();
    Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
    Int_t NclsDeDx = track->GetTPCNcls();

    // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
    if (NclsDeDx < 60) continue;     
    if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
    if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
    
    // Get seeds
    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
    if (!friendTrack) continue;
    TObject *calibObject;
    AliTPCseed *seed = 0;
    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
    }    

    if (seed) { 
      Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
      fHistDeDxTotal->Fill(meanP, TPCsignalMax);
      //
      if (fLowMemoryConsumption) {
	if (meanP < 20) continue;
	meanP = 30; // set all momenta to one in order to save memory
      }
      //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
      Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP}; 
      fHistGainTime->Fill(vec);
    }
    
  }

}


void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {

  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
  if (!ESDfriend) {
   Printf("ERROR: ESDfriend not available");
   return;
  }
  //
  UInt_t time = event->GetTimeStamp();
  Int_t ntracks = event->GetNumberOfTracks();
  //
  // track loop
  //
  for (Int_t i=0;i<ntracks;++i) {

    AliESDtrack *track = event->GetTrack(i);
    if (!track) continue;
        
    const AliExternalTrackParam * trackIn = track->GetInnerParam();
    const AliExternalTrackParam * trackOut = track->GetOuterParam();
    if (!trackIn) continue;
    if (!trackOut) continue;

    // calculate necessary track parameters
    Double_t meanP = trackIn->GetP();
    Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
    Int_t NclsDeDx = track->GetTPCNcls();

    // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
    if (NclsDeDx < 60) continue;     
    if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
    if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
    
    // Get seeds
    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
    if (!friendTrack) continue;
    TObject *calibObject;
    AliTPCseed *seed = 0;
    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
    }    

    if (seed) { 
      Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
      fHistDeDxTotal->Fill(meanP, TPCsignalMax);
      //
      if (fLowMemoryConsumption) {
	if (meanP > 0.5 || meanP < 0.4) continue;
	meanP = 0.45; // set all momenta to one in order to save memory
      }
      //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
      Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP}; 
      fHistGainTime->Fill(vec);
    }
    
  }

}


void AliTPCcalibTimeGain::Analyze() {
  //
  //
  //
  if (fIsCosmic) {
    fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
    fHistGainTime->GetAxis(4)->SetRangeUser(20,100);    // only Fermi-Plateau muons
  } else {
    fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
    fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
  }
  //
  fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,2000,10);
  //
  return;
}


Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {

  TIterator* iter = li->MakeIterator();
  AliTPCcalibTimeGain* cal = 0;

  while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
    if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
      Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
      return -1;
    }

    // add histograms here...
    if (cal->GetHistGainTime()) fHistGainTime->Add(cal->GetHistGainTime());
    if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());

  }
  
  return 0;
  
}


void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {

  // Method for the correct logarithmic binning of histograms

  TAxis *axis = h->GetAxis(axisDim);
  int bins = axis->GetNbins();

  Double_t from = axis->GetXmin();
  Double_t to = axis->GetXmax();
  Double_t *new_bins = new Double_t[bins + 1];
   
  new_bins[0] = from;
  Double_t factor = pow(to/from, 1./bins);
  
  for (int i = 1; i <= bins; i++) {
   new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete new_bins;
  
}


void AliTPCcalibTimeGain::BinLogX(TH1 *h) {

  // Method for the correct logarithmic binning of histograms

  TAxis *axis = h->GetXaxis();
  int bins = axis->GetNbins();

  Double_t from = axis->GetXmin();
  Double_t to = axis->GetXmax();
  Double_t *new_bins = new Double_t[bins + 1];
   
  new_bins[0] = from;
  Double_t factor = pow(to/from, 1./bins);
  
  for (int i = 1; i <= bins; i++) {
   new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete new_bins;
  
}


void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
  //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
  const Double_t sigma = 0.06;

  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());
  BinLogX(histBG);

  TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
  TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
  TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
  TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
  TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
  foElectron->SetParameters(ini);
  foMuon->SetParameters(ini);
  foPion->SetParameters(ini);
  foKaon->SetParameters(ini);
  foProton->SetParameters(ini);
  
  TCanvas *CanvCheck1 = new TCanvas();
  hist->Draw("colz");
  foElectron->Draw("same");
  foMuon->Draw("same");
  foPion->Draw("same");
  foKaon->Draw("same");  
  foProton->Draw("same");
 
  // Loop over all points of the input histogram
  
  for(Int_t i=1; i < hist->GetNbinsX(); i++) {
    Double_t x = hist->GetXaxis()->GetBinCenter(i);   
    for(Int_t j=1; j < hist->GetNbinsY(); j++) {
      Long64_t n = hist->GetBin(i, j);
      Double_t y = hist->GetYaxis()->GetBinCenter(j);
      Double_t entries = hist->GetBinContent(n);
      Double_t mass = 0;

      // 1. identify protons
      mass = 0.938;
      if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
	for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
      }

      // 2. identify electrons
      mass = 0.000511;
      if (fIsCosmic) {
	if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
	  for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	}
      } else {
	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))) {
	  for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	}
      }
      
      // 3. identify either muons or pions depending on cosmic or not
      if (fIsCosmic) {
	mass = 0.1056;
	if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
	  for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	}
      } else {
	mass = 0.1396;
	if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
	  for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	}
      }
      
      // 4. for pp also Kaons must be included
      if (!fIsCosmic) {
	mass = 0.4936;
	if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
	  for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	}
      }
    }
  }

  // Fit Aleph-Parameters to the obtained profile
  TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
  funcAlephD->SetParameters(ini);

  TCanvas *CanvCheck2 = new TCanvas();
  histBG->Draw();
  
  //FitSlices
  TObjArray * arr = new TObjArray();
  histBG->FitSlicesY(0,0,-1,0,"QN",arr);
  TH1D * fitMean = (TH1D*) arr->At(1);
  fitMean->Draw("same");

  funcAlephD->SetParLimits(2,1e-3,1e-7);
  funcAlephD->SetParLimits(3,0.5,3.5);
  funcAlephD->SetParLimits(4,0.5,3.5);
  fitMean->Fit(funcAlephD, "QNR");
  funcAlephD->Draw("same");

  for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);

  foElectron->SetParameters(ini);
  foMuon->SetParameters(ini);
  foPion->SetParameters(ini);
  foKaon->SetParameters(ini);
  foProton->SetParameters(ini);
  
  TCanvas *CanvCheck3 = new TCanvas();
  hist->Draw("colz");
  foElectron->Draw("same");
  foMuon->Draw("same");
  foPion->Draw("same");
  foKaon->Draw("same");  
  foProton->Draw("same");
  
  CanvCheck1->Print();
  CanvCheck2->Print();
  CanvCheck3->Print();

  return;


}
Commit	Line	Data
489dce1b	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	//0. Libraries to lod
	19	gSystem->Load("libANALYSIS");
	20	gSystem->Load("libSTAT");
	21	gSystem->Load("libTPCcalib");
	22
	23
	24	//1. Do calibration ...
	25	//
	26	// compare reference
	27
	28	//
d0b31f57	29	//2. Visualize results
489dce1b	30	//
	31	TFile fcalib("CalibObjects.root");
	32	TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
	33	AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
d0b31f57	34	TGraphErrors * gr = AliTPCcalibBase::FitSlices(gain->GetHistGainTime(),0,1,2000,10,0.2,0.8);
	35	gain->GetHistGainTime()->Projection(0,1)->Draw("colz")
	36	gr->SetMarkerStyle(25);
	37	gr->Draw("lp")
489dce1b	38
	39
	40	//
	41	// MakeSlineFit example
	42	//
	43	AliSplineFit fit;
	44	fit.SetGraph(gr)
	45	fit->SetMinPoints(gr->GetN()+1);
	46	fit->InitKnots(gr,2,0,0.001)
	47	fit.SplineFit(0)
	48	fit.MakeDiffHisto(gr)->Draw();
	49	TGraph * grfit = fit.MakeGraph(gr->GetX()[0],gr->GetX()[gr->GetN()-1],50000,0);
	50
	51	gr->SetMarkerStyle(25);
	52	gr->Draw("alp");
	53	grfit->SetLineColor(2);
	54	grfit->Draw("lu");
	55
	56	*/
	57
	58
	59	#include "Riostream.h"
	60	#include "TChain.h"
	61	#include "TTree.h"
	62	#include "TH1F.h"
	63	#include "TH2F.h"
	64	#include "TH3F.h"
	65	#include "THnSparse.h"
d0b31f57	66	#include "TGraphErrors.h"
489dce1b	67	#include "TList.h"
	68	#include "TMath.h"
	69	#include "TCanvas.h"
	70	#include "TFile.h"
	71	#include "TF1.h"
	72	#include "TVectorD.h"
	73	#include "TProfile.h"
	74	#include "TGraphErrors.h"
	75	#include "TCanvas.h"
	76
	77	#include "AliTPCclusterMI.h"
	78	#include "AliTPCseed.h"
	79	#include "AliESDVertex.h"
	80	#include "AliESDEvent.h"
	81	#include "AliESDfriend.h"
	82	#include "AliESDInputHandler.h"
	83	#include "AliAnalysisManager.h"
	84
	85	#include "AliTracker.h"
	86	#include "AliMagF.h"
	87	#include "AliTPCCalROC.h"
	88
	89	#include "AliLog.h"
	90
	91	#include "AliTPCcalibTimeGain.h"
	92
	93	#include "TTreeStream.h"
	94	#include "AliTPCTracklet.h"
	95	#include "TTimeStamp.h"
	96	#include "AliTPCcalibDB.h"
	97	#include "AliTPCcalibLaser.h"
	98	#include "AliDCSSensorArray.h"
	99	#include "AliDCSSensor.h"
	100
	101	ClassImp(AliTPCcalibTimeGain)
	102
	103
	104	AliTPCcalibTimeGain::AliTPCcalibTimeGain()
	105	:AliTPCcalibBase(),
	106	fHistGainTime(0),
	107	fGainVsTime(0),
	108	fHistDeDxTotal(0),
	109	fIntegrationTimeDeDx(0),
	110	fMIP(0),
	111	fLowerTrunc(0),
	112	fUpperTrunc(0),
	113	fUseShapeNorm(0),
	114	fUsePosNorm(0),
	115	fUsePadNorm(0),
d0b31f57	116	fIsCosmic(0),
d0b31f57	117	fLowMemoryConsumption(0)
489dce1b	118	{
	119	AliInfo("Default Constructor");
	120	}
	121
	122
	123	AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t name, const Text_t title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain)
	124	:AliTPCcalibBase(),
	125	fHistGainTime(0),
	126	fGainVsTime(0),
	127	fHistDeDxTotal(0),
	128	fIntegrationTimeDeDx(0),
	129	fMIP(0),
	130	fLowerTrunc(0),
	131	fUpperTrunc(0),
	132	fUseShapeNorm(0),
	133	fUsePosNorm(0),
	134	fUsePadNorm(0),
d0b31f57	135	fIsCosmic(0),
d0b31f57	136	fLowMemoryConsumption(0)
489dce1b	137	{
	138
	139	SetName(name);
	140	SetTitle(title);
	141
	142	AliInfo("Non Default Constructor");
	143
	144	fIntegrationTimeDeDx = deltaIntegrationTimeGain;
	145
	146	Double_t deltaTime = EndTime - StartTime;
	147
	148
d0b31f57	149	// main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available)
	150	Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
	151	Int_t binsGainTime[5] = {100, timeBins, 2, 25, 200};
	152	Double_t xminGainTime[5] = {0.5, StartTime, 0.5, 0, 0.1};
489dce1b	153	Double_t xmaxGainTime[5] = { 4, EndTime, 2.5, 250, 50};
d0b31f57	154	fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx;dEdx,time,type,driftlength,momenta",5,binsGainTime,xminGainTime,xmaxGainTime);
489dce1b	155	BinLogX(fHistGainTime, 4);
	156	//
	157	fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,1000);
	158	BinLogX(fHistDeDxTotal);
	159
	160	// default values for dE/dx
	161	fMIP = 50.;
	162	fLowerTrunc = 0.0;
	163	fUpperTrunc = 0.7;
	164	fUseShapeNorm = kTRUE;
	165	fUsePosNorm = kFALSE;
	166	fUsePadNorm = kFALSE;
	167	//
	168	fIsCosmic = kTRUE;
d0b31f57	169	fLowMemoryConsumption = kTRUE;
489dce1b	170	//
	171
	172	}
	173
	174
	175
	176	AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
	177	//
	178	//
	179	//
	180	}
	181
	182
	183	void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
	184	//
	185	// main track loop
	186	//
	187	if (!event) {
	188	Printf("ERROR: ESD not available");
	189	return;
d0b31f57	190	}
	191
	192	if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
	193	ProcessCosmicEvent(event);
	194	} else {
	195	ProcessBeamEvent(event);
	196	}
	197
489dce1b	198
489dce1b	199
d0b31f57	200
	201	}
	202
	203
	204	void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
	205
489dce1b	206	AliESDfriend ESDfriend=static_cast<AliESDfriend>(event->FindListObject("AliESDfriend"));
	207	if (!ESDfriend) {
	208	Printf("ERROR: ESDfriend not available");
	209	return;
	210	}
	211	//
	212	UInt_t time = event->GetTimeStamp();
d0b31f57	213	Int_t ntracks = event->GetNumberOfTracks();
489dce1b	214	//
	215	// track loop
	216	//
	217	for (Int_t i=0;i<ntracks;++i) {
	218
	219	AliESDtrack *track = event->GetTrack(i);
	220	if (!track) continue;
	221
	222	const AliExternalTrackParam * trackIn = track->GetInnerParam();
	223	const AliExternalTrackParam * trackOut = track->GetOuterParam();
	224	if (!trackIn) continue;
	225	if (!trackOut) continue;
	226
	227	// calculate necessary track parameters
489dce1b	228	Double_t meanP = trackIn->GetP();
489dce1b	229	Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
489dce1b	230	Int_t NclsDeDx = track->GetTPCNcls();
489dce1b	231
489dce1b	232	// exclude tracks which do not look like primaries or are simply too short or on wrong sectors
d0b31f57	233	if (NclsDeDx < 60) continue;
489dce1b	234	if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
	235	if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
	236
	237	// Get seeds
	238	AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
	239	if (!friendTrack) continue;
	240	TObject *calibObject;
	241	AliTPCseed *seed = 0;
	242	for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
	243	if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
	244	}
	245
	246	if (seed) {
	247	Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
489dce1b	248	fHistDeDxTotal->Fill(meanP, TPCsignalMax);
489dce1b	249	//
d0b31f57	250	if (fLowMemoryConsumption) {
	251	if (meanP < 20) continue;
	252	meanP = 30; // set all momenta to one in order to save memory
489dce1b	253	}
d0b31f57	254	//dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
	255	Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP};
	256	fHistGainTime->Fill(vec);
	257	}
489dce1b	258
d0b31f57	259	}
489dce1b	260
d0b31f57	261	}
	262
	263
	264
	265	void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
	266
	267	AliESDfriend ESDfriend=static_cast<AliESDfriend>(event->FindListObject("AliESDfriend"));
	268	if (!ESDfriend) {
	269	Printf("ERROR: ESDfriend not available");
	270	return;
	271	}
	272	//
	273	UInt_t time = event->GetTimeStamp();
	274	Int_t ntracks = event->GetNumberOfTracks();
	275	//
	276	// track loop
	277	//
	278	for (Int_t i=0;i<ntracks;++i) {
	279
	280	AliESDtrack *track = event->GetTrack(i);
	281	if (!track) continue;
	282
	283	const AliExternalTrackParam * trackIn = track->GetInnerParam();
	284	const AliExternalTrackParam * trackOut = track->GetOuterParam();
	285	if (!trackIn) continue;
	286	if (!trackOut) continue;
	287
	288	// calculate necessary track parameters
	289	Double_t meanP = trackIn->GetP();
	290	Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
	291	Int_t NclsDeDx = track->GetTPCNcls();
	292
	293	// exclude tracks which do not look like primaries or are simply too short or on wrong sectors
	294	if (NclsDeDx < 60) continue;
	295	if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
	296	if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
	297
	298	// Get seeds
	299	AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
	300	if (!friendTrack) continue;
	301	TObject *calibObject;
	302	AliTPCseed *seed = 0;
	303	for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
	304	if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
	305	}
	306
	307	if (seed) {
	308	Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
	309	fHistDeDxTotal->Fill(meanP, TPCsignalMax);
	310	//
	311	if (fLowMemoryConsumption) {
	312	if (meanP > 0.5 \|\| meanP < 0.4) continue;
	313	meanP = 0.45; // set all momenta to one in order to save memory
	314	}
	315	//dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
	316	Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP};
	317	fHistGainTime->Fill(vec);
489dce1b	318	}
	319
	320	}
d0b31f57	321
489dce1b	322	}
	323
	324
	325	void AliTPCcalibTimeGain::Analyze() {
	326	//
	327	//
	328	//
489dce1b	329	if (fIsCosmic) {
d0b31f57	330	fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
d0b31f57	331	fHistGainTime->GetAxis(4)->SetRangeUser(20,100); // only Fermi-Plateau muons
489dce1b	332	} else {
d0b31f57	333	fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
d0b31f57	334	fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
489dce1b	335	}
489dce1b	336	//
d0b31f57	337	fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,2000,10);
489dce1b	338	//
	339	return;
	340	}
	341
	342
	343	Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
	344
	345	TIterator* iter = li->MakeIterator();
	346	AliTPCcalibTimeGain* cal = 0;
	347
	348	while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
	349	if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
	350	Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
	351	return -1;
	352	}
	353
	354	// add histograms here...
	355	if (cal->GetHistGainTime()) fHistGainTime->Add(cal->GetHistGainTime());
	356	if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
	357
	358	}
	359
	360	return 0;
	361
	362	}
	363
	364
489dce1b	365
	366	void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
	367
	368	// Method for the correct logarithmic binning of histograms
	369
	370	TAxis *axis = h->GetAxis(axisDim);
	371	int bins = axis->GetNbins();
	372
	373	Double_t from = axis->GetXmin();
	374	Double_t to = axis->GetXmax();
	375	Double_t *new_bins = new Double_t[bins + 1];
	376
	377	new_bins[0] = from;
	378	Double_t factor = pow(to/from, 1./bins);
	379
	380	for (int i = 1; i <= bins; i++) {
	381	new_bins[i] = factor * new_bins[i-1];
	382	}
	383	axis->Set(bins, new_bins);
	384	delete new_bins;
	385
	386	}
	387
	388
	389	void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
	390
	391	// Method for the correct logarithmic binning of histograms
	392
	393	TAxis *axis = h->GetXaxis();
	394	int bins = axis->GetNbins();
	395
	396	Double_t from = axis->GetXmin();
	397	Double_t to = axis->GetXmax();
	398	Double_t *new_bins = new Double_t[bins + 1];
	399
	400	new_bins[0] = from;
	401	Double_t factor = pow(to/from, 1./bins);
	402
	403	for (int i = 1; i <= bins; i++) {
	404	new_bins[i] = factor * new_bins[i-1];
	405	}
	406	axis->Set(bins, new_bins);
	407	delete new_bins;
	408
	409	}
	410
	411
	412
	413	void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F hist, Double_t ini) {
	414	//{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
	415	const Double_t sigma = 0.06;
	416
d0b31f57	417	TH2F histBG = new TH2F("histBG","dEdxBg; #beta #gamma; TPC signal (a.u.)",TMath::Nint(0.5hist->GetNbinsX()),0.1,5000.,hist->GetNbinsY(),hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
489dce1b	418	BinLogX(histBG);
489dce1b	419
d0b31f57	420	TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
489dce1b	421	TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
	422	TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
	423	TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
	424	TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
	425	foElectron->SetParameters(ini);
	426	foMuon->SetParameters(ini);
	427	foPion->SetParameters(ini);
	428	foKaon->SetParameters(ini);
	429	foProton->SetParameters(ini);
	430
	431	TCanvas *CanvCheck1 = new TCanvas();
	432	hist->Draw("colz");
	433	foElectron->Draw("same");
	434	foMuon->Draw("same");
	435	foPion->Draw("same");
	436	foKaon->Draw("same");
	437	foProton->Draw("same");
	438
	439	// Loop over all points of the input histogram
	440
	441	for(Int_t i=1; i < hist->GetNbinsX(); i++) {
	442	Double_t x = hist->GetXaxis()->GetBinCenter(i);
	443	for(Int_t j=1; j < hist->GetNbinsY(); j++) {
	444	Long64_t n = hist->GetBin(i, j);
	445	Double_t y = hist->GetYaxis()->GetBinCenter(j);
	446	Double_t entries = hist->GetBinContent(n);
	447	Double_t mass = 0;
	448
	449	// 1. identify protons
	450	mass = 0.938;
	451	if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
	452	for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	453	}
	454
	455	// 2. identify electrons
	456	mass = 0.000511;
	457	if (fIsCosmic) {
	458	if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
	459	for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	460	}
	461	} else {
d0b31f57	462	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))) {
489dce1b	463	for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	464	}
	465	}
	466
	467	// 3. identify either muons or pions depending on cosmic or not
	468	if (fIsCosmic) {
	469	mass = 0.1056;
	470	if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
	471	for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	472	}
	473	} else {
	474	mass = 0.1396;
	475	if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
	476	for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	477	}
	478	}
	479
	480	// 4. for pp also Kaons must be included
	481	if (!fIsCosmic) {
	482	mass = 0.4936;
	483	if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
	484	for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
	485	}
	486	}
	487	}
	488	}
	489
	490	// Fit Aleph-Parameters to the obtained profile
	491	TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
	492	funcAlephD->SetParameters(ini);
	493
	494	TCanvas *CanvCheck2 = new TCanvas();
	495	histBG->Draw();
	496
	497	//FitSlices
	498	TObjArray * arr = new TObjArray();
	499	histBG->FitSlicesY(0,0,-1,0,"QN",arr);
	500	TH1D * fitMean = (TH1D*) arr->At(1);
	501	fitMean->Draw("same");
	502
	503	funcAlephD->SetParLimits(2,1e-3,1e-7);
	504	funcAlephD->SetParLimits(3,0.5,3.5);
	505	funcAlephD->SetParLimits(4,0.5,3.5);
	506	fitMean->Fit(funcAlephD, "QNR");
	507	funcAlephD->Draw("same");
	508
	509	for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
	510
	511	foElectron->SetParameters(ini);
	512	foMuon->SetParameters(ini);
	513	foPion->SetParameters(ini);
	514	foKaon->SetParameters(ini);
	515	foProton->SetParameters(ini);
	516
	517	TCanvas *CanvCheck3 = new TCanvas();
	518	hist->Draw("colz");
	519	foElectron->Draw("same");
	520	foMuon->Draw("same");
	521	foPion->Draw("same");
	522	foKaon->Draw("same");
	523	foProton->Draw("same");
	524
	525	CanvCheck1->Print();
	526	CanvCheck2->Print();
527	CanvCheck3->Print();
528
529	return;
530
531
532	}