[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.                  *
 **************************************************************************/

/*

This class provides the calibration of the time dependence of the TPC gain due to pressure and temperature changes etc.


//0.  Libraries to load
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 = gain->GetGraphGainVsTime(0,500)

TH2D * GainTime = gain->GetHistGainTime()->Projection(0,1)
GainTime->GetXaxis()->SetTimeDisplay(kTRUE)
GainTime->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}")
GainTime->Draw("colz")

//
// MakeSlineFit example
//
AliSplineFit *fit = AliTPCcalibTimeGain::MakeSplineFit(gr)

TGraph * grfit = fit.MakeGraph(gr->GetX()[0],gr->GetX()[gr->GetN()-1],50000,0);

gr->SetMarkerStyle(25);
gr->Draw("lp");
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),
   fUseMax(0),
   fLowerTrunc(0),
   fUpperTrunc(0),
   fUseShapeNorm(0),
   fUsePosNorm(0),
   fUsePadNorm(0),
   fUseCookAnalytical(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),
   fUseMax(0),
   fLowerTrunc(0),
   fUpperTrunc(0),
   fUseShapeNorm(0),
   fUsePosNorm(0),
   fUsePadNorm(0),
   fUseCookAnalytical(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), run number
  Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
  Int_t binsGainTime[6]    = {100,  timeBins,    2,  25, 200, 10000000};
  Double_t xminGainTime[6] = {0.5, StartTime,  0.5,   0, 0.1,    -0.5};
  Double_t xmaxGainTime[6] = {  4,   EndTime,  2.5, 250,  50, 9999999.5};
  fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number;dEdx",6,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.;
  fUseMax = kTRUE;
  fLowerTrunc = 0.0;
  fUpperTrunc = 0.7;
  fUseShapeNorm = kTRUE;
  fUsePosNorm = kFALSE;
  fUsePadNorm = kFALSE;
  fUseCookAnalytical = kFALSE;
  //
  fIsCosmic = kTRUE;
  fLowMemoryConsumption = kFALSE;
  //
  
 }


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();
  Int_t runNumber = event->GetRunNumber();
  //
  // 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 TPCsignal = GetTPCdEdx(seed);
      fHistDeDxTotal->Fill(meanP, TPCsignal);
      //
      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[6] = {TPCsignal,time,1,meanDrift,meanP,runNumber};
      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();
  Int_t runNumber = event->GetRunNumber();
  //
  // 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 TPCsignal = GetTPCdEdx(seed);
      fHistDeDxTotal->Fill(meanP, TPCsignal);
      //
      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[6] = {TPCsignal,time,2,meanDrift,meanP,runNumber};
      fHistGainTime->Fill(vec);
    }
    
  }

}


Float_t AliTPCcalibTimeGain::GetTPCdEdx(AliTPCseed * seed) {

  Double_t signal = 0;
  //
  if (!fUseCookAnalytical) {
    signal = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,fUseMax,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
  } else {
    signal = (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax);
  }
  //
  return signal;
}


void AliTPCcalibTimeGain::AnalyzeRun(Int_t minEntries) {
  //
  //
  //
  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,minEntries,10);
  //
  return;
}


TGraphErrors * AliTPCcalibTimeGain::GetGraphGainVsTime(Int_t runNumber, Int_t minEntries) {
  //
  //
  //
  if (runNumber == 0) {
    if (!fGainVsTime) {
      AnalyzeRun(minEntries);
    }
  } else {
    // 1st check if the current run was cosmic or beam event
    fHistGainTime->GetAxis(5)->SetRangeUser(runNumber,runNumber);
    AnalyzeRun(minEntries);
  }
  if (fGainVsTime->GetN() == 0) return 0;
  return fGainVsTime;
}

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;
  
}


AliSplineFit * AliTPCcalibTimeGain::MakeSplineFit(TGraphErrors * graph) {
  //
  //
  //
  AliSplineFit *fit = new AliSplineFit();
  fit->SetGraph(graph);
  fit->SetMinPoints(graph->GetN()+1);
  fit->InitKnots(graph,2,0,0.001);
  fit->SplineFit(0);
  return fit;
  
}


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