[u/mrichter/AliRoot.git] / PWGPP / EvTrkSelection / AliAnalysisTrackingUncertainties.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, proviyaded 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 purapose. It is         *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
// Analysis task for the systematic study of the uncertainties related to   //
// the tracking and ITS-TPC matching efficiency for different particle      //
// species.                                                                 //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////


#include "Riostream.h"
#include "TH1F.h"
#include "TH2D.h"
#include "TH3F.h"
#include "THn.h"
#include "TList.h"
#include "TMath.h"
#include "TParticlePDG.h"
//
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisManager.h"
#include "AliPID.h"
#include "AliESDtrackCuts.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDtrack.h"
#include "AliESDpid.h"
#include "AliESDUtils.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliLog.h"
//
#include "AliAnalysisTrackingUncertainties.h"


ClassImp(AliAnalysisTrackingUncertainties)

// histogram constants
const Int_t kNumberOfAxes = 5;

//________________________________________________________________________
AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties() 
  : AliAnalysisTaskSE("TaskTestPA"), 
    fESD(0),
    fESDpid(0),
    fUtils(0),
    fMCtrue(0),
    fListHist(0),
    fESDtrackCuts(0),
    fMatchTr(),
    fMatchChi()
{
  // default Constructor
  /* fast compilation test

     gSystem->Load("libANALYSIS");
     gSystem->Load("libANALYSISalice");
     .L AliAnalysisTrackingUncertainties.cxx++
   */
}


//________________________________________________________________________
AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties(const char *name) 
  : AliAnalysisTaskSE(name),
    fESD(0),
    fESDpid(0),
    fUtils(0),
    fMCtrue(0),
    fListHist(0),
    fESDtrackCuts(0),
    fMatchTr(),
    fMatchChi()
{
  //
  // standard constructur which should be used
  //
  fMCtrue = kTRUE; 
  //
  // create track cuts
  //
  fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
  fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
  fESDtrackCuts->SetEtaRange(-1., 1.);
  //
  // analysis utils if needed
  //
  fUtils = new AliAnalysisUtils();
  //

  // Output slot #0 writes into a TList container
  DefineOutput(1, TList::Class());

}


//________________________________________________________________________
void AliAnalysisTrackingUncertainties::UserCreateOutputObjects() 
{
  //
  // Create histograms
  // Called once
  //
  fListHist = new TList();
  fListHist->SetOwner(kTRUE);
  //
  // basic QA and statistics histograms
  //
  TH2F * histVertexSelection = new TH2F("histVertexSelection", "vertex selection; vertex z (cm); accepted/rejected", 100, -50., 50., 2, -0.5, 1.5);
  fListHist->Add(histVertexSelection);
  //
  // track cut variation histograms
  InitializeTrackCutHistograms();
  //
  //
  // matching histograms
  //
  const int nbPt=40;
  const double ptMax=5;
  //
  TH2F * hNMatch    = new TH2F("hNMatch","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
  TH2F * hBestMatch = new TH2F("hBestMatch","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
  TH2F * hAllMatch  = new TH2F("hAllMatch","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
  TH2F * hAllMatchGlo  = new TH2F("hAllMatchGlo","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
  fListHist->Add(hNMatch);
  fListHist->Add(hBestMatch);
  fListHist->Add(hAllMatch);
  fListHist->Add(hAllMatchGlo);
  //
  TH2F * hNMatchBg    = new TH2F("hNMatchBg","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
  TH2F * hBestMatchBg = new TH2F("hBestMatchBg","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
  TH2F * hAllMatchBg  = new TH2F("hAllMatchBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
  TH2F * hAllMatchGloBg  = new TH2F("hAllMatchGloBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
  fListHist->Add(hNMatchBg);
  fListHist->Add(hBestMatchBg);
  fListHist->Add(hAllMatchBg);
  fListHist->Add(hAllMatchGloBg);
  //
  // post data
  //
  PostData(1, fListHist);


}


//________________________________________________________________________
void AliAnalysisTrackingUncertainties::UserExec(Option_t *) 
{
  //
  // main event loop
  //
  fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
  if (!fESD) {
    PostData(1, fListHist);
    return;
  }
  //
  if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
  //
  // Check Monte Carlo information and other access first:
  //
  AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  if (!eventHandler) {
    fMCtrue = kFALSE;
  }
  //
  // extract generated particles information
  //
  AliMCEvent* mcEvent = 0x0;
  AliStack* stack = 0x0;
  if (eventHandler) mcEvent = eventHandler->MCEvent();
  if (!mcEvent) {
    if (fMCtrue) return;
  }
  if (fMCtrue) {
    //
    stack = mcEvent->Stack();
    if (!stack) return;
    //
    for(Int_t i = 0; i < stack->GetNtrack(); i++) {
      /* at the moment nothing is needed here
      TParticle * trackMC = stack->Particle(i);
      Double_t rap = trackMC->Eta();
      Double_t y = trackMC->Y();
      Double_t pT  = trackMC->Pt();
      */
      //
    }

  }
  //
  if (!fESDtrackCuts) {
    PostData(1, fListHist);
    return;
  }
  //
  // monitor vertex position and selection
  //
  TH2F * histVertexSelection = (TH2F *) fListHist->FindObject("histVertexSelection");
  //
  Float_t vertexZ = 0.;
  if (IsVertexAccepted(fESD, vertexZ)) {
    histVertexSelection->Fill(vertexZ, 0);
  } else {
    histVertexSelection->Fill(vertexZ, 1);
    return;
  }
  //
  // fill track cut variation histograms
  //
  ProcessTrackCutVariation();
  //
  // fill ITS->TPC matching histograms
  //
  ProcessItsTpcMatching();
  
  // Post output data  
  PostData(1, fListHist);
  
}      


//________________________________________________________________________
void  AliAnalysisTrackingUncertainties::ProcessItsTpcMatching(){
  //
  // check how many its-sa tracks get matched to TPC
  //
  int ntr = fESD->GetNumberOfTracks();
  //
  // initialize histograms
  //
  TH2F * hNMatch      = (TH2F*) fListHist->FindObject("hNMatch");
  TH2F * hBestMatch   = (TH2F*) fListHist->FindObject("hBestMatch");
  TH2F * hAllMatch    = (TH2F*) fListHist->FindObject("hAllMatch");
  TH2F * hAllMatchGlo = (TH2F*) fListHist->FindObject("hAllMatchGlo");  
  //
  TH2F * hNMatchBg      = (TH2F*) fListHist->FindObject("hNMatchBg");
  TH2F * hBestMatchBg   = (TH2F*) fListHist->FindObject("hBestMatchBg");
  TH2F * hAllMatchBg    = (TH2F*) fListHist->FindObject("hAllMatchBg");
  TH2F * hAllMatchGloBg = (TH2F*) fListHist->FindObject("hAllMatchGloBg");    
  //
  for (int it=0;it<ntr;it++) {
    AliESDtrack* trSA = fESD->GetTrack(it);
    if (!trSA->IsOn(AliESDtrack::kITSpureSA) || !trSA->IsOn(AliESDtrack::kITSrefit)) continue;
    double pt = trSA->Pt();
    //
    Int_t nmatch = 0;
    for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
    for (int it1=0;it1<ntr;it1++) {
      if (it1==it) continue;
      AliESDtrack* trESD = fESD->GetTrack(it1);
      if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
      Match(trSA,trESD, nmatch);
    }
    //
    hNMatch->Fill(pt,nmatch);
    if (nmatch>0) hBestMatch->Fill(pt,fMatchChi[0]);
    for (int imt=nmatch;imt--;) {
      hAllMatch->Fill(pt,fMatchChi[imt]);
      if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGlo->Fill(pt,fMatchChi[imt]);
    }
    //
    nmatch = 0;
    for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
    for (int it1=0;it1<ntr;it1++) {
      if (it1==it) continue;
      AliESDtrack* trESD = fESD->GetTrack(it1);
      if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
      Match(trSA,trESD,TMath::Pi(), nmatch);
    }
    //
    hNMatchBg->Fill(pt,nmatch);
    if (nmatch>0) hBestMatchBg->Fill(pt,fMatchChi[0]);
    for (int imt=nmatch;imt--;) {
      hAllMatchBg->Fill(pt,fMatchChi[imt]);
      if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGloBg->Fill(pt,fMatchChi[imt]);
    }
    //
  }


}


void AliAnalysisTrackingUncertainties::Match(const AliESDtrack* tr0, const AliESDtrack* tr1,Double_t rotate, Int_t nmatch) {
  //
  // check if two tracks are matching, possible rotation for combinatoric backgr.
  // 
  Float_t bField = fESD->GetMagneticField();
  //
  const AliExternalTrackParam* trtpc0 = tr1->GetInnerParam();
  if (!trtpc0) return;
  AliExternalTrackParam trtpc(*trtpc0);
  //
  if (TMath::Abs(rotate)>1e-5) {
    const double *par = trtpc.GetParameter();
    const double *cov = trtpc.GetCovariance();
    double alp = trtpc.GetAlpha() + rotate;
    trtpc.Set(trtpc.GetX(),alp,par,cov);
  }
  //
  if (!trtpc.Rotate(tr0->GetAlpha())) return;
  if (!trtpc.PropagateTo(tr0->GetX(),bField)) return;
  double chi2 = tr0->GetPredictedChi2(&trtpc);
  if (chi2>kMaxChi2) return;
  //
  int ins;
  for (ins=0;ins<nmatch;ins++) if (chi2<fMatchChi[ins]) break;
  if (ins>=kMaxMatch) return;
  
  for (int imv=nmatch;imv>ins;imv--) {
    if (imv>=kMaxMatch) continue;
    fMatchTr[imv]  = fMatchTr[imv-1];
    fMatchChi[imv] = fMatchChi[imv-1];
  }
  fMatchTr[ins] = tr1;
  fMatchChi[ins] = chi2;
  nmatch++;
  if (nmatch>=kMaxMatch) nmatch = kMaxMatch;
  //
}


//________________________________________________________________________
void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
  //
  // fill track cut variation histograms - undo cuts step-by-step and fill histograms
  //
  //
  // initialize histograms
  //
  THnF * histNcl     = (THnF *) fListHist->FindObject("histNcl");
  THnF * histChi2Tpc = (THnF *) fListHist->FindObject("histChi2Tpc");
  THnF * histDcaZ    = (THnF *) fListHist->FindObject("histDcaZ");
  THnF * histSpd     = (THnF *) fListHist->FindObject("histSpd");
  //
  Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
  //
  for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
    //
    AliESDtrack *track =fESD->GetTrack(i);
    //
    // relevant variables
    //
    Int_t nclsTPC = track->GetTPCncls();
    Float_t pT  = track->Pt();
    Float_t eta = track->Eta();
    Float_t phi = track->Phi();
    Float_t chi2TPC = track->GetTPCchi2();
    if (nclsTPC != 0) {
      chi2TPC /= nclsTPC; 
    } else {
      chi2TPC = 999.;
    }
    //
    track->GetImpactParameters(dca, cov);
    //
    // (1.) fill number of clusters histogram
    //
    Int_t minNclsTPC = fESDtrackCuts->GetMinNClusterTPC();
    fESDtrackCuts->SetMinNClustersTPC(0);
    if (fESDtrackCuts->AcceptTrack(track)) {
      Double_t vecHistNcl[kNumberOfAxes] = {nclsTPC, pT, eta, phi, 0.};
      histNcl->Fill(vecHistNcl);
    }
    fESDtrackCuts->SetMinNClustersTPC(minNclsTPC);
    //
    // (2.) fill chi2 TPC histogram
    //
    Float_t maxChi2 = fESDtrackCuts->GetMaxChi2PerClusterTPC();
    fESDtrackCuts->SetMaxChi2PerClusterTPC(999.);
    if (fESDtrackCuts->AcceptTrack(track)) {
      Double_t vecHistChi2Tpc[kNumberOfAxes] = {chi2TPC, pT, eta, phi, 0.};
      histChi2Tpc->Fill(vecHistChi2Tpc);
    }
    fESDtrackCuts->SetMaxChi2PerClusterTPC(maxChi2);
    //
    // (3.) fill dca_z histogram
    //
    Float_t maxDcaZ = fESDtrackCuts->GetMaxDCAToVertexZ();
    fESDtrackCuts->SetMaxDCAToVertexZ(999.);
    if (fESDtrackCuts->AcceptTrack(track)) {
      Double_t vecHistDcaZ[kNumberOfAxes] = {dca[1], pT, eta, phi, 0.};
      histDcaZ->Fill(vecHistDcaZ);
    }
    fESDtrackCuts->SetMaxDCAToVertexZ(maxDcaZ);
    //
    // (4.) fill hit in SPD histogram
    //
    fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
    if (fESDtrackCuts->AcceptTrack(track)) {
      Int_t hasPoint = 0;
      if (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) hasPoint = 1;
      Double_t vecHistSpd[kNumberOfAxes] = {hasPoint, pT, eta, phi, 0.};
      histSpd->Fill(vecHistSpd);
    }
    fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);


  } // end of track loop


}


//________________________________________________________________________
void AliAnalysisTrackingUncertainties::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query


}


//________________________________________________________________________
void AliAnalysisTrackingUncertainties::BinLogAxis(const THn *h, Int_t axisNumber) {
  //
  // Method for the correct logarithmic binning of histograms
  //
  TAxis *axis = h->GetAxis(axisNumber);
  int bins = axis->GetNbins();

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


//________________________________________________________________________
Bool_t AliAnalysisTrackingUncertainties::IsVertexAccepted(AliESDEvent * esd, Float_t &vertexZ) {
  //
  // function to check if a proper vertex is reconstructed and write z-position in vertexZ
  //
  vertexZ = -999.;
  Bool_t vertexOkay = kFALSE;
  const AliESDVertex *vertex = esd->GetPrimaryVertexTracks();
  if (vertex->GetNContributors() < 1) {
    //
    vertex = esd->GetPrimaryVertexSPD();
    if (vertex->GetNContributors() < 1) {
      vertexOkay = kFALSE; }
    else {
      vertexOkay = kTRUE;
    }
    //
    TString vtxTyp = vertex->GetTitle();
    Double_t cov[6]={0};
    vertex->GetCovarianceMatrix(cov);
    Double_t zRes = TMath::Sqrt(cov[5]);
    if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) vertexOkay = kFALSE;
  }
  else {
    vertexOkay = kTRUE;
  }

  vertexZ = vertex->GetZ();  
  return vertexOkay;

}


//________________________________________________________________________
void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
  //
  // create histograms for the track cut studies
  //
  //
  // (1.) number of clusters       
  //                               0-ncl, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
  Int_t    binsNcl[kNumberOfAxes] = { 40,   50,    20,            18,  6};
  Double_t minNcl[kNumberOfAxes]  = { 20,  0.1,    -1,             0, -0.5};
  Double_t maxNcl[kNumberOfAxes]  = {160,   20,    +1, 2*TMath::Pi(),  5.5};
  //
  TString axisNameNcl[kNumberOfAxes]  = {"ncl","pT","eta","phi","pid"};
  TString axisTitleNcl[kNumberOfAxes] = {"ncl","pT","eta","phi","pid"};
  //
  THnF * histNcl = new THnF("histNcl","number of clusters histogram",kNumberOfAxes, binsNcl, minNcl, maxNcl);
  BinLogAxis(histNcl, 1);
  fListHist->Add(histNcl);
  //
  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
    histNcl->GetAxis(iaxis)->SetName(axisNameNcl[iaxis]);
    histNcl->GetAxis(iaxis)->SetTitle(axisTitleNcl[iaxis]);
  }
  //
  // (2.) chi2/cls-TPC            
  //                                  0-chi2, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
  Int_t    binsChi2Tpc[kNumberOfAxes] = { 40,   50,    20,            18,  6};
  Double_t minChi2Tpc[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
  Double_t maxChi2Tpc[kNumberOfAxes]  = {  8,   20,    +1, 2*TMath::Pi(),  5.5};
  //
  TString axisNameChi2Tpc[kNumberOfAxes]  = {"chi2tpc","pT","eta","phi","pid"};
  TString axisTitleChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
  //
  THnF * histChi2Tpc = new THnF("histChi2Tpc","number of clusters histogram",kNumberOfAxes, binsChi2Tpc, minChi2Tpc, maxChi2Tpc);
  BinLogAxis(histChi2Tpc, 1);
  fListHist->Add(histChi2Tpc);
  //
  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
    histChi2Tpc->GetAxis(iaxis)->SetName(axisNameChi2Tpc[iaxis]);
    histChi2Tpc->GetAxis(iaxis)->SetTitle(axisTitleChi2Tpc[iaxis]);
  }
  //
  // (3.) dca_z
  //                               0-dcaZ, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
  Int_t    binsDcaZ[kNumberOfAxes] = { 10,   50,    20,            18,  6};
  Double_t minDcaZ[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
  Double_t maxDcaZ[kNumberOfAxes]  = {  5,   20,    +1, 2*TMath::Pi(),  5.5};
  //
  TString axisNameDcaZ[kNumberOfAxes]  = {"dcaZ","pT","eta","phi","pid"};
  TString axisTitleDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
  //
  THnF * histDcaZ = new THnF("histDcaZ","number of clusters histogram",kNumberOfAxes, binsDcaZ, minDcaZ, maxDcaZ);
  BinLogAxis(histDcaZ, 1);
  fListHist->Add(histDcaZ);
  //
  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
    histDcaZ->GetAxis(iaxis)->SetName(axisNameDcaZ[iaxis]);
    histDcaZ->GetAxis(iaxis)->SetTitle(axisTitleDcaZ[iaxis]);
  }
  //
  // (4.) hit in SPD layer
  //                              0-spdHit, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
  Int_t    binsSpd[kNumberOfAxes] = {    2,   50,    20,            18,  6};
  Double_t minSpd[kNumberOfAxes]  = { -0.5,  0.1,    -1,             0, -0.5};
  Double_t maxSpd[kNumberOfAxes]  = {  1.5,   20,    +1, 2*TMath::Pi(),  5.5};
  //
  TString axisNameSpd[kNumberOfAxes]  = {"spdHit","pT","eta","phi","pid"};
  TString axisTitleSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
  //
  THnF * histSpd = new THnF("histSpd","number of clusters histogram",kNumberOfAxes, binsSpd, minSpd, maxSpd);
  BinLogAxis(histSpd, 1);
  fListHist->Add(histSpd);
  //
  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
    histSpd->GetAxis(iaxis)->SetName(axisNameSpd[iaxis]);
    histSpd->GetAxis(iaxis)->SetTitle(axisTitleSpd[iaxis]);
  }


}
Commit	Line	Data
b938b7fc	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, proviyaded 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 purapose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//////////////////////////////////////////////////////////////////////////////
	17	// //
	18	// Analysis task for the systematic study of the uncertainties related to //
	19	// the tracking and ITS-TPC matching efficiency for different particle //
	20	// species. //
	21	// //
	22	//////////////////////////////////////////////////////////////////////////////
	23
	24
	25	#include "Riostream.h"
	26	#include "TH1F.h"
	27	#include "TH2D.h"
	28	#include "TH3F.h"
	29	#include "THn.h"
	30	#include "TList.h"
	31	#include "TMath.h"
	32	#include "TParticlePDG.h"
	33	//
	34	#include "AliAnalysisTaskSE.h"
	35	#include "AliAnalysisManager.h"
	36	#include "AliPID.h"
	37	#include "AliESDtrackCuts.h"
	38	#include "AliESDVertex.h"
	39	#include "AliESDEvent.h"
	40	#include "AliESDInputHandler.h"
	41	#include "AliESDtrack.h"
	42	#include "AliESDpid.h"
	43	#include "AliESDUtils.h"
	44	#include "AliMCEventHandler.h"
	45	#include "AliMCEvent.h"
	46	#include "AliStack.h"
	47	#include "AliLog.h"
	48	//
	49	#include "AliAnalysisTrackingUncertainties.h"
	50
	51
	52	ClassImp(AliAnalysisTrackingUncertainties)
	53
	54	// histogram constants
	55	const Int_t kNumberOfAxes = 5;
	56
	57	//________________________________________________________________________
	58	AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties()
	59	: AliAnalysisTaskSE("TaskTestPA"),
	60	fESD(0),
	61	fESDpid(0),
	62	fUtils(0),
	63	fMCtrue(0),
	64	fListHist(0),
65	fESDtrackCuts(0),
66	fMatchTr(),
67	fMatchChi()
68	{
69	// default Constructor
70	/* fast compilation test
71
72	gSystem->Load("libANALYSIS");
73	gSystem->Load("libANALYSISalice");
74	.L AliAnalysisTrackingUncertainties.cxx++
75	*/
76	}
77
78
79	//________________________________________________________________________
80	AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties(const char *name)
81	: AliAnalysisTaskSE(name),
82	fESD(0),
83	fESDpid(0),
84	fUtils(0),
85	fMCtrue(0),
86	fListHist(0),
87	fESDtrackCuts(0),
88	fMatchTr(),
89	fMatchChi()
90	{
91	//
92	// standard constructur which should be used
93	//
94	fMCtrue = kTRUE;
95	//
96	// create track cuts
97	//
98	fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
99	fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
100	fESDtrackCuts->SetEtaRange(-1., 1.);
101	//
102	// analysis utils if needed
103	//
104	fUtils = new AliAnalysisUtils();
105	//
106
107	// Output slot #0 writes into a TList container
108	DefineOutput(1, TList::Class());
109
110	}
111
112
113	//________________________________________________________________________
114	void AliAnalysisTrackingUncertainties::UserCreateOutputObjects()
115	{
116	//
117	// Create histograms
118	// Called once
119	//
120	fListHist = new TList();
121	fListHist->SetOwner(kTRUE);
122	//
123	// basic QA and statistics histograms
124	//
125	TH2F * histVertexSelection = new TH2F("histVertexSelection", "vertex selection; vertex z (cm); accepted/rejected", 100, -50., 50., 2, -0.5, 1.5);
126	fListHist->Add(histVertexSelection);
127	//
128	// track cut variation histograms
129	InitializeTrackCutHistograms();
130	//
131	//
132	// matching histograms
133	//
134	const int nbPt=40;
135	const double ptMax=5;
136	//
137	TH2F * hNMatch = new TH2F("hNMatch","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
138	TH2F * hBestMatch = new TH2F("hBestMatch","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
139	TH2F * hAllMatch = new TH2F("hAllMatch","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
140	TH2F * hAllMatchGlo = new TH2F("hAllMatchGlo","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
141	fListHist->Add(hNMatch);
142	fListHist->Add(hBestMatch);
143	fListHist->Add(hAllMatch);
144	fListHist->Add(hAllMatchGlo);
145	//
146	TH2F * hNMatchBg = new TH2F("hNMatchBg","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
147	TH2F * hBestMatchBg = new TH2F("hBestMatchBg","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
148	TH2F * hAllMatchBg = new TH2F("hAllMatchBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
149	TH2F * hAllMatchGloBg = new TH2F("hAllMatchGloBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
150	fListHist->Add(hNMatchBg);
151	fListHist->Add(hBestMatchBg);
152	fListHist->Add(hAllMatchBg);
153	fListHist->Add(hAllMatchGloBg);
154	//
155	// post data
156	//
157	PostData(1, fListHist);
158
159
160	}
161
162
163
164	//________________________________________________________________________
165	void AliAnalysisTrackingUncertainties::UserExec(Option_t *)
166	{
167	//
168	// main event loop
169	//
170	fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
171	if (!fESD) {
172	PostData(1, fListHist);
173	return;
174	}
175	//
176	if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
177	//
178	// Check Monte Carlo information and other access first:
179	//
180	AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
181	if (!eventHandler) {
182	fMCtrue = kFALSE;
183	}
184	//
185	// extract generated particles information
186	//
187	AliMCEvent* mcEvent = 0x0;
188	AliStack* stack = 0x0;
189	if (eventHandler) mcEvent = eventHandler->MCEvent();
190	if (!mcEvent) {
191	if (fMCtrue) return;
192	}
193	if (fMCtrue) {
194	//
195	stack = mcEvent->Stack();
196	if (!stack) return;
197	//
198	for(Int_t i = 0; i < stack->GetNtrack(); i++) {
199	/* at the moment nothing is needed here
200	TParticle * trackMC = stack->Particle(i);
201	Double_t rap = trackMC->Eta();
202	Double_t y = trackMC->Y();
203	Double_t pT = trackMC->Pt();
204	*/
205	//
206	}
207
208	}
209	//
210	if (!fESDtrackCuts) {
211	PostData(1, fListHist);
212	return;
213	}
214	//
215	// monitor vertex position and selection
216	//
217	TH2F * histVertexSelection = (TH2F *) fListHist->FindObject("histVertexSelection");
218	//
219	Float_t vertexZ = 0.;
220	if (IsVertexAccepted(fESD, vertexZ)) {
221	histVertexSelection->Fill(vertexZ, 0);
222	} else {
223	histVertexSelection->Fill(vertexZ, 1);
224	return;
225	}
226	//
227	// fill track cut variation histograms
228	//
229	ProcessTrackCutVariation();
230	//
231	// fill ITS->TPC matching histograms
232	//
233	ProcessItsTpcMatching();
234
235	// Post output data
236	PostData(1, fListHist);
237
238	}
239
240
241	//________________________________________________________________________
242	void AliAnalysisTrackingUncertainties::ProcessItsTpcMatching(){
243	//
244	// check how many its-sa tracks get matched to TPC
245	//
246	int ntr = fESD->GetNumberOfTracks();
247	//
248	// initialize histograms
249	//
250	TH2F * hNMatch = (TH2F*) fListHist->FindObject("hNMatch");
251	TH2F * hBestMatch = (TH2F*) fListHist->FindObject("hBestMatch");
252	TH2F * hAllMatch = (TH2F*) fListHist->FindObject("hAllMatch");
253	TH2F * hAllMatchGlo = (TH2F*) fListHist->FindObject("hAllMatchGlo");
254	//
255	TH2F * hNMatchBg = (TH2F*) fListHist->FindObject("hNMatchBg");
256	TH2F * hBestMatchBg = (TH2F*) fListHist->FindObject("hBestMatchBg");
257	TH2F * hAllMatchBg = (TH2F*) fListHist->FindObject("hAllMatchBg");
258	TH2F * hAllMatchGloBg = (TH2F*) fListHist->FindObject("hAllMatchGloBg");
259	//
260	for (int it=0;it<ntr;it++) {
261	AliESDtrack* trSA = fESD->GetTrack(it);
262	if (!trSA->IsOn(AliESDtrack::kITSpureSA) \|\| !trSA->IsOn(AliESDtrack::kITSrefit)) continue;
263	double pt = trSA->Pt();
264	//
265	Int_t nmatch = 0;
266	for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
267	for (int it1=0;it1<ntr;it1++) {
268	if (it1==it) continue;
269	AliESDtrack* trESD = fESD->GetTrack(it1);
270	if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
271	Match(trSA,trESD, nmatch);
272	}
273	//
274	hNMatch->Fill(pt,nmatch);
275	if (nmatch>0) hBestMatch->Fill(pt,fMatchChi[0]);
276	for (int imt=nmatch;imt--;) {
277	hAllMatch->Fill(pt,fMatchChi[imt]);
278	if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGlo->Fill(pt,fMatchChi[imt]);
279	}
280	//
281	nmatch = 0;
282	for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
283	for (int it1=0;it1<ntr;it1++) {
284	if (it1==it) continue;
285	AliESDtrack* trESD = fESD->GetTrack(it1);
286	if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
287	Match(trSA,trESD,TMath::Pi(), nmatch);
288	}
289	//
290	hNMatchBg->Fill(pt,nmatch);
291	if (nmatch>0) hBestMatchBg->Fill(pt,fMatchChi[0]);
292	for (int imt=nmatch;imt--;) {
293	hAllMatchBg->Fill(pt,fMatchChi[imt]);
294	if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGloBg->Fill(pt,fMatchChi[imt]);
295	}
296	//
297	}
298
299
300	}
301
302
303	void AliAnalysisTrackingUncertainties::Match(const AliESDtrack* tr0, const AliESDtrack* tr1,Double_t rotate, Int_t nmatch) {
304	//
305	// check if two tracks are matching, possible rotation for combinatoric backgr.
306	//
307	Float_t bField = fESD->GetMagneticField();
308	//
309	const AliExternalTrackParam* trtpc0 = tr1->GetInnerParam();
310	if (!trtpc0) return;
311	AliExternalTrackParam trtpc(*trtpc0);
312	//
313	if (TMath::Abs(rotate)>1e-5) {
314	const double *par = trtpc.GetParameter();
315	const double *cov = trtpc.GetCovariance();
316	double alp = trtpc.GetAlpha() + rotate;
317	trtpc.Set(trtpc.GetX(),alp,par,cov);
318	}
319	//
320	if (!trtpc.Rotate(tr0->GetAlpha())) return;
321	if (!trtpc.PropagateTo(tr0->GetX(),bField)) return;
322	double chi2 = tr0->GetPredictedChi2(&trtpc);
323	if (chi2>kMaxChi2) return;
324	//
325	int ins;
326	for (ins=0;ins<nmatch;ins++) if (chi2<fMatchChi[ins]) break;
327	if (ins>=kMaxMatch) return;
328
329	for (int imv=nmatch;imv>ins;imv--) {
330	if (imv>=kMaxMatch) continue;
331	fMatchTr[imv] = fMatchTr[imv-1];
332	fMatchChi[imv] = fMatchChi[imv-1];
333	}
334	fMatchTr[ins] = tr1;
335	fMatchChi[ins] = chi2;
336	nmatch++;
337	if (nmatch>=kMaxMatch) nmatch = kMaxMatch;
338	//
339	}
340
341
342	//________________________________________________________________________
343	void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
344	//
345	// fill track cut variation histograms - undo cuts step-by-step and fill histograms
346	//
347	//
348	// initialize histograms
349	//
350	THnF * histNcl = (THnF *) fListHist->FindObject("histNcl");
351	THnF * histChi2Tpc = (THnF *) fListHist->FindObject("histChi2Tpc");
352	THnF * histDcaZ = (THnF *) fListHist->FindObject("histDcaZ");
353	THnF * histSpd = (THnF *) fListHist->FindObject("histSpd");
354	//
355	Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
356	//
357	for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
358	//
359	AliESDtrack *track =fESD->GetTrack(i);
360	//
361	// relevant variables
362	//
363	Int_t nclsTPC = track->GetTPCncls();
364	Float_t pT = track->Pt();
365	Float_t eta = track->Eta();
366	Float_t phi = track->Phi();
367	Float_t chi2TPC = track->GetTPCchi2();
368	if (nclsTPC != 0) {
369	chi2TPC /= nclsTPC;
370	} else {
371	chi2TPC = 999.;
372	}
373	//
374	track->GetImpactParameters(dca, cov);
375	//
376	// (1.) fill number of clusters histogram
377	//
378	Int_t minNclsTPC = fESDtrackCuts->GetMinNClusterTPC();
379	fESDtrackCuts->SetMinNClustersTPC(0);
380	if (fESDtrackCuts->AcceptTrack(track)) {
381	Double_t vecHistNcl[kNumberOfAxes] = {nclsTPC, pT, eta, phi, 0.};
382	histNcl->Fill(vecHistNcl);
383	}
384	fESDtrackCuts->SetMinNClustersTPC(minNclsTPC);
385	//
386	// (2.) fill chi2 TPC histogram
387	//
388	Float_t maxChi2 = fESDtrackCuts->GetMaxChi2PerClusterTPC();
389	fESDtrackCuts->SetMaxChi2PerClusterTPC(999.);
390	if (fESDtrackCuts->AcceptTrack(track)) {
391	Double_t vecHistChi2Tpc[kNumberOfAxes] = {chi2TPC, pT, eta, phi, 0.};
392	histChi2Tpc->Fill(vecHistChi2Tpc);
393	}
394	fESDtrackCuts->SetMaxChi2PerClusterTPC(maxChi2);
395	//
396	// (3.) fill dca_z histogram
397	//
398	Float_t maxDcaZ = fESDtrackCuts->GetMaxDCAToVertexZ();
399	fESDtrackCuts->SetMaxDCAToVertexZ(999.);
400	if (fESDtrackCuts->AcceptTrack(track)) {
401	Double_t vecHistDcaZ[kNumberOfAxes] = {dca[1], pT, eta, phi, 0.};
402	histDcaZ->Fill(vecHistDcaZ);
403	}
404	fESDtrackCuts->SetMaxDCAToVertexZ(maxDcaZ);
405	//
406	// (4.) fill hit in SPD histogram
407	//
408	fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
409	if (fESDtrackCuts->AcceptTrack(track)) {
410	Int_t hasPoint = 0;
411	if (track->HasPointOnITSLayer(0) \|\| track->HasPointOnITSLayer(1)) hasPoint = 1;
412	Double_t vecHistSpd[kNumberOfAxes] = {hasPoint, pT, eta, phi, 0.};
413	histSpd->Fill(vecHistSpd);
414	}
415	fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
416
417
418	} // end of track loop
419
420
421	}
422
423
424
425	//________________________________________________________________________
426	void AliAnalysisTrackingUncertainties::Terminate(Option_t *)
427	{
428	// Draw result to the screen
429	// Called once at the end of the query
430
431
432	}
433
434
435	//________________________________________________________________________
436	void AliAnalysisTrackingUncertainties::BinLogAxis(const THn *h, Int_t axisNumber) {
437	//
438	// Method for the correct logarithmic binning of histograms
439	//
440	TAxis *axis = h->GetAxis(axisNumber);
441	int bins = axis->GetNbins();
442
443	Double_t from = axis->GetXmin();
444	Double_t to = axis->GetXmax();
445	Double_t *newBins = new Double_t[bins + 1];
446
447	newBins[0] = from;
448	Double_t factor = pow(to/from, 1./bins);
449
450	for (int i = 1; i <= bins; i++) {
451	newBins[i] = factor * newBins[i-1];
452	}
453	axis->Set(bins, newBins);
454	delete [] newBins;
455
456	}
457
458
459	//________________________________________________________________________
460	Bool_t AliAnalysisTrackingUncertainties::IsVertexAccepted(AliESDEvent * esd, Float_t &vertexZ) {
461	//
462	// function to check if a proper vertex is reconstructed and write z-position in vertexZ
463	//
464	vertexZ = -999.;
465	Bool_t vertexOkay = kFALSE;
466	const AliESDVertex *vertex = esd->GetPrimaryVertexTracks();
467	if (vertex->GetNContributors() < 1) {
468	//
469	vertex = esd->GetPrimaryVertexSPD();
470	if (vertex->GetNContributors() < 1) {
471	vertexOkay = kFALSE; }
472	else {
473	vertexOkay = kTRUE;
474	}
475	//
476	TString vtxTyp = vertex->GetTitle();
477	Double_t cov[6]={0};
478	vertex->GetCovarianceMatrix(cov);
479	Double_t zRes = TMath::Sqrt(cov[5]);
480	if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) vertexOkay = kFALSE;
481	}
482	else {
483	vertexOkay = kTRUE;
484	}
485
486	vertexZ = vertex->GetZ();
487	return vertexOkay;
488
489	}
490
491
492	//________________________________________________________________________
493	void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
494	//
495	// create histograms for the track cut studies
496	//
497	//
498	// (1.) number of clusters
499	// 0-ncl, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
500	Int_t binsNcl[kNumberOfAxes] = { 40, 50, 20, 18, 6};
501	Double_t minNcl[kNumberOfAxes] = { 20, 0.1, -1, 0, -0.5};
502	Double_t maxNcl[kNumberOfAxes] = {160, 20, +1, 2*TMath::Pi(), 5.5};
503	//
504	TString axisNameNcl[kNumberOfAxes] = {"ncl","pT","eta","phi","pid"};
505	TString axisTitleNcl[kNumberOfAxes] = {"ncl","pT","eta","phi","pid"};
506	//
507	THnF * histNcl = new THnF("histNcl","number of clusters histogram",kNumberOfAxes, binsNcl, minNcl, maxNcl);
508	BinLogAxis(histNcl, 1);
509	fListHist->Add(histNcl);
510	//
511	for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
512	histNcl->GetAxis(iaxis)->SetName(axisNameNcl[iaxis]);
513	histNcl->GetAxis(iaxis)->SetTitle(axisTitleNcl[iaxis]);
514	}
515	//
516	// (2.) chi2/cls-TPC
517	// 0-chi2, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
518	Int_t binsChi2Tpc[kNumberOfAxes] = { 40, 50, 20, 18, 6};
519	Double_t minChi2Tpc[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
520	Double_t maxChi2Tpc[kNumberOfAxes] = { 8, 20, +1, 2*TMath::Pi(), 5.5};
521	//
522	TString axisNameChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
523	TString axisTitleChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
524	//
525	THnF * histChi2Tpc = new THnF("histChi2Tpc","number of clusters histogram",kNumberOfAxes, binsChi2Tpc, minChi2Tpc, maxChi2Tpc);
526	BinLogAxis(histChi2Tpc, 1);
527	fListHist->Add(histChi2Tpc);
528	//
529	for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
530	histChi2Tpc->GetAxis(iaxis)->SetName(axisNameChi2Tpc[iaxis]);
531	histChi2Tpc->GetAxis(iaxis)->SetTitle(axisTitleChi2Tpc[iaxis]);
532	}
533	//
534	// (3.) dca_z
535	// 0-dcaZ, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
536	Int_t binsDcaZ[kNumberOfAxes] = { 10, 50, 20, 18, 6};
537	Double_t minDcaZ[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
538	Double_t maxDcaZ[kNumberOfAxes] = { 5, 20, +1, 2*TMath::Pi(), 5.5};
539	//
540	TString axisNameDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
541	TString axisTitleDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
542	//
543	THnF * histDcaZ = new THnF("histDcaZ","number of clusters histogram",kNumberOfAxes, binsDcaZ, minDcaZ, maxDcaZ);
544	BinLogAxis(histDcaZ, 1);
545	fListHist->Add(histDcaZ);
546	//
547	for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
548	histDcaZ->GetAxis(iaxis)->SetName(axisNameDcaZ[iaxis]);
549	histDcaZ->GetAxis(iaxis)->SetTitle(axisTitleDcaZ[iaxis]);
550	}
551	//
552	// (4.) hit in SPD layer
553	// 0-spdHit, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
554	Int_t binsSpd[kNumberOfAxes] = { 2, 50, 20, 18, 6};
555	Double_t minSpd[kNumberOfAxes] = { -0.5, 0.1, -1, 0, -0.5};
556	Double_t maxSpd[kNumberOfAxes] = { 1.5, 20, +1, 2*TMath::Pi(), 5.5};
557	//
558	TString axisNameSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
559	TString axisTitleSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
560	//
561	THnF * histSpd = new THnF("histSpd","number of clusters histogram",kNumberOfAxes, binsSpd, minSpd, maxSpd);
562	BinLogAxis(histSpd, 1);
563	fListHist->Add(histSpd);
564	//
565	for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
566	histSpd->GetAxis(iaxis)->SetName(axisNameSpd[iaxis]);
567	histSpd->GetAxis(iaxis)->SetTitle(axisTitleSpd[iaxis]);
568	}
569
570
571
572
573	}
574