[u/mrichter/AliRoot.git] / PWGPP / TRD / TRDPID / AliTRDPIDTree.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>             *
* 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.                  *
**************************************************************************/
//
// The task:
// stores TRD PID quantities in a Tree
// output can then be used for e.g. parameter creation
//
// Author:
// Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>
//


#include "AliTRDPIDTree.h"
#include "TChain.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliPID.h"
#include "AliESDtrack.h"
#include "AliPIDResponse.h"
#include "AliInputEventHandler.h"
#include "AliESDInputHandler.h"
#include "AliESDv0KineCuts.h"
#include "AliESDv0.h"
#include "AliCentrality.h"
#include "AliTRDgeometry.h"
#include "TTree.h"

class TCanvas;
class TAxis;
class TFile;
class TStyle;
class TString;
class TH1F;
class TH2D;
class THnSparse;
class TLegend;
class TVirtualFitter;
class AliESDtrackCuts;
class AliStack;
class AliMCParticle;


using namespace std;

ClassImp(AliTRDPIDTree)

Double_t AliTRDPIDTree::fgMinLayer = 1.;

//________________________________________________________________________
AliTRDPIDTree::AliTRDPIDTree(const char *name)
    : AliAnalysisTaskSE(name), fV0tags(0x0), fV0cuts(0x0), fV0electrons(0x0), fV0pions(0x0), fV0protons(0x0),
    fESDEvent(0), fMCEvent(0), fMCStack(0), fTreeTRDPID(0), fPIDResponse(0), fOutputContainer(0), fESDtrackCuts(0),
    fESDtrackCutsV0(0), fListQATRD(0x0), fListQATRDV0(0x0),
    fNumTagsStored(0), fCollisionSystem(3),
    fpdg(0), frun(0), frunnumber(0), fcentrality(0), fTRDNtracklets(0), fTRDNcls(0), fTRDntracklets(0), fTRDntrackletsPID(0),
    fTRDtheta(0), fTRDsignal(0), fTRDnclsdEdx(0), fTRDnch(0), fPDG(0), fPDGTRUE(0), fChi2(0), fhtrackCuts(0), fhArmenteros(0)
{

  //
  // Constructor
  //
  DefineInput(0, TChain::Class());
  DefineOutput(1, TTree::Class());
  DefineOutput(2, TList::Class());

}


//_________________________________________________
AliTRDPIDTree::~AliTRDPIDTree()
{

  //
  // Destructor
  //

    delete fV0cuts;
    delete fV0electrons;
    delete fV0pions;
    delete fV0protons;
    delete fV0tags;
    fV0tags = 0;
    fNumTagsStored = 0;
}


//________________________________________________________________________
void AliTRDPIDTree::UserCreateOutputObjects()
{
    //
    // Definition of user output ntuple and histogram file
    //


    AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  if (!inputHandler)
    printf("Inputhandler not available \n");
  else
    fPIDResponse = inputHandler->GetPIDResponse();

   // V0 Kine cuts 
  fV0cuts = new AliESDv0KineCuts();

  // V0 PID Obj arrays
  fV0electrons = new TObjArray;
  fV0pions     = new TObjArray;
  fV0protons   = new TObjArray;

 
  OpenFile(1);
    fTreeTRDPID = new TTree("TreeTRDPID", "TRD PID");
    fTreeTRDPID->Branch("run", &frunnumber);
    fTreeTRDPID->Branch("centrality", &fcentrality);
    fTreeTRDPID->Branch("TRDNtracklets",&fTRDNtracklets);
    fTreeTRDPID->Branch("TRDslices[48]",fTRDslices);
    fTreeTRDPID->Branch("TRDMomentum[6]",fTRDMomentum);
    fTreeTRDPID->Branch("TRDNcls",&fTRDNcls);
    fTreeTRDPID->Branch("TRDntracklets",&fTRDntracklets);
    fTreeTRDPID->Branch("TRDntrackletsPID",&fTRDntrackletsPID);
    fTreeTRDPID->Branch("TRDphi[6]",fTRDphi);
    fTreeTRDPID->Branch("TRDY[6]",fTRDY);
    fTreeTRDPID->Branch("TRDtheta",&fTRDtheta);
    fTreeTRDPID->Branch("TRDsignal",&fTRDsignal);
    fTreeTRDPID->Branch("TRDnclsdEdx",&fTRDnclsdEdx);
    fTreeTRDPID->Branch("TRDnch",&fTRDnch);
    fTreeTRDPID->Branch("NSigmaTPC[3]",fNSigmaTPC);
    fTreeTRDPID->Branch("NSigmaTOF[3]",fNSigmaTOF);
    fTreeTRDPID->Branch("PDG",&fPDG);
    fTreeTRDPID->Branch("PDGTRUE",&fPDGTRUE);
    fTreeTRDPID->Branch("DCA[2]",fDCA);
    fTreeTRDPID->Branch("Chi2",&fChi2);
    fTreeTRDPID->Branch("TRDsigma[5]",fsigmaTRD);
    fTreeTRDPID->Branch("TRDdelta[5]",fdeltaTRD);
    fTreeTRDPID->Branch("TRDratio[5]",fratioTRD);

    PostData(1,fTreeTRDPID);


  OpenFile(2);
  fListQATRD=new TList;
  fListQATRD->SetOwner();
  fListQATRDV0=new TList;
  fListQATRDV0->SetOwner();
  fListQATRDV0->SetName("V0decay");
  fListQATRD->Add(fListQATRDV0);

  SetupV0qa();

  fhtrackCuts  = new TH1F("fhtrackCuts","TrackEventCuts QA",10,-0.5,9.5);
  fListQATRD->Add(fhtrackCuts);

  PostData(2,fListQATRD);


}


//_____________________________________________________________________________
void AliTRDPIDTree::UserExec(Option_t *)
{
    //
    //calls the Process function
    //

    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());

    if (!esdH) {
      printf("ERROR: Could not get ESDInputHandler \n");
    }
    else fESDEvent = esdH->GetEvent();
    
    // If MC, forward MCevent
    fMCEvent = dynamic_cast<AliMCEvent*>(MCEvent());
    fMCStack=NULL;
    if(fMCEvent){
	fMCStack = fMCEvent->Stack();
    }

    FillV0PIDlist();
    //
    Process(fESDEvent, fMCEvent);
    //

    PostData(1,fTreeTRDPID);
    PostData(2,fListQATRD);

    // Clear the V0 PID arrays
    ClearV0PIDlist();

}


//________________________________________________________________________
void AliTRDPIDTree::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
{
    //
    //called for each event
    //

  if (!esdEvent) {
    Printf("ERROR: esdEvent not available"); 
    return;
  }

  if (!mcEvent) {
//    Printf("ERROR: mcEvent not available");
//    return;
  }

  if (!fPIDResponse ) {
    Printf("ERROR: No PIDresponse and/or v0KineCuts!");
    return;
  }

  Float_t centralityFper=99;

  AliCentrality *esdCentrality = esdEvent->GetCentrality();
  centralityFper = esdCentrality->GetCentralityPercentile("V0M");


  const AliESDVertex* fESDEventvertex = esdEvent->GetPrimaryVertexTracks(); 
  if (!fESDEventvertex)
    return;
  
  Int_t ncontr = fESDEventvertex->GetNContributors();

  if (ncontr <= 0) return;
  frunnumber = fESDEvent->GetRunNumber();

  // - Begin: track loop for electrons from V0 -
    for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
    //  AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack);
      AliESDtrack *track=(AliESDtrack*)fV0electrons->At(itrack);

      fpdg=11;
      FillTree(track, fpdg, frun, centralityFper);

    } // - End: track loop for electrons from V0 -

  // - Begin: track loop for pions from V0 -
    for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
	//       AliVTrack *track=(AliVTrack*)fV0pions->At(itrack);
	AliESDtrack *track=(AliESDtrack*)fV0pions->At(itrack);
	fpdg=211;
	FillTree(track, fpdg, frun, centralityFper);
    }

    // - Begin: track loop for protons from V0 -
    for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
	//       AliVTrack *track=(AliVTrack*)fV0protons->At(itrack);
	AliESDtrack *track=(AliESDtrack*)fV0protons->At(itrack);
	fpdg=2212;
	FillTree(track, fpdg, frun, centralityFper);
    }

   
    PostData(1,fTreeTRDPID);
    PostData(2,fListQATRD);

}      

//________________________________________________________________________
void AliTRDPIDTree::FillTree(AliESDtrack *track, Int_t pdgfromv0, Int_t runnumber, Int_t centralityvalue)
{
    //
    // Fill tree with track and event variables for e, pi, proton from V0 candidates
    //

    // global event and track cuts
    if(!PassTrackCuts(track)) return;

    // track cuts:
    Int_t ntrackletstracking=track->GetTRDntracklets();
    if(ntrackletstracking<fgMinLayer) return;
    Int_t trdnch  = track->GetTRDNchamberdEdx(); // only available for p-Pb up to now
//    if(trdnch<1) return;


    // check if TOF can provide PID (if not true return)
    UInt_t status = track->GetStatus();
    Bool_t hasTOFout  = status&AliESDtrack::kTOFout;
    Bool_t hasTOFtime = status&AliESDtrack::kTIME;
    if(!(hasTOFout && hasTOFtime)) return;

    // add kaons?
    Double_t nSigmaTPC = 0;
    Double_t nSigmaTOF = 0;

    if(TMath::Abs(pdgfromv0)==11){
	nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron);
	nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron);
    }
    if(TMath::Abs(pdgfromv0)==211){
	nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
	nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
    }
    if(TMath::Abs(pdgfromv0)==2212){
	nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
	nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);
    }

    if(nSigmaTPC>3||nSigmaTPC<-3){
	return;
    }

    if(nSigmaTOF>3||nSigmaTOF<-3){
	return;
    }

    frun=runnumber;

    Int_t charge = track->Charge();
    fPDG=pdgfromv0*charge;
    fcentrality = centralityvalue;
    fTRDntracklets = ntrackletstracking;
    fTRDntrackletsPID = track->GetTRDntrackletsPID();

    fTRDNcls=track->GetTRDncls();
    fTRDtheta=track->Theta()-0.5*TMath::Pi();
    fChi2=track->GetTRDchi2();

    fTRDsignal=track->GetTRDsignal(); // truncated mean signal
    fTRDnclsdEdx=track->GetTRDNclusterdEdx(); // number of clusters dedx
    fTRDnch=trdnch;  // number of chambers dedx

    track->GetImpactParameters(fDCA[0],fDCA[1]);

    fNSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron);
    fNSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
    fNSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
    fNSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron);
    fNSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
    fNSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);

    Int_t ntracklets=0;
    for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){

	fTRDphi[iPl]=GetPhi(track,iPl,fTRDY[iPl]);
	Float_t dEdx=0;
	fTRDMomentum[iPl]= track->GetTRDmomentum(iPl);
	for(int isl= 0; isl<= 7;isl++){
	    fTRDslices[iPl*8+isl]=track->GetTRDslice(iPl,isl);
	    Double_t sig=track->GetTRDslice(iPl,isl);
	    if(sig>0){
		dEdx+=sig;
	    }
	}
	if(dEdx>0){
	    ntracklets++;
	}
    }
    fTRDNtracklets=ntracklets;

    fPDGTRUE=0;
    if(fMCStack){
	Int_t label=track->GetLabel();
	if(label>=0){
	    TParticle *part=fMCStack->Particle(label);
	    if(part){
		fPDGTRUE=part->GetPdgCode();
		if(TMath::Abs(fPDGTRUE)==11)fPDGTRUE=-fPDGTRUE; // switch signs (particle code says 11 for e- and -11 for e+)
	    }
	}
    }

    AliPID::EParticleType types[]={AliPID::kElectron, AliPID::kMuon, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
    for(Int_t itrdpid=0; itrdpid<5; itrdpid++){
	fsigmaTRD[itrdpid]= fPIDResponse->NumberOfSigmas(AliPIDResponse::kTRD, track, types[itrdpid]); //  (signal-expsig)/(sigma), where sigma = mean*resolution
	fdeltaTRD[itrdpid]= fPIDResponse->GetSignalDelta(AliPIDResponse::kTRD, track, types[itrdpid]);   // signal/(expsig)
	fratioTRD[itrdpid]= fPIDResponse->GetSignalDelta(AliPIDResponse::kTRD, track, types[itrdpid], kTRUE);  // difference signal to theor. value divided by resolution
    }


    fTreeTRDPID->Fill();

}

//________________________________________________________________________
Int_t  AliTRDPIDTree::CompareFloat(Float_t f1, Float_t f2) const
{
    //
    //compares if the Float_t f1 is equal with f2 and returns 1 if true and 0 if false
    //
    Float_t precision = 0.00001;
    if (((f1 - precision) < f2) &&
	((f1 + precision) > f2))
    {
	return 1;
    }
    else
    {
	return 0;
    }
}


//________________________________________________________________________
void AliTRDPIDTree::Terminate(const Option_t *)
{
    //
    // Terminate function
    //


}


//______________________________________________________________________________
void AliTRDPIDTree::FillV0PIDlist(){

  //
  // Fill the PID object arrays holding the pointers to identified particle tracks
  //

  // Dynamic cast to ESD events (DO NOTHING for AOD events)
  AliESDEvent *event = dynamic_cast<AliESDEvent *>(InputEvent());
  if ( !event )  return;

  if(IsPbPb()) {
      fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPbPb);
  }
  else {
      fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPP);
  }


  // V0 selection
  // set event
  fV0cuts->SetEvent(event);

  const Int_t numTracks = event->GetNumberOfTracks();
  fV0tags = new Int_t[numTracks];
  for (Int_t i = 0; i < numTracks; i++)
    fV0tags[i] = 0;

  fNumTagsStored = numTracks;

  // loop over V0 particles
  for(Int_t iv0=0; iv0<event->GetNumberOfV0s();iv0++){


    AliESDv0 *v0 = (AliESDv0 *) event->GetV0(iv0);
 
    if(!v0) continue;
    if(v0->GetOnFlyStatus()) continue; 
    // Get the particle selection 
    Bool_t foundV0 = kFALSE;
    Int_t pdgV0, pdgP, pdgN;
    foundV0 = fV0cuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
    if(!foundV0) continue;
    Int_t iTrackP = v0->GetPindex();  // positive track
    Int_t iTrackN = v0->GetNindex();  // negative track

    // v0 Armenteros plot (QA)
    Float_t armVar[2] = {0.0,0.0};
    fV0cuts->Armenteros(v0, armVar);

    if(fListQATRDV0&&fhArmenteros)fhArmenteros->Fill(armVar[0],armVar[1]);

    // fill the Object arrays
    // positive particles
    if( pdgP == -11){
	fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackP));
        fV0tags[iTrackP] = 11;
    }
    else if( pdgP == 211){
	fV0pions->Add((AliVTrack*)event->GetTrack(iTrackP));
        fV0tags[iTrackP] = 211;
    }
    else if( pdgP == 2212){
	fV0protons->Add((AliVTrack*)event->GetTrack(iTrackP));
	fV0tags[iTrackP] = 2212;
    }

    // negative particles
    if( pdgN == 11){
	fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackN));
        fV0tags[iTrackN] = -11;
    }
    else if( pdgN == -211){
	fV0pions->Add((AliVTrack*)event->GetTrack(iTrackN));
        fV0tags[iTrackN] = -211;
    }
    else if( pdgN == -2212){
	fV0protons->Add((AliVTrack*)event->GetTrack(iTrackN));
        fV0tags[iTrackN] = -2212;
    }


  }
}

//______________________________________________________________________________
Int_t AliTRDPIDTree::GetV0tag(Int_t trackIndex) const
{
  //
  // Get the tag for the corresponding trackIndex. Returns -99 in case of invalid index/tag list.
  //


  if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0tags) return -99;
  else
  {
      return fV0tags[trackIndex];
  }
}

//______________________________________________________________________________
void AliTRDPIDTree::ClearV0PIDlist(){

  //
  // Clear the PID object arrays
  //

  fV0electrons->Clear();
  fV0pions->Clear();
  fV0protons->Clear();

  delete[] fV0tags;
  fV0tags = 0;

  fNumTagsStored = 0;
}


//______________________________________________________________________________
void AliTRDPIDTree::SetupV0qa()
{
  //
  // Create the qa objects for V0 Kine cuts
  //
  
  fhArmenteros  = new TH2F("fhArmenteros","Armenteros plot",200,-1.,1.,200,0.,0.4);
  fListQATRDV0->Add(fhArmenteros);
 
}

//________________________________________________________________________
Double_t AliTRDPIDTree::GetPhi(AliESDtrack *const fTrack,Int_t iPl, Double_t &yposlayer)
{
    //
    // extrapolate track to TRD radii and convert global phi angle to local coordinate system
    //

    Double_t phi=-999;
    yposlayer=-999;
    // Phi at entrance of TRD
    Double_t xtrdbeg=AliTRDgeometry::GetXtrdBeg();
    Double_t xtrdend=AliTRDgeometry::GetXtrdEnd();
    Double_t x=xtrdbeg+iPl*(xtrdend-xtrdbeg)/6;
    if(fTrack->GetOuterParam()){
	AliExternalTrackParam param(*fTrack->GetOuterParam());
	param.PropagateTo(x,fESDEvent->GetMagneticField());
	phi=param.Phi()-param.GetAlpha();
        yposlayer= param.GetY();
    }
    if(phi<-TMath::Pi())phi+=2*TMath::Pi();
    if(phi>TMath::Pi())phi-=2*TMath::Pi();

    return phi;
}


//________________________________________________________________________
Bool_t AliTRDPIDTree::PassTrackCuts(AliESDtrack *fESDTrack)
{
    //
    // check if tracks pass minimum quality critieria
    //

    if(!fESDTrack) return kFALSE;

    if(fhtrackCuts) fhtrackCuts->Fill(0);


    // DCA to PV
    Float_t dca[2];
    fESDTrack->GetImpactParameters(dca[0],dca[1]);
    if(dca[0]>5||dca[1]>10){
	if(fhtrackCuts) fhtrackCuts->Fill(1);
	return kFALSE;
    }

    // remove kinks
    if(fESDTrack->GetKinkIndex(0)>0){
      if(fhtrackCuts) fhtrackCuts->Fill(2);
      return kFALSE;
    }

    // eta cut
    if((TMath::Abs(fESDTrack->Eta()))>1.0){
        if(fhtrackCuts) fhtrackCuts->Fill(3);
	return kFALSE;
    }

    // Chi2
    Float_t tpcchi2=99;
    Int_t tpcnclusF=fESDTrack->GetTPCNclsF();
    if(tpcnclusF!=0) tpcchi2=(Float_t)fESDTrack->GetTPCchi2()/tpcnclusF;
    else tpcchi2=1000;
    if(tpcchi2 > 4){
	if(fhtrackCuts) fhtrackCuts->Fill(4);
	return kFALSE;
    }

    //TRD out
    if((fESDTrack->GetStatus()&AliVTrack::kTRDout)==0){
	if(fhtrackCuts) fhtrackCuts->Fill(5);
	return kFALSE;
    }

    // QA #TRD PID tracklets (no cut, just for QA)
    if(fESDTrack->GetTRDntrackletsPID()<4){
	if(fhtrackCuts) fhtrackCuts->Fill(6);
    }

    // QA Missing layers (no cut, just for QA)
    if(HasMissingLayer(fESDTrack)){
	if(fhtrackCuts) fhtrackCuts->Fill(7);
    }

    return kTRUE;
}

//_______________________________________________________________________
Bool_t AliTRDPIDTree::HasMissingLayer(const AliVTrack *fESDTrack){

    //
    // returns if slices missing
    //

    Int_t ntrl=0;
    for(Int_t jPl=0;jPl<AliVTrack::kTRDnPlanes;jPl++){
	Double_t signal=0;
	for(int isl= 0; isl<= 8;isl++){
	    Double_t sigsl=fESDTrack->GetTRDslice(jPl,isl);
	    if(sigsl>0)signal+=sigsl;
	}
        // if signal is missing, stop counting
	if(signal<=0||fESDTrack->GetTRDmomentum(jPl)<=0)break;
	ntrl++;
    }
    //compare with number of layers
    if(ntrl!=GetNTrackletsPID(fESDTrack))return kTRUE;

    return kFALSE;
}

//_______________________________________________________________________
Int_t AliTRDPIDTree::GetNTrackletsPID(const AliVTrack *fESDTrack) const
{
    //
    // returns the number of PID tracklets, re-calculated from dE/dx stored in slices
    //

    Int_t ntracklets=0;

    for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
	Double_t signal=0;
	for(int isl= 0; isl<= 8;isl++){
	    Double_t sigsl=fESDTrack->GetTRDslice(iPl,isl);
	    if(sigsl>0)signal+=sigsl;
	}
	if(signal>0){
	    ntracklets++;
	}
    }
    return ntracklets;
}
Commit	Line	Data
eb10d583	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de> *
	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	// The task:
	17	// stores TRD PID quantities in a Tree
	18	// output can then be used for e.g. parameter creation
	19	//
	20	// Author:
	21	// Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>
	22	//
	23
	24
	25	#include "AliTRDPIDTree.h"
	26	#include "TChain.h"
	27	#include "AliAnalysisManager.h"
	28	#include "AliESDEvent.h"
	29	#include "AliMCEvent.h"
	30	#include "AliPID.h"
	31	#include "AliESDtrack.h"
	32	#include "AliPIDResponse.h"
	33	#include "AliInputEventHandler.h"
	34	#include "AliESDInputHandler.h"
	35	#include "AliESDv0KineCuts.h"
	36	#include "AliESDv0.h"
	37	#include "AliCentrality.h"
	38	#include "AliTRDgeometry.h"
	39	#include "TTree.h"
	40
	41	class TCanvas;
	42	class TAxis;
	43	class TFile;
	44	class TStyle;
	45	class TString;
	46	class TH1F;
	47	class TH2D;
	48	class THnSparse;
	49	class TLegend;
	50	class TVirtualFitter;
	51	class AliESDtrackCuts;
	52	class AliStack;
	53	class AliMCParticle;
	54
	55
	56	using namespace std;
	57
	58	ClassImp(AliTRDPIDTree)
	59
	60	Double_t AliTRDPIDTree::fgMinLayer = 1.;
	61
	62	//________________________________________________________________________
	63	AliTRDPIDTree::AliTRDPIDTree(const char *name)
	64	: AliAnalysisTaskSE(name), fV0tags(0x0), fV0cuts(0x0), fV0electrons(0x0), fV0pions(0x0), fV0protons(0x0),
65	fESDEvent(0), fMCEvent(0), fMCStack(0), fTreeTRDPID(0), fPIDResponse(0), fOutputContainer(0), fESDtrackCuts(0),
66	fESDtrackCutsV0(0), fListQATRD(0x0), fListQATRDV0(0x0),
67	fNumTagsStored(0), fCollisionSystem(3),
68	fpdg(0), frun(0), frunnumber(0), fcentrality(0), fTRDNtracklets(0), fTRDNcls(0), fTRDntracklets(0), fTRDntrackletsPID(0),
370f5346	69	fTRDtheta(0), fTRDsignal(0), fTRDnclsdEdx(0), fTRDnch(0), fPDG(0), fPDGTRUE(0), fChi2(0), fhtrackCuts(0), fhArmenteros(0)
eb10d583	70	{
	71
	72	//
	73	// Constructor
	74	//
	75	DefineInput(0, TChain::Class());
	76	DefineOutput(1, TTree::Class());
	77	DefineOutput(2, TList::Class());
	78
	79	}
	80
	81
	82	//_________________________________________________
	83	AliTRDPIDTree::~AliTRDPIDTree()
	84	{
	85
	86	//
	87	// Destructor
	88	//
	89
	90	delete fV0cuts;
	91	delete fV0electrons;
	92	delete fV0pions;
	93	delete fV0protons;
	94	delete fV0tags;
	95	fV0tags = 0;
	96	fNumTagsStored = 0;
	97	}
	98
	99
	100	//________________________________________________________________________
	101	void AliTRDPIDTree::UserCreateOutputObjects()
	102	{
	103	//
	104	// Definition of user output ntuple and histogram file
	105	//
	106
	107
	108	AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	109	if (!inputHandler)
	110	printf("Inputhandler not available \n");
	111	else
	112	fPIDResponse = inputHandler->GetPIDResponse();
	113
	114	// V0 Kine cuts
	115	fV0cuts = new AliESDv0KineCuts();
	116
	117	// V0 PID Obj arrays
	118	fV0electrons = new TObjArray;
	119	fV0pions = new TObjArray;
	120	fV0protons = new TObjArray;
	121
	122
	123
	124	OpenFile(1);
	125	fTreeTRDPID = new TTree("TreeTRDPID", "TRD PID");
	126	fTreeTRDPID->Branch("run", &frunnumber);
	127	fTreeTRDPID->Branch("centrality", &fcentrality);
	128	fTreeTRDPID->Branch("TRDNtracklets",&fTRDNtracklets);
	129	fTreeTRDPID->Branch("TRDslices[48]",fTRDslices);
	130	fTreeTRDPID->Branch("TRDMomentum[6]",fTRDMomentum);
	131	fTreeTRDPID->Branch("TRDNcls",&fTRDNcls);
	132	fTreeTRDPID->Branch("TRDntracklets",&fTRDntracklets);
	133	fTreeTRDPID->Branch("TRDntrackletsPID",&fTRDntrackletsPID);
134	fTreeTRDPID->Branch("TRDphi[6]",fTRDphi);
135	fTreeTRDPID->Branch("TRDY[6]",fTRDY);
136	fTreeTRDPID->Branch("TRDtheta",&fTRDtheta);
137	fTreeTRDPID->Branch("TRDsignal",&fTRDsignal);
138	fTreeTRDPID->Branch("TRDnclsdEdx",&fTRDnclsdEdx);
139	fTreeTRDPID->Branch("TRDnch",&fTRDnch);
140	fTreeTRDPID->Branch("NSigmaTPC[3]",fNSigmaTPC);
141	fTreeTRDPID->Branch("NSigmaTOF[3]",fNSigmaTOF);
142	fTreeTRDPID->Branch("PDG",&fPDG);
143	fTreeTRDPID->Branch("PDGTRUE",&fPDGTRUE);
144	fTreeTRDPID->Branch("DCA[2]",fDCA);
145	fTreeTRDPID->Branch("Chi2",&fChi2);
370f5346	146	fTreeTRDPID->Branch("TRDsigma[5]",fsigmaTRD);
	147	fTreeTRDPID->Branch("TRDdelta[5]",fdeltaTRD);
	148	fTreeTRDPID->Branch("TRDratio[5]",fratioTRD);
eb10d583	149
	150	PostData(1,fTreeTRDPID);
	151
	152
	153	OpenFile(2);
	154	fListQATRD=new TList;
	155	fListQATRD->SetOwner();
	156	fListQATRDV0=new TList;
	157	fListQATRDV0->SetOwner();
	158	fListQATRDV0->SetName("V0decay");
	159	fListQATRD->Add(fListQATRDV0);
	160
	161	SetupV0qa();
	162
	163	fhtrackCuts = new TH1F("fhtrackCuts","TrackEventCuts QA",10,-0.5,9.5);
	164	fListQATRD->Add(fhtrackCuts);
	165
	166	PostData(2,fListQATRD);
	167
	168
	169	}
	170
	171
	172	//_____________________________________________________________________________
	173	void AliTRDPIDTree::UserExec(Option_t *)
	174	{
	175	//
	176	//calls the Process function
	177	//
	178
	179	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	180
	181	if (!esdH) {
	182	printf("ERROR: Could not get ESDInputHandler \n");
	183	}
	184	else fESDEvent = esdH->GetEvent();
	185
	186	// If MC, forward MCevent
	187	fMCEvent = dynamic_cast<AliMCEvent*>(MCEvent());
	188	fMCStack=NULL;
	189	if(fMCEvent){
	190	fMCStack = fMCEvent->Stack();
	191	}
	192
	193	FillV0PIDlist();
	194	//
	195	Process(fESDEvent, fMCEvent);
	196	//
	197
	198	PostData(1,fTreeTRDPID);
	199	PostData(2,fListQATRD);
	200
	201	// Clear the V0 PID arrays
	202	ClearV0PIDlist();
	203
	204	}
	205
	206
	207
	208	//________________________________________________________________________
	209	void AliTRDPIDTree::Process(AliESDEvent const esdEvent, AliMCEvent const mcEvent)
	210	{
	211	//
	212	//called for each event
213	//
214
215	if (!esdEvent) {
216	Printf("ERROR: esdEvent not available");
217	return;
218	}
219
220	if (!mcEvent) {
221	// Printf("ERROR: mcEvent not available");
222	// return;
223	}
224
225	if (!fPIDResponse ) {
226	Printf("ERROR: No PIDresponse and/or v0KineCuts!");
227	return;
228	}
229
230	Float_t centralityFper=99;
231
232	AliCentrality *esdCentrality = esdEvent->GetCentrality();
233	centralityFper = esdCentrality->GetCentralityPercentile("V0M");
234
235
236	const AliESDVertex* fESDEventvertex = esdEvent->GetPrimaryVertexTracks();
237	if (!fESDEventvertex)
238	return;
239
240	Int_t ncontr = fESDEventvertex->GetNContributors();
241
242	if (ncontr <= 0) return;
243	frunnumber = fESDEvent->GetRunNumber();
244
245	// - Begin: track loop for electrons from V0 -
246	for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
247	// AliVTrack track=(AliVTrack)fV0electrons->At(itrack);
248	AliESDtrack track=(AliESDtrack)fV0electrons->At(itrack);
249
250	fpdg=11;
251	FillTree(track, fpdg, frun, centralityFper);
252
253	} // - End: track loop for electrons from V0 -
254
255	// - Begin: track loop for pions from V0 -
256	for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
257	// AliVTrack track=(AliVTrack)fV0pions->At(itrack);
258	AliESDtrack track=(AliESDtrack)fV0pions->At(itrack);
259	fpdg=211;
260	FillTree(track, fpdg, frun, centralityFper);
261	}
262
263	// - Begin: track loop for protons from V0 -
264	for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
265	// AliVTrack track=(AliVTrack)fV0protons->At(itrack);
266	AliESDtrack track=(AliESDtrack)fV0protons->At(itrack);
267	fpdg=2212;
268	FillTree(track, fpdg, frun, centralityFper);
269	}
270
271
272
273	PostData(1,fTreeTRDPID);
274	PostData(2,fListQATRD);
275
276	}
277
278	//________________________________________________________________________
279	void AliTRDPIDTree::FillTree(AliESDtrack *track, Int_t pdgfromv0, Int_t runnumber, Int_t centralityvalue)
280	{
281	//
282	// Fill tree with track and event variables for e, pi, proton from V0 candidates
283	//
284
285	// global event and track cuts
286	if(!PassTrackCuts(track)) return;
287
288	// track cuts:
289	Int_t ntrackletstracking=track->GetTRDntracklets();
290	if(ntrackletstracking<fgMinLayer) return;
291	Int_t trdnch = track->GetTRDNchamberdEdx(); // only available for p-Pb up to now
292	// if(trdnch<1) return;
293
294
295
296	// check if TOF can provide PID (if not true return)
297	UInt_t status = track->GetStatus();
298	Bool_t hasTOFout = status&AliESDtrack::kTOFout;
299	Bool_t hasTOFtime = status&AliESDtrack::kTIME;
370f5346	300	if(!(hasTOFout && hasTOFtime)) return;
eb10d583	301
eb10d583	302	// add kaons?
41a6c5e0	303	Double_t nSigmaTPC = 0;
41a6c5e0	304	Double_t nSigmaTOF = 0;
eb10d583	305
	306	if(TMath::Abs(pdgfromv0)==11){
	307	nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron);
	308	nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron);
	309	}
	310	if(TMath::Abs(pdgfromv0)==211){
	311	nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
	312	nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
	313	}
	314	if(TMath::Abs(pdgfromv0)==2212){
	315	nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
	316	nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);
	317	}
	318
	319	if(nSigmaTPC>3\|\|nSigmaTPC<-3){
	320	return;
	321	}
	322
	323	if(nSigmaTOF>3\|\|nSigmaTOF<-3){
	324	return;
	325	}
	326
	327	frun=runnumber;
	328
	329	Int_t charge = track->Charge();
	330	fPDG=pdgfromv0*charge;
	331	fcentrality = centralityvalue;
	332	fTRDntracklets = ntrackletstracking;
	333	fTRDntrackletsPID = track->GetTRDntrackletsPID();
	334
	335	fTRDNcls=track->GetTRDncls();
	336	fTRDtheta=track->Theta()-0.5*TMath::Pi();
	337	fChi2=track->GetTRDchi2();
	338
	339	fTRDsignal=track->GetTRDsignal(); // truncated mean signal
	340	fTRDnclsdEdx=track->GetTRDNclusterdEdx(); // number of clusters dedx
	341	fTRDnch=trdnch; // number of chambers dedx
	342
	343	track->GetImpactParameters(fDCA[0],fDCA[1]);
	344
	345	fNSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron);
	346	fNSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
	347	fNSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
	348	fNSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron);
	349	fNSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
	350	fNSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);
	351
	352	Int_t ntracklets=0;
	353	for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
	354
	355	fTRDphi[iPl]=GetPhi(track,iPl,fTRDY[iPl]);
	356	Float_t dEdx=0;
	357	fTRDMomentum[iPl]= track->GetTRDmomentum(iPl);
	358	for(int isl= 0; isl<= 7;isl++){
	359	fTRDslices[iPl*8+isl]=track->GetTRDslice(iPl,isl);
	360	Double_t sig=track->GetTRDslice(iPl,isl);
	361	if(sig>0){
	362	dEdx+=sig;
	363	}
	364	}
	365	if(dEdx>0){
	366	ntracklets++;
	367	}
	368	}
369	fTRDNtracklets=ntracklets;
370
371	fPDGTRUE=0;
372	if(fMCStack){
373	Int_t label=track->GetLabel();
374	if(label>=0){
375	TParticle *part=fMCStack->Particle(label);
376	if(part){
377	fPDGTRUE=part->GetPdgCode();
378	if(TMath::Abs(fPDGTRUE)==11)fPDGTRUE=-fPDGTRUE; // switch signs (particle code says 11 for e- and -11 for e+)
379	}
380	}
381	}
382
eb10d583	383	AliPID::EParticleType types[]={AliPID::kElectron, AliPID::kMuon, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
eb10d583	384	for(Int_t itrdpid=0; itrdpid<5; itrdpid++){
370f5346	385	fsigmaTRD[itrdpid]= fPIDResponse->NumberOfSigmas(AliPIDResponse::kTRD, track, types[itrdpid]); // (signal-expsig)/(sigma), where sigma = mean*resolution
	386	fdeltaTRD[itrdpid]= fPIDResponse->GetSignalDelta(AliPIDResponse::kTRD, track, types[itrdpid]); // signal/(expsig)
	387	fratioTRD[itrdpid]= fPIDResponse->GetSignalDelta(AliPIDResponse::kTRD, track, types[itrdpid], kTRUE); // difference signal to theor. value divided by resolution
eb10d583	388	}
	389
	390
	391	fTreeTRDPID->Fill();
	392
	393	}
	394
	395	//________________________________________________________________________
	396	Int_t AliTRDPIDTree::CompareFloat(Float_t f1, Float_t f2) const
	397	{
	398	//
	399	//compares if the Float_t f1 is equal with f2 and returns 1 if true and 0 if false
	400	//
	401	Float_t precision = 0.00001;
	402	if (((f1 - precision) < f2) &&
	403	((f1 + precision) > f2))
	404	{
	405	return 1;
	406	}
	407	else
	408	{
	409	return 0;
	410	}
	411	}
	412
	413
	414	//________________________________________________________________________
	415	void AliTRDPIDTree::Terminate(const Option_t *)
	416	{
	417	//
	418	// Terminate function
	419	//
	420
	421
	422	}
	423
	424
	425	//______________________________________________________________________________
	426	void AliTRDPIDTree::FillV0PIDlist(){
	427
	428	//
	429	// Fill the PID object arrays holding the pointers to identified particle tracks
	430	//
	431
	432	// Dynamic cast to ESD events (DO NOTHING for AOD events)
	433	AliESDEvent event = dynamic_cast<AliESDEvent >(InputEvent());
	434	if ( !event ) return;
	435
	436	if(IsPbPb()) {
	437	fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPbPb);
	438	}
	439	else {
	440	fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPP);
	441	}
	442
	443
	444	// V0 selection
	445	// set event
	446	fV0cuts->SetEvent(event);
	447
	448	const Int_t numTracks = event->GetNumberOfTracks();
	449	fV0tags = new Int_t[numTracks];
	450	for (Int_t i = 0; i < numTracks; i++)
	451	fV0tags[i] = 0;
452
453	fNumTagsStored = numTracks;
454
455	// loop over V0 particles
456	for(Int_t iv0=0; iv0<event->GetNumberOfV0s();iv0++){
457
458
459	AliESDv0 v0 = (AliESDv0 ) event->GetV0(iv0);
460
461	if(!v0) continue;
462	if(v0->GetOnFlyStatus()) continue;
463	// Get the particle selection
464	Bool_t foundV0 = kFALSE;
465	Int_t pdgV0, pdgP, pdgN;
466	foundV0 = fV0cuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
467	if(!foundV0) continue;
468	Int_t iTrackP = v0->GetPindex(); // positive track
469	Int_t iTrackN = v0->GetNindex(); // negative track
470
471	// v0 Armenteros plot (QA)
472	Float_t armVar[2] = {0.0,0.0};
473	fV0cuts->Armenteros(v0, armVar);
474
475	if(fListQATRDV0&&fhArmenteros)fhArmenteros->Fill(armVar[0],armVar[1]);
476
477	// fill the Object arrays
478	// positive particles
479	if( pdgP == -11){
480	fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackP));
481	fV0tags[iTrackP] = 11;
482	}
483	else if( pdgP == 211){
484	fV0pions->Add((AliVTrack*)event->GetTrack(iTrackP));
485	fV0tags[iTrackP] = 211;
486	}
487	else if( pdgP == 2212){
488	fV0protons->Add((AliVTrack*)event->GetTrack(iTrackP));
489	fV0tags[iTrackP] = 2212;
490	}
491
492	// negative particles
493	if( pdgN == 11){
494	fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackN));
495	fV0tags[iTrackN] = -11;
496	}
497	else if( pdgN == -211){
498	fV0pions->Add((AliVTrack*)event->GetTrack(iTrackN));
499	fV0tags[iTrackN] = -211;
500	}
501	else if( pdgN == -2212){
502	fV0protons->Add((AliVTrack*)event->GetTrack(iTrackN));
503	fV0tags[iTrackN] = -2212;
504	}
505
506
507	}
508	}
509
510	//______________________________________________________________________________
511	Int_t AliTRDPIDTree::GetV0tag(Int_t trackIndex) const
512	{
513	//
514	// Get the tag for the corresponding trackIndex. Returns -99 in case of invalid index/tag list.
515	//
516
517
518
519	if (trackIndex < 0 \|\| trackIndex >= fNumTagsStored \|\| !fV0tags) return -99;
520	else
521	{
522	return fV0tags[trackIndex];
523	}
524	}
525
526	//______________________________________________________________________________
527	void AliTRDPIDTree::ClearV0PIDlist(){
528
529	//
530	// Clear the PID object arrays
531	//
532
533	fV0electrons->Clear();
534	fV0pions->Clear();
535	fV0protons->Clear();
536
41a6c5e0	537	delete[] fV0tags;
eb10d583	538	fV0tags = 0;
	539
	540	fNumTagsStored = 0;
	541	}
	542
	543
	544	//______________________________________________________________________________
	545	void AliTRDPIDTree::SetupV0qa()
	546	{
	547	//
	548	// Create the qa objects for V0 Kine cuts
	549	//
	550
	551	fhArmenteros = new TH2F("fhArmenteros","Armenteros plot",200,-1.,1.,200,0.,0.4);
	552	fListQATRDV0->Add(fhArmenteros);
	553
	554	}
	555
	556	//________________________________________________________________________
	557	Double_t AliTRDPIDTree::GetPhi(AliESDtrack *const fTrack,Int_t iPl, Double_t &yposlayer)
	558	{
	559	//
	560	// extrapolate track to TRD radii and convert global phi angle to local coordinate system
	561	//
	562
	563	Double_t phi=-999;
	564	yposlayer=-999;
	565	// Phi at entrance of TRD
	566	Double_t xtrdbeg=AliTRDgeometry::GetXtrdBeg();
	567	Double_t xtrdend=AliTRDgeometry::GetXtrdEnd();
	568	Double_t x=xtrdbeg+iPl*(xtrdend-xtrdbeg)/6;
	569	if(fTrack->GetOuterParam()){
	570	AliExternalTrackParam param(*fTrack->GetOuterParam());
	571	param.PropagateTo(x,fESDEvent->GetMagneticField());
	572	phi=param.Phi()-param.GetAlpha();
	573	yposlayer= param.GetY();
	574	}
	575	if(phi<-TMath::Pi())phi+=2*TMath::Pi();
	576	if(phi>TMath::Pi())phi-=2*TMath::Pi();
	577
	578	return phi;
	579	}
	580
	581
	582	//________________________________________________________________________
	583	Bool_t AliTRDPIDTree::PassTrackCuts(AliESDtrack *fESDTrack)
	584	{
	585	//
	586	// check if tracks pass minimum quality critieria
	587	//
	588
	589	if(!fESDTrack) return kFALSE;
	590
	591	if(fhtrackCuts) fhtrackCuts->Fill(0);
	592
	593
	594	// DCA to PV
	595	Float_t dca[2];
	596	fESDTrack->GetImpactParameters(dca[0],dca[1]);
	597	if(dca[0]>5\|\|dca[1]>10){
	598	if(fhtrackCuts) fhtrackCuts->Fill(1);
	599	return kFALSE;
	600	}
	601
602	// remove kinks
603	if(fESDTrack->GetKinkIndex(0)>0){
604	if(fhtrackCuts) fhtrackCuts->Fill(2);
605	return kFALSE;
606	}
607
608	// eta cut
609	if((TMath::Abs(fESDTrack->Eta()))>1.0){
610	if(fhtrackCuts) fhtrackCuts->Fill(3);
611	return kFALSE;
612	}
613
614	// Chi2
615	Float_t tpcchi2=99;
616	Int_t tpcnclusF=fESDTrack->GetTPCNclsF();
617	if(tpcnclusF!=0) tpcchi2=(Float_t)fESDTrack->GetTPCchi2()/tpcnclusF;
618	else tpcchi2=1000;
619	if(tpcchi2 > 4){
620	if(fhtrackCuts) fhtrackCuts->Fill(4);
621	return kFALSE;
622	}
623
624	//TRD out
625	if((fESDTrack->GetStatus()&AliVTrack::kTRDout)==0){
626	if(fhtrackCuts) fhtrackCuts->Fill(5);
627	return kFALSE;
628	}
629
630	// QA #TRD PID tracklets (no cut, just for QA)
631	if(fESDTrack->GetTRDntrackletsPID()<4){
632	if(fhtrackCuts) fhtrackCuts->Fill(6);
633	}
634
635	// QA Missing layers (no cut, just for QA)
636	if(HasMissingLayer(fESDTrack)){
637	if(fhtrackCuts) fhtrackCuts->Fill(7);
638	}
639
640	return kTRUE;
641	}
642
643	//_______________________________________________________________________
644	Bool_t AliTRDPIDTree::HasMissingLayer(const AliVTrack *fESDTrack){
645
646	//
647	// returns if slices missing
648	//
649
650	Int_t ntrl=0;
651	for(Int_t jPl=0;jPl<AliVTrack::kTRDnPlanes;jPl++){
652	Double_t signal=0;
653	for(int isl= 0; isl<= 8;isl++){
654	Double_t sigsl=fESDTrack->GetTRDslice(jPl,isl);
655	if(sigsl>0)signal+=sigsl;
656	}
657	// if signal is missing, stop counting
658	if(signal<=0\|\|fESDTrack->GetTRDmomentum(jPl)<=0)break;
659	ntrl++;
660	}
661	//compare with number of layers
662	if(ntrl!=GetNTrackletsPID(fESDTrack))return kTRUE;
663
664	return kFALSE;
665	}
666
667	//_______________________________________________________________________
668	Int_t AliTRDPIDTree::GetNTrackletsPID(const AliVTrack *fESDTrack) const
669	{
670	//
671	// returns the number of PID tracklets, re-calculated from dE/dx stored in slices
672	//
673
674	Int_t ntracklets=0;
675
676	for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
677	Double_t signal=0;
678	for(int isl= 0; isl<= 8;isl++){
679	Double_t sigsl=fESDTrack->GetTRDslice(iPl,isl);
680	if(sigsl>0)signal+=sigsl;
681	}
682	if(signal>0){
683	ntracklets++;
684	}
685	}
686	return ntracklets;
687	}