[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckESD.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.                  *
**************************************************************************/

/////////////////////////////////////////////////////
//
// Check basic detector results at ESD level
//   - Geometrical efficiency  
//   - Tracking efficiency  
//   - PID efficiency  
//   - Refit efficiency  
//
// Author
//   Alex Bercuci <A.Bercuci@gsi.de>
//
//////////////////////////////////////////////////////

#include <TClonesArray.h>
#include <TObjArray.h>
#include <TObject.h>
#include <TH2I.h>
#include <TGraphErrors.h>
#include <TGraphAsymmErrors.h>
#include <TFile.h>
#include <TTree.h>
#include <TROOT.h>
#include <TChain.h>
#include <TParticle.h>

#include "AliLog.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliESDkink.h"
#include "AliMCEvent.h"
#include "AliESDInputHandler.h"
#include "AliMCEventHandler.h"

#include "AliESDtrack.h"
#include "AliMCParticle.h"
#include "AliPID.h"
#include "AliStack.h"
#include "AliTrackReference.h"

#include "AliTRDcheckESD.h"

ClassImp(AliTRDcheckESD)

const Float_t AliTRDcheckESD::fgkxTPC = 290.;
const Float_t AliTRDcheckESD::fgkxTOF = 365.;
FILE* AliTRDcheckESD::fgFile = 0x0;

//____________________________________________________________________
AliTRDcheckESD::AliTRDcheckESD():
  AliAnalysisTask("checkESD", "ESD checker for TRD info")
  ,fStatus(0)
  ,fESD(0x0)
  ,fMC(0x0)
  ,fHistos(0x0)
  ,fResults(0x0)
{
  //
  // Default constructor
  //
  SetMC(kTRUE);
  DefineInput(0, TChain::Class());
  DefineOutput(0, TObjArray::Class());
}

//____________________________________________________________________
AliTRDcheckESD::~AliTRDcheckESD()
{
// Destructor
  if(fHistos){
    //fHistos->Delete();
    delete fHistos;
  }
  if(fResults){
    fResults->Delete();
    delete fResults;
  }
}

//____________________________________________________________________
void AliTRDcheckESD::ConnectInputData(Option_t *)
{
  //
  // Link the Input Data
  //
  TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
  if(tree) tree->SetBranchStatus("Tracks", 1);

  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  fESD = esdH ? esdH->GetEvent() : 0x0;

  if(!HasMC()) return;
  AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  fMC = mcH ? mcH->MCEvent() : 0x0;
}

//____________________________________________________________________
void AliTRDcheckESD::CreateOutputObjects()
{	
  //
  // Create Output Containers (TObjectArray containing 1D histograms)
  //
  OpenFile(0, "RECREATE");  
  Histos();
}

//____________________________________________________________________
TGraph* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
{
// Retrieve graph with "id"
// Possible options are :
//   "b" - build graph if none found
//   "c" - clear existing graph

  Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
         kCLEAR = strstr(opt, "c"); // clear existing graph

  const Char_t *name[] = {
    "Geo", "Trk", "Pid", "Ref", "Max06", "Mean09"
  };
  const Char_t *title[] = {
    "TRD geometrical efficiency (TRDin/TPCout)"
    ,"TRD tracking efficiency (TRDout/TRDin)"
    ,"TRD PID efficiency (TRDpid/TRDin)"
    ,"TRD refit efficiency (TRDrefit/TRDin)"
    ,"TRD Eloss (Max/90% quantile)"
    ,"TRD Eloss (Mean/60% quantile)"
  };
  const Int_t ngr = sizeof(name)/sizeof(Char_t*);
  if(ngr != kNgraphs){
    AliWarning("No of graphs defined different from definition");
    return 0x0;
  }

  if(!fResults){
    fResults = new TObjArray(kNgraphs);
    fResults->SetOwner();
    fResults->SetName("results");
  }

  TGraph *g = 0x0;
  if((g = dynamic_cast<TGraph*>(fResults->At(id)))){
    if(kCLEAR){ 
      for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
    } else {
      PutTrendValue(name[id], g->GetMean(2));
      PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
    }
  } else {
    if(kBUILD){
      switch(id){
      case 0:
        g = new TGraphErrors();
        break;
      case 1:
        g = new TGraphErrors();
        break;
      case 2:
        g = new TGraphErrors();
        break;
      case 3:
        g = new TGraphErrors();
        break;
      case 4:
        g = new TGraphAsymmErrors(6);
        g->SetMarkerStyle(22);g->SetMarkerColor(kRed);
        g->SetLineColor(kBlack);g->SetLineWidth(2);
        break;
      case 5:
        g = new TGraphAsymmErrors(6);
        g->SetMarkerStyle(21);
        g->SetLineColor(kRed);g->SetLineWidth(2);
        break;
      default:
        AliWarning(Form("Graph index[%d] missing/not defined.", id));
        return 0x0;
      }
      g->SetNameTitle(name[id], title[id]);
      fResults->AddAt(g, id);
    }
  }
  return g;
}

//____________________________________________________________________
void AliTRDcheckESD::Exec(Option_t *){
  //
  // Run the Analysis
  //
  if(!fESD){
    AliError("ESD event missing.");
    return;
  }

  // Get MC information if available
  AliStack * fStack = 0x0;
  if(HasMC()){
    if(!fMC){ 
      AliWarning("MC event missing");
      SetMC(kFALSE);
    } else {
      if(!(fStack = fMC->Stack())){
        AliWarning("MC stack missing");
        SetMC(kFALSE);
      }
    }
  }
  Bool_t bTRDin(0), bTRDout(0), bTRDpid(0);

  AliESDtrack *esdTrack = 0x0;
  for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
    bTRDin=0;bTRDout=0;bTRDpid=0;
    esdTrack = fESD->GetTrack(itrk);

//     if(esdTrack->GetNcls(1)) nTPC++;
//     if(esdTrack->GetNcls(2)) nTRD++;

    // track status
    ULong_t status = esdTrack->GetStatus();
    //PrintStatus(status);

    // define TPC out tracks
    if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
    if(esdTrack->GetKinkIndex(0) > 0) continue;

    // TRD PID
    Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
    // pid quality
    //esdTrack->GetTRDntrackletsPID();

    // look at external track param
    const AliExternalTrackParam *op = esdTrack->GetOuterParam();
    const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
    Double_t xyz[3];
    if(op){
      op->GetXYZ(xyz);
      op->Global2LocalPosition(xyz, op->GetAlpha());
      //printf("op @ X[%7.3f]\n", xyz[0]);
    }

    // read MC info
    if(!HasMC()) continue;

    Int_t fLabel = esdTrack->GetLabel();
    if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue; 
    
    // read MC particle
    AliMCParticle *mcParticle = 0x0; 
    if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
      AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
      continue;
    }

    AliTrackReference *ref = 0x0; 
    Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
    if(!nRefs){
      AliWarning(Form("Track refs missing. Label[%d].", fLabel));
      continue;
    }
    Int_t iref = 0;
    while(iref<nRefs){
      ref = mcParticle->GetTrackReference(iref);
      if(ref->LocalX() > fgkxTPC) break;
      ref=0x0; iref++;
    }

    // read TParticle
    //TParticle *tParticle = mcParticle->Particle(); 
    //Int_t fPdg = tParticle->GetPdgCode();
    // reject secondaries
    //if(!tParticle->IsPrimary()) continue;

    if(ref){ 
      if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
        ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
      } else {
        bTRDin=1;
        if(esdTrack->GetNcls(2)) bTRDout=1;
        if(esdTrack->GetTRDntrackletsPID()>=4) bTRDpid=1;
      }
    } else { // track stopped in TPC 
      ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
    }
    // get the MC pt !!
    Float_t pt = ref->Pt();

    TH2 *h = (TH2I*)fHistos->At(kTRDstat);
    if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
    if(/*status & AliESDtrack::k*/bTRDin) h->Fill(pt, kTRDin);
    if(/*status & AliESDtrack::k*/bTRDout){ 
      ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
      h->Fill(pt, kTRDout);
    }
    if(/*status & AliESDtrack::k*/bTRDpid) h->Fill(pt, kTRDpid);
    if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);

    if(ip){
      h = (TH2I*)fHistos->At(kTRDmom);
      Float_t p(0.);
      for(Int_t ily=6; ily--;){
        if((p=esdTrack->GetTRDmomentum(ily))<0.) continue;
        h->Fill(ip->GetP()-p, ily);
      }
    }
  }  
  PostData(0, fHistos);
}

//____________________________________________________________________
TObjArray* AliTRDcheckESD::Histos()
{
// Retrieve histograms array if already build or build it

  if(fHistos) return fHistos;

  fHistos = new TObjArray(kNhistos);
  //fHistos->SetOwner(kTRUE);
  
  TH1 *h = 0x0;

  // clusters per tracklet
  if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
    h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
    h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
    h->GetYaxis()->SetTitle("entries");
  } else h->Reset();
  fHistos->AddAt(h, kNCl);

  // status bits histogram
  if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
    h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
    h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    h->GetYaxis()->SetTitle("status bits");
    h->GetZaxis()->SetTitle("entries");
  } else h->Reset();
  fHistos->AddAt(h, kTRDstat);

  // energy loss
  if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
    h = new TH2I("hTRDmom", "TRD energy loss", 100, -1., 2., 6, -0.5, 5.5);
    h->GetXaxis()->SetTitle("p_{inner} - p_{ly} [GeV/c]");
    h->GetYaxis()->SetTitle("layer");
    h->GetZaxis()->SetTitle("entries");
  } else h->Reset();
  fHistos->AddAt(h, kTRDmom);

  return fHistos;
}

//____________________________________________________________________
Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
{
// Load data from performance file

  if(!TFile::Open(filename)){
    AliWarning(Form("Couldn't open file %s.", filename));
    return kFALSE;
  }
  TObjArray *o = 0x0;
  if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
    AliWarning("Missing histogram container.");
    return kFALSE;
  }
  fHistos = (TObjArray*)o->Clone(GetName());
  gFile->Close();
  return kTRUE;
}

//_______________________________________________________
Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
{
// Dump trending value to default file

  if(!fgFile){
    fgFile = fopen("TRD.Performance.txt", "at");
  }
  fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
  return kTRUE;
}

//____________________________________________________________________
void AliTRDcheckESD::Terminate(Option_t *)
{
// Steer post-processing 


  // geometrical efficiency
  TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
  TH1 *h1[2] = {0x0, 0x0};
  h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
  h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
  Process(h1, (TGraphErrors*)GetGraph(0));
  delete h1[0];delete h1[1];

  // tracking efficiency
  h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
  h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
  Process(h1, (TGraphErrors*)GetGraph(1));
  delete h1[1];

  // PID efficiency
  h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
  Process(h1, (TGraphErrors*)GetGraph(2));
  delete h1[1];

  // Refit efficiency
  h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
  Process(h1, (TGraphErrors*)GetGraph(3));
  delete h1[1];
  if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
 
  TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)GetGraph(4), *g09 = (TGraphAsymmErrors*)GetGraph(5);
  TAxis *ax=h2->GetXaxis();
  const Int_t nq(4);
  const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
  Double_t yq[nq];
  for(Int_t ily=6; ily--;){
    h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
    h1[0]->GetQuantiles(nq,yq,xq);
    g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
    g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
    g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
    g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);

    //printf(" max[%f] mean[%f] q[%f %f %f %f]\n", ax->GetBinCenter(h1[0]->GetMaximumBin()), h1[0]->GetMean(), yq[0], yq[1], yq[2], yq[3]);
    delete h1[0];
  }
}

//____________________________________________________________________
void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
{
// Generic function to process one reference plot

  Int_t n1 = 0, n2 = 0, ip=0;
  Double_t eff = 0.;

  TAxis *ax = h1[0]->GetXaxis();
  for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
    if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
    n2 = (Int_t)h1[1]->GetBinContent(ib);
    eff = n2/Float_t(n1);

    ip=g->GetN();
    g->SetPoint(ip, ax->GetBinCenter(ib), eff);
    g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
  }
}  

//____________________________________________________________________
void AliTRDcheckESD::PrintStatus(ULong_t status)
{
// Dump track status to stdout

  printf("ITS[i(%d) o(%d) r(%d)] TPC[i(%d) o(%d) r(%d) p(%d)] TRD[i(%d) o(%d) r(%d) p(%d) s(%d)] HMPID[o(%d) p(%d)]\n"
    ,Bool_t(status & AliESDtrack::kITSin)
    ,Bool_t(status & AliESDtrack::kITSout)
    ,Bool_t(status & AliESDtrack::kITSrefit)
    ,Bool_t(status & AliESDtrack::kTPCin)
    ,Bool_t(status & AliESDtrack::kTPCout)
    ,Bool_t(status & AliESDtrack::kTPCrefit)
    ,Bool_t(status & AliESDtrack::kTPCpid)
    ,Bool_t(status & AliESDtrack::kTRDin)
    ,Bool_t(status & AliESDtrack::kTRDout)
    ,Bool_t(status & AliESDtrack::kTRDrefit)
    ,Bool_t(status & AliESDtrack::kTRDpid)
    ,Bool_t(status & AliESDtrack::kTRDStop)
    ,Bool_t(status & AliESDtrack::kHMPIDout)
    ,Bool_t(status & AliESDtrack::kHMPIDpid)
  );
}
Commit	Line	Data
1ee39b3a	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	// Check basic detector results at ESD level
	19	// - Geometrical efficiency
	20	// - Tracking efficiency
	21	// - PID efficiency
	22	// - Refit efficiency
	23	//
	24	// Author
	25	// Alex Bercuci <A.Bercuci@gsi.de>
	26	//
	27	//////////////////////////////////////////////////////
	28
	29	#include <TClonesArray.h>
	30	#include <TObjArray.h>
	31	#include <TObject.h>
	32	#include <TH2I.h>
	33	#include <TGraphErrors.h>
	34	#include <TGraphAsymmErrors.h>
	35	#include <TFile.h>
	36	#include <TTree.h>
	37	#include <TROOT.h>
	38	#include <TChain.h>
	39	#include <TParticle.h>
	40
	41	#include "AliLog.h"
	42	#include "AliAnalysisManager.h"
	43	#include "AliESDEvent.h"
	44	#include "AliESDkink.h"
	45	#include "AliMCEvent.h"
	46	#include "AliESDInputHandler.h"
	47	#include "AliMCEventHandler.h"
	48
	49	#include "AliESDtrack.h"
	50	#include "AliMCParticle.h"
	51	#include "AliPID.h"
	52	#include "AliStack.h"
	53	#include "AliTrackReference.h"
	54
	55	#include "AliTRDcheckESD.h"
	56
	57	ClassImp(AliTRDcheckESD)
	58
	59	const Float_t AliTRDcheckESD::fgkxTPC = 290.;
	60	const Float_t AliTRDcheckESD::fgkxTOF = 365.;
	61	FILE* AliTRDcheckESD::fgFile = 0x0;
	62
	63	//____________________________________________________________________
	64	AliTRDcheckESD::AliTRDcheckESD():
65	AliAnalysisTask("checkESD", "ESD checker for TRD info")
66	,fStatus(0)
67	,fESD(0x0)
68	,fMC(0x0)
69	,fHistos(0x0)
70	,fResults(0x0)
71	{
72	//
73	// Default constructor
74	//
75	SetMC(kTRUE);
76	DefineInput(0, TChain::Class());
77	DefineOutput(0, TObjArray::Class());
78	}
79
80	//____________________________________________________________________
81	AliTRDcheckESD::~AliTRDcheckESD()
82	{
83	// Destructor
84	if(fHistos){
85	//fHistos->Delete();
86	delete fHistos;
87	}
88	if(fResults){
89	fResults->Delete();
90	delete fResults;
91	}
92	}
93
94	//____________________________________________________________________
95	void AliTRDcheckESD::ConnectInputData(Option_t *)
96	{
97	//
98	// Link the Input Data
99	//
100	TTree tree = dynamic_cast<TChain>(GetInputData(0));
101	if(tree) tree->SetBranchStatus("Tracks", 1);
102
103	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
104	fESD = esdH ? esdH->GetEvent() : 0x0;
105
106	if(!HasMC()) return;
107	AliMCEventHandler mcH = dynamic_cast<AliMCEventHandler>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
108	fMC = mcH ? mcH->MCEvent() : 0x0;
109	}
110
111	//____________________________________________________________________
112	void AliTRDcheckESD::CreateOutputObjects()
113	{
114	//
115	// Create Output Containers (TObjectArray containing 1D histograms)
116	//
117	OpenFile(0, "RECREATE");
118	Histos();
119	}
120
121	//____________________________________________________________________
122	TGraph* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
123	{
124	// Retrieve graph with "id"
125	// Possible options are :
126	// "b" - build graph if none found
127	// "c" - clear existing graph
128
129	Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
130	kCLEAR = strstr(opt, "c"); // clear existing graph
131
132	const Char_t *name[] = {
133	"Geo", "Trk", "Pid", "Ref", "Max06", "Mean09"
134	};
135	const Char_t *title[] = {
136	"TRD geometrical efficiency (TRDin/TPCout)"
137	,"TRD tracking efficiency (TRDout/TRDin)"
138	,"TRD PID efficiency (TRDpid/TRDin)"
139	,"TRD refit efficiency (TRDrefit/TRDin)"
140	,"TRD Eloss (Max/90% quantile)"
141	,"TRD Eloss (Mean/60% quantile)"
142	};
143	const Int_t ngr = sizeof(name)/sizeof(Char_t*);
144	if(ngr != kNgraphs){
145	AliWarning("No of graphs defined different from definition");
146	return 0x0;
147	}
148
149	if(!fResults){
150	fResults = new TObjArray(kNgraphs);
151	fResults->SetOwner();
152	fResults->SetName("results");
153	}
154
155	TGraph *g = 0x0;
156	if((g = dynamic_cast<TGraph*>(fResults->At(id)))){
157	if(kCLEAR){
158	for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
159	} else {
160	PutTrendValue(name[id], g->GetMean(2));
161	PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
162	}
163	} else {
164	if(kBUILD){
165	switch(id){
166	case 0:
167	g = new TGraphErrors();
168	break;
169	case 1:
170	g = new TGraphErrors();
171	break;
172	case 2:
173	g = new TGraphErrors();
174	break;
175	case 3:
176	g = new TGraphErrors();
177	break;
178	case 4:
179	g = new TGraphAsymmErrors(6);
180	g->SetMarkerStyle(22);g->SetMarkerColor(kRed);
181	g->SetLineColor(kBlack);g->SetLineWidth(2);
182	break;
183	case 5:
184	g = new TGraphAsymmErrors(6);
185	g->SetMarkerStyle(21);
186	g->SetLineColor(kRed);g->SetLineWidth(2);
187	break;
188	default:
189	AliWarning(Form("Graph index[%d] missing/not defined.", id));
190	return 0x0;
191	}
192	g->SetNameTitle(name[id], title[id]);
193	fResults->AddAt(g, id);
194	}
195	}
196	return g;
197	}
198
199	//____________________________________________________________________
200	void AliTRDcheckESD::Exec(Option_t *){
201	//
202	// Run the Analysis
203	//
204	if(!fESD){
205	AliError("ESD event missing.");
206	return;
207	}
208
209	// Get MC information if available
210	AliStack * fStack = 0x0;
211	if(HasMC()){
212	if(!fMC){
213	AliWarning("MC event missing");
214	SetMC(kFALSE);
215	} else {
216	if(!(fStack = fMC->Stack())){
217	AliWarning("MC stack missing");
218	SetMC(kFALSE);
219	}
220	}
221	}
222	Bool_t bTRDin(0), bTRDout(0), bTRDpid(0);
223
224	AliESDtrack *esdTrack = 0x0;
225	for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
226	bTRDin=0;bTRDout=0;bTRDpid=0;
227	esdTrack = fESD->GetTrack(itrk);
228
229	// if(esdTrack->GetNcls(1)) nTPC++;
230	// if(esdTrack->GetNcls(2)) nTRD++;
231
232	// track status
233	ULong_t status = esdTrack->GetStatus();
234	//PrintStatus(status);
235
236	// define TPC out tracks
237	if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
238	if(esdTrack->GetKinkIndex(0) > 0) continue;
239
240	// TRD PID
241	Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
242	// pid quality
243	//esdTrack->GetTRDntrackletsPID();
244
245	// look at external track param
246	const AliExternalTrackParam *op = esdTrack->GetOuterParam();
247	const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
248	Double_t xyz[3];
249	if(op){
250	op->GetXYZ(xyz);
251	op->Global2LocalPosition(xyz, op->GetAlpha());
252	//printf("op @ X[%7.3f]\n", xyz[0]);
253	}
254
255	// read MC info
256	if(!HasMC()) continue;
257
258	Int_t fLabel = esdTrack->GetLabel();
259	if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
260
261	// read MC particle
262	AliMCParticle *mcParticle = 0x0;
263	if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
264	AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
265	continue;
266	}
267
268	AliTrackReference *ref = 0x0;
269	Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
270	if(!nRefs){
271	AliWarning(Form("Track refs missing. Label[%d].", fLabel));
272	continue;
273	}
274	Int_t iref = 0;
275	while(iref<nRefs){
276	ref = mcParticle->GetTrackReference(iref);
277	if(ref->LocalX() > fgkxTPC) break;
278	ref=0x0; iref++;
279	}
280
281	// read TParticle
282	//TParticle *tParticle = mcParticle->Particle();
283	//Int_t fPdg = tParticle->GetPdgCode();
284	// reject secondaries
285	//if(!tParticle->IsPrimary()) continue;
286
287	if(ref){
288	if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
289	ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
290	} else {
291	bTRDin=1;
292	if(esdTrack->GetNcls(2)) bTRDout=1;
293	if(esdTrack->GetTRDntrackletsPID()>=4) bTRDpid=1;
294	}
295	} else { // track stopped in TPC
296	ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
297	}
298	// get the MC pt !!
299	Float_t pt = ref->Pt();
300
301	TH2 h = (TH2I)fHistos->At(kTRDstat);
302	if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
303	if(/status & AliESDtrack::k/bTRDin) h->Fill(pt, kTRDin);
304	if(/status & AliESDtrack::k/bTRDout){
305	((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
306	h->Fill(pt, kTRDout);
307	}
308	if(/status & AliESDtrack::k/bTRDpid) h->Fill(pt, kTRDpid);
309	if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
310
311	if(ip){
312	h = (TH2I*)fHistos->At(kTRDmom);
313	Float_t p(0.);
314	for(Int_t ily=6; ily--;){
315	if((p=esdTrack->GetTRDmomentum(ily))<0.) continue;
316	h->Fill(ip->GetP()-p, ily);
317	}
318	}
319	}
320	PostData(0, fHistos);
321	}
322
323	//____________________________________________________________________
324	TObjArray* AliTRDcheckESD::Histos()
325	{
326	// Retrieve histograms array if already build or build it
327
328	if(fHistos) return fHistos;
329
330	fHistos = new TObjArray(kNhistos);
331	//fHistos->SetOwner(kTRUE);
332
333	TH1 *h = 0x0;
334
335	// clusters per tracklet
336	if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
337	h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
338	h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
339	h->GetYaxis()->SetTitle("entries");
340	} else h->Reset();
341	fHistos->AddAt(h, kNCl);
342
343	// status bits histogram
344	if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
345	h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
346	h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
347	h->GetYaxis()->SetTitle("status bits");
348	h->GetZaxis()->SetTitle("entries");
349	} else h->Reset();
350	fHistos->AddAt(h, kTRDstat);
351
352	// energy loss
353	if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
354	h = new TH2I("hTRDmom", "TRD energy loss", 100, -1., 2., 6, -0.5, 5.5);
355	h->GetXaxis()->SetTitle("p_{inner} - p_{ly} [GeV/c]");
356	h->GetYaxis()->SetTitle("layer");
357	h->GetZaxis()->SetTitle("entries");
358	} else h->Reset();
359	fHistos->AddAt(h, kTRDmom);
360
361	return fHistos;
362	}
363
364	//____________________________________________________________________
365	Bool_t AliTRDcheckESD::Load(const Char_t filename, const Char_t name)
366	{
367	// Load data from performance file
368
369	if(!TFile::Open(filename)){
370	AliWarning(Form("Couldn't open file %s.", filename));
371	return kFALSE;
372	}
373	TObjArray *o = 0x0;
374	if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
375	AliWarning("Missing histogram container.");
376	return kFALSE;
377	}
378	fHistos = (TObjArray*)o->Clone(GetName());
379	gFile->Close();
380	return kTRUE;
381	}
382
383	//_______________________________________________________
384	Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
385	{
386	// Dump trending value to default file
387
388	if(!fgFile){
389	fgFile = fopen("TRD.Performance.txt", "at");
390	}
391	fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
392	return kTRUE;
393	}
394
395	//____________________________________________________________________
396	void AliTRDcheckESD::Terminate(Option_t *)
397	{
398	// Steer post-processing
399
400
401	// geometrical efficiency
402	TH2I h2 = (TH2I)fHistos->At(kTRDstat);
403	TH1 *h1[2] = {0x0, 0x0};
404	h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
405	h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
406	Process(h1, (TGraphErrors*)GetGraph(0));
407	delete h1[0];delete h1[1];
408
409	// tracking efficiency
410	h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
411	h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
412	Process(h1, (TGraphErrors*)GetGraph(1));
413	delete h1[1];
414
415	// PID efficiency
416	h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
417	Process(h1, (TGraphErrors*)GetGraph(2));
418	delete h1[1];
419
420	// Refit efficiency
421	h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
422	Process(h1, (TGraphErrors*)GetGraph(3));
423	delete h1[1];
424	if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
425
426	TGraphAsymmErrors g06 = (TGraphAsymmErrors)GetGraph(4), g09 = (TGraphAsymmErrors)GetGraph(5);
427	TAxis *ax=h2->GetXaxis();
428	const Int_t nq(4);
429	const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
430	Double_t yq[nq];
431	for(Int_t ily=6; ily--;){
432	h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
433	h1[0]->GetQuantiles(nq,yq,xq);
434	g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
435	g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
436	g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
437	g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
438
439	//printf(" max[%f] mean[%f] q[%f %f %f %f]\n", ax->GetBinCenter(h1[0]->GetMaximumBin()), h1[0]->GetMean(), yq[0], yq[1], yq[2], yq[3]);
440	delete h1[0];
441	}
442	}
443
444	//____________________________________________________________________
445	void AliTRDcheckESD::Process(TH1 *h1, TGraphErrors g)
446	{
447	// Generic function to process one reference plot
448
449	Int_t n1 = 0, n2 = 0, ip=0;
450	Double_t eff = 0.;
451
452	TAxis *ax = h1[0]->GetXaxis();
453	for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
454	if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
455	n2 = (Int_t)h1[1]->GetBinContent(ib);
456	eff = n2/Float_t(n1);
457
458	ip=g->GetN();
459	g->SetPoint(ip, ax->GetBinCenter(ib), eff);
460	g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
461	}
462	}
463
464	//____________________________________________________________________
465	void AliTRDcheckESD::PrintStatus(ULong_t status)
466	{
467	// Dump track status to stdout
468
469	printf("ITS[i(%d) o(%d) r(%d)] TPC[i(%d) o(%d) r(%d) p(%d)] TRD[i(%d) o(%d) r(%d) p(%d) s(%d)] HMPID[o(%d) p(%d)]\n"
470	,Bool_t(status & AliESDtrack::kITSin)
471	,Bool_t(status & AliESDtrack::kITSout)
472	,Bool_t(status & AliESDtrack::kITSrefit)
473	,Bool_t(status & AliESDtrack::kTPCin)
474	,Bool_t(status & AliESDtrack::kTPCout)
475	,Bool_t(status & AliESDtrack::kTPCrefit)
476	,Bool_t(status & AliESDtrack::kTPCpid)
477	,Bool_t(status & AliESDtrack::kTRDin)
478	,Bool_t(status & AliESDtrack::kTRDout)
479	,Bool_t(status & AliESDtrack::kTRDrefit)
480	,Bool_t(status & AliESDtrack::kTRDpid)
481	,Bool_t(status & AliESDtrack::kTRDStop)
482	,Bool_t(status & AliESDtrack::kHMPIDout)
483	,Bool_t(status & AliESDtrack::kHMPIDpid)
484	);
485	}
486