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


#include <TClonesArray.h>
#include <TObjArray.h>
#include <TObject.h>
#include <TH2I.h>
#include <TGraphErrors.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 Int_t   AliTRDcheckESD::fgkNgraphs = 4;
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)
{
  //
  // Default constructor
  //
  SetMC(kTRUE);
  DefineInput(0, TChain::Class());
  DefineOutput(0, TObjArray::Class());
}

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

//____________________________________________________________________
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();
}

//____________________________________________________________________
TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
{
  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"
  };
  const Char_t *title[] = {
    "TRD geometrical efficiency (TRDin/TPCout)"
    ,"TRD tracking efficiency (TRDout/TRDin)"
    ,"TRD PID efficiency (TRDpid/TRDin)"
    ,"TRD refit efficiency (TRDrefit/TRDin)"
  };
  const Int_t ngr = sizeof(name)/sizeof(Char_t*);
  if(ngr != fgkNgraphs){
    AliWarning("No of graphs defined different from definition");
    return 0x0;
  }

  TObjArray *res = 0x0;
  if(!(res = (TObjArray*)fHistos->At(kResults)) ||
      (id < 0 || id >= ngr)){
    AliWarning("Graph array missing.");
    return 0x0;
  }

  TGraphErrors *g = 0x0;
  if((g = dynamic_cast<TGraphErrors*>(res->At(id)))){
    if(kCLEAR){ 
      for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
    } else {
      PutTrendValue(name[id], g->GetMean(2), g->GetRMS(2));
    }
  } else {
    if(kBUILD){
      g = new TGraphErrors();
      g->SetNameTitle(name[id], title[id]);
      res->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 TRDin(0), TRDout(0), TRDpid(0);

  AliESDtrack *esdTrack = 0x0;
  for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
    TRDin=0;TRDout=0;TRDpid=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();
    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 {
        TRDin=1;
        if(esdTrack->GetNcls(2)) TRDout=1;
        if(esdTrack->GetTRDntrackletsPID()>=4) TRDpid=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*/TRDin) h->Fill(pt, kTRDin);
    if(/*status & AliESDtrack::k*/TRDout){ 
      ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
      h->Fill(pt, kTRDout);
    }
    if(/*status & AliESDtrack::k*/TRDpid) h->Fill(pt, kTRDpid);
    if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
  }  
  PostData(0, fHistos);
}

//____________________________________________________________________
TObjArray* AliTRDcheckESD::Histos()
{
  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);

  // results array
  TObjArray *res = new TObjArray();
  res->SetName("Results");
  fHistos->AddAt(res, kResults);
  return fHistos;
}

//____________________________________________________________________
Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
{
  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, Double_t err)
{
  if(!fgFile){
    fgFile = fopen("TRD.Performance.txt", "at");
  }
  fprintf(fgFile, "%s_%s %f %f\n", GetName(), name, val, err);
  return kTRUE;
}

//____________________________________________________________________
void AliTRDcheckESD::Terminate(Option_t *)
{
  TObjArray *res = 0x0;
  if(!(res = (TObjArray*)fHistos->At(kResults))){
    AliWarning("Graph container missing.");
    return;
  }
  if(!res->GetEntriesFast()) res->Expand(fgkNgraphs);
  
  // geometrical efficiency
  TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
  TH1 *h1[2] = {0x0, 0x0};
  h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
  h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
  Process(h1, GetGraph(0));
  delete h1[0];delete h1[1];

  // tracking efficiency
  h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
  h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
  Process(h1, GetGraph(1));
  delete h1[1];

  // PID efficiency
  //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
  h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
  Process(h1, GetGraph(2));
  delete h1[1];

  // Refit efficiency
  //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
  h1[1] = h2->ProjectionX("px1", kTRDref, kTRDref);
  Process(h1, GetGraph(3));
  delete h1[1];
}

//____________________________________________________________________
void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
{
  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)
{
  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
9a281e83	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	#include <TClonesArray.h>
	18	#include <TObjArray.h>
	19	#include <TObject.h>
	20	#include <TH2I.h>
48797466	21	#include <TGraphErrors.h>
9a281e83	22	#include <TFile.h>
	23	#include <TTree.h>
	24	#include <TROOT.h>
	25	#include <TChain.h>
	26	#include <TParticle.h>
	27
	28	#include "AliLog.h"
	29	#include "AliAnalysisManager.h"
	30	#include "AliESDEvent.h"
d12237d6	31	#include "AliESDkink.h"
9a281e83	32	#include "AliMCEvent.h"
	33	#include "AliESDInputHandler.h"
	34	#include "AliMCEventHandler.h"
	35
	36	#include "AliESDtrack.h"
	37	#include "AliMCParticle.h"
	38	#include "AliPID.h"
	39	#include "AliStack.h"
	40	#include "AliTrackReference.h"
9a281e83	41
	42	#include "AliTRDcheckESD.h"
	43
	44	ClassImp(AliTRDcheckESD)
	45
965e895b	46	const Int_t AliTRDcheckESD::fgkNgraphs = 4;
	47	const Float_t AliTRDcheckESD::fgkxTPC = 290.;
	48	const Float_t AliTRDcheckESD::fgkxTOF = 365.;
59589fc9	49	FILE* AliTRDcheckESD::fgFile = 0x0;
9a281e83	50
	51	//____________________________________________________________________
	52	AliTRDcheckESD::AliTRDcheckESD():
6da3eee3	53	AliAnalysisTask("checkESD", "ESD checker for TRD info")
9a281e83	54	,fStatus(0)
	55	,fESD(0x0)
	56	,fMC(0x0)
	57	,fHistos(0x0)
	58	{
	59	//
	60	// Default constructor
	61	//
49c7f0a5	62	SetMC(kTRUE);
9a281e83	63	DefineInput(0, TChain::Class());
	64	DefineOutput(0, TObjArray::Class());
	65	}
	66
	67	//____________________________________________________________________
	68	AliTRDcheckESD::~AliTRDcheckESD()
	69	{
	70	if(fHistos){
	71	//fHistos->Delete();
	72	delete fHistos;
	73	}
	74	}
	75
	76	//____________________________________________________________________
	77	void AliTRDcheckESD::ConnectInputData(Option_t *)
	78	{
	79	//
	80	// Link the Input Data
	81	//
	82	TTree tree = dynamic_cast<TChain>(GetInputData(0));
	83	if(tree) tree->SetBranchStatus("Tracks", 1);
	84
	85	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	86	fESD = esdH ? esdH->GetEvent() : 0x0;
	87
	88	if(!HasMC()) return;
	89	AliMCEventHandler mcH = dynamic_cast<AliMCEventHandler>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
	90	fMC = mcH ? mcH->MCEvent() : 0x0;
	91	}
	92
	93	//____________________________________________________________________
	94	void AliTRDcheckESD::CreateOutputObjects()
	95	{
	96	//
	97	// Create Output Containers (TObjectArray containing 1D histograms)
	98	//
	99	OpenFile(0, "RECREATE");
c2fd5953	100	Histos();
9a281e83	101	}
9a281e83	102
965e895b	103	//____________________________________________________________________
ea4502f6	104	TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
965e895b	105	{
ea4502f6	106	Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
ea4502f6	107	kCLEAR = strstr(opt, "c"); // clear existing graph
965e895b	108
	109	const Char_t *name[] = {
	110	"Geo", "Trk", "Pid", "Ref"
	111	};
	112	const Char_t *title[] = {
	113	"TRD geometrical efficiency (TRDin/TPCout)"
	114	,"TRD tracking efficiency (TRDout/TRDin)"
	115	,"TRD PID efficiency (TRDpid/TRDin)"
	116	,"TRD refit efficiency (TRDrefit/TRDin)"
	117	};
	118	const Int_t ngr = sizeof(name)/sizeof(Char_t*);
	119	if(ngr != fgkNgraphs){
	120	AliWarning("No of graphs defined different from definition");
	121	return 0x0;
	122	}
	123
	124	TObjArray *res = 0x0;
	125	if(!(res = (TObjArray*)fHistos->At(kResults)) \|\|
	126	(id < 0 \|\| id >= ngr)){
ea4502f6	127	AliWarning("Graph array missing.");
965e895b	128	return 0x0;
	129	}
	130
	131	TGraphErrors *g = 0x0;
	132	if((g = dynamic_cast<TGraphErrors*>(res->At(id)))){
59589fc9	133	if(kCLEAR){
	134	for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
	135	} else {
	136	PutTrendValue(name[id], g->GetMean(2), g->GetRMS(2));
	137	}
965e895b	138	} else {
	139	if(kBUILD){
	140	g = new TGraphErrors();
	141	g->SetNameTitle(name[id], title[id]);
	142	res->AddAt(g, id);
	143	}
	144	}
	145	return g;
	146	}
	147
9a281e83	148	//____________________________________________________________________
	149	void AliTRDcheckESD::Exec(Option_t *){
	150	//
	151	// Run the Analysis
	152	//
	153	if(!fESD){
48797466	154	AliError("ESD event missing.");
9a281e83	155	return;
	156	}
	157
	158	// Get MC information if available
	159	AliStack * fStack = 0x0;
2b104c26	160	if(HasMC()){
	161	if(!fMC){
	162	AliWarning("MC event missing");
9a281e83	163	SetMC(kFALSE);
2b104c26	164	} else {
	165	if(!(fStack = fMC->Stack())){
	166	AliWarning("MC stack missing");
	167	SetMC(kFALSE);
	168	}
9a281e83	169	}
9a281e83	170	}
d12237d6	171	Bool_t TRDin(0), TRDout(0), TRDpid(0);
9a281e83	172
9a281e83	173	AliESDtrack *esdTrack = 0x0;
9a281e83	174	for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
d12237d6	175	TRDin=0;TRDout=0;TRDpid=0;
9a281e83	176	esdTrack = fESD->GetTrack(itrk);
d12237d6	177
	178	// if(esdTrack->GetNcls(1)) nTPC++;
	179	// if(esdTrack->GetNcls(2)) nTRD++;
9a281e83	180
9a281e83	181	// track status
d12237d6	182	ULong_t status = esdTrack->GetStatus();
2b104c26	183	//PrintStatus(status);
d12237d6	184
	185	// define TPC out tracks
	186	if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
	187	if(esdTrack->GetKinkIndex(0) > 0) continue;
9a281e83	188
	189	// TRD PID
	190	Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
	191	// pid quality
ed15ef4f	192	//esdTrack->GetTRDntrackletsPID();
9a281e83	193
	194	// look at external track param
	195	const AliExternalTrackParam *op = esdTrack->GetOuterParam();
	196	Double_t xyz[3];
	197	if(op){
	198	op->GetXYZ(xyz);
	199	op->Global2LocalPosition(xyz, op->GetAlpha());
	200	//printf("op @ X[%7.3f]\n", xyz[0]);
	201	}
	202
	203	// read MC info
	204	if(!HasMC()) continue;
	205
	206	Int_t fLabel = esdTrack->GetLabel();
	207	if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
	208
	209	// read MC particle
	210	AliMCParticle *mcParticle = 0x0;
7aad0c47	211	if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
48797466	212	AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
9a281e83	213	continue;
	214	}
	215
	216	AliTrackReference *ref = 0x0;
	217	Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
48797466	218	if(!nRefs){
	219	AliWarning(Form("Track refs missing. Label[%d].", fLabel));
	220	continue;
	221	}
9a281e83	222	Int_t iref = 0;
	223	while(iref<nRefs){
	224	ref = mcParticle->GetTrackReference(iref);
965e895b	225	if(ref->LocalX() > fgkxTPC) break;
9a281e83	226	ref=0x0; iref++;
	227	}
	228
	229	// read TParticle
965e895b	230	//TParticle *tParticle = mcParticle->Particle();
48797466	231	//Int_t fPdg = tParticle->GetPdgCode();
d12237d6	232	// reject secondaries
d12237d6	233	//if(!tParticle->IsPrimary()) continue;
9a281e83	234
d12237d6	235	if(ref){
965e895b	236	if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
48797466	237	ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
9a281e83	238	} else {
d12237d6	239	TRDin=1;
d12237d6	240	if(esdTrack->GetNcls(2)) TRDout=1;
ed15ef4f	241	if(esdTrack->GetTRDntrackletsPID()>=4) TRDpid=1;
9a281e83	242	}
d12237d6	243	} else { // track stopped in TPC
48797466	244	ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
48797466	245	}
d12237d6	246	// get the MC pt !!
9a281e83	247	Float_t pt = ref->Pt();
d12237d6	248
9a281e83	249	TH2 h = (TH2I)fHistos->At(kTRDstat);
d12237d6	250	if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
d12237d6	251	if(/status & AliESDtrack::k/TRDin) h->Fill(pt, kTRDin);
9a281e83	252	if(/status & AliESDtrack::k/TRDout){
9a281e83	253	((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
d12237d6	254	h->Fill(pt, kTRDout);
9a281e83	255	}
965e895b	256	if(/status & AliESDtrack::k/TRDpid) h->Fill(pt, kTRDpid);
965e895b	257	if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
9a281e83	258	}
9a281e83	259	PostData(0, fHistos);
	260	}
	261
c2fd5953	262	//____________________________________________________________________
	263	TObjArray* AliTRDcheckESD::Histos()
	264	{
	265	if(fHistos) return fHistos;
	266
	267	fHistos = new TObjArray(kNhistos);
	268	//fHistos->SetOwner(kTRUE);
	269
	270	TH1 *h = 0x0;
	271
	272	// clusters per tracklet
	273	if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
	274	h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
	275	h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
	276	h->GetYaxis()->SetTitle("entries");
	277	} else h->Reset();
	278	fHistos->AddAt(h, kNCl);
	279
	280	// status bits histogram
	281	if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
	282	h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
	283	h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	284	h->GetYaxis()->SetTitle("status bits");
	285	h->GetZaxis()->SetTitle("entries");
	286	} else h->Reset();
	287	fHistos->AddAt(h, kTRDstat);
	288
	289	// results array
	290	TObjArray *res = new TObjArray();
	291	res->SetName("Results");
	292	fHistos->AddAt(res, kResults);
	293	return fHistos;
	294	}
	295
965e895b	296	//____________________________________________________________________
	297	Bool_t AliTRDcheckESD::Load(const Char_t filename, const Char_t name)
	298	{
	299	if(!TFile::Open(filename)){
	300	AliWarning(Form("Couldn't open file %s.", filename));
	301	return kFALSE;
	302	}
	303	TObjArray *o = 0x0;
	304	if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
	305	AliWarning("Missing histogram container.");
	306	return kFALSE;
	307	}
	308	fHistos = (TObjArray*)o->Clone(GetName());
	309	gFile->Close();
	310	return kTRUE;
	311	}
9a281e83	312
59589fc9	313	//_______________________________________________________
	314	Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val, Double_t err)
	315	{
	316	if(!fgFile){
	317	fgFile = fopen("TRD.Performance.txt", "at");
	318	}
	319	fprintf(fgFile, "%s_%s %f %f\n", GetName(), name, val, err);
	320	return kTRUE;
	321	}
	322
9a281e83	323	//____________________________________________________________________
	324	void AliTRDcheckESD::Terminate(Option_t *)
	325	{
965e895b	326	TObjArray *res = 0x0;
	327	if(!(res = (TObjArray*)fHistos->At(kResults))){
	328	AliWarning("Graph container missing.");
	329	return;
	330	}
	331	if(!res->GetEntriesFast()) res->Expand(fgkNgraphs);
48797466	332
48797466	333	// geometrical efficiency
965e895b	334	TH2I h2 = (TH2I)fHistos->At(kTRDstat);
965e895b	335	TH1 *h1[2] = {0x0, 0x0};
d12237d6	336	h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
d12237d6	337	h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
965e895b	338	Process(h1, GetGraph(0));
ea4502f6	339	delete h1[0];delete h1[1];
48797466	340
48797466	341	// tracking efficiency
d12237d6	342	h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
d12237d6	343	h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
965e895b	344	Process(h1, GetGraph(1));
ea4502f6	345	delete h1[1];
48797466	346
48797466	347	// PID efficiency
965e895b	348	//h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
d12237d6	349	h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
965e895b	350	Process(h1, GetGraph(2));
ea4502f6	351	delete h1[1];
965e895b	352
	353	// Refit efficiency
	354	//h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
	355	h1[1] = h2->ProjectionX("px1", kTRDref, kTRDref);
	356	Process(h1, GetGraph(3));
ea4502f6	357	delete h1[1];
965e895b	358	}
	359
	360	//____________________________________________________________________
	361	void AliTRDcheckESD::Process(TH1 *h1, TGraphErrors g)
	362	{
	363	Int_t n1 = 0, n2 = 0, ip=0;
	364	Double_t eff = 0.;
	365
	366	TAxis *ax = h1[0]->GetXaxis();
48797466	367	for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
	368	if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
	369	n2 = (Int_t)h1[1]->GetBinContent(ib);
	370	eff = n2/Float_t(n1);
	371
	372	ip=g->GetN();
	373	g->SetPoint(ip, ax->GetBinCenter(ib), eff);
	374	g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
	375	}
2b104c26	376	}
	377
	378	//____________________________________________________________________
	379	void AliTRDcheckESD::PrintStatus(ULong_t status)
	380	{
	381	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"
	382	,Bool_t(status & AliESDtrack::kITSin)
	383	,Bool_t(status & AliESDtrack::kITSout)
	384	,Bool_t(status & AliESDtrack::kITSrefit)
	385	,Bool_t(status & AliESDtrack::kTPCin)
	386	,Bool_t(status & AliESDtrack::kTPCout)
	387	,Bool_t(status & AliESDtrack::kTPCrefit)
	388	,Bool_t(status & AliESDtrack::kTPCpid)
	389	,Bool_t(status & AliESDtrack::kTRDin)
	390	,Bool_t(status & AliESDtrack::kTRDout)
	391	,Bool_t(status & AliESDtrack::kTRDrefit)
	392	,Bool_t(status & AliESDtrack::kTRDpid)
	393	,Bool_t(status & AliESDtrack::kTRDStop)
	394	,Bool_t(status & AliESDtrack::kHMPIDout)
	395	,Bool_t(status & AliESDtrack::kHMPIDpid)
	396	);
9a281e83	397	}
2b104c26	398