[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, nSigmaTOF;

    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
	302	// add kaons?
	303	Double_t nSigmaTPC, nSigmaTOF;
	304
	305	if(TMath::Abs(pdgfromv0)==11){
	306	nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron);
	307	nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron);
	308	}
	309	if(TMath::Abs(pdgfromv0)==211){
	310	nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
	311	nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
	312	}
	313	if(TMath::Abs(pdgfromv0)==2212){
	314	nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
	315	nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);
	316	}
	317
	318	if(nSigmaTPC>3\|\|nSigmaTPC<-3){
	319	return;
	320	}
	321
	322	if(nSigmaTOF>3\|\|nSigmaTOF<-3){
	323	return;
	324	}
	325
	326	frun=runnumber;
	327
	328	Int_t charge = track->Charge();
	329	fPDG=pdgfromv0*charge;
	330	fcentrality = centralityvalue;
	331	fTRDntracklets = ntrackletstracking;
	332	fTRDntrackletsPID = track->GetTRDntrackletsPID();
	333
	334	fTRDNcls=track->GetTRDncls();
	335	fTRDtheta=track->Theta()-0.5*TMath::Pi();
	336	fChi2=track->GetTRDchi2();
	337
	338	fTRDsignal=track->GetTRDsignal(); // truncated mean signal
	339	fTRDnclsdEdx=track->GetTRDNclusterdEdx(); // number of clusters dedx
	340	fTRDnch=trdnch; // number of chambers dedx
	341
	342	track->GetImpactParameters(fDCA[0],fDCA[1]);
	343
	344	fNSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron);
	345	fNSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
	346	fNSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
	347	fNSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron);
	348	fNSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
	349	fNSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);
	350
	351	Int_t ntracklets=0;
	352	for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
	353
	354	fTRDphi[iPl]=GetPhi(track,iPl,fTRDY[iPl]);
	355	Float_t dEdx=0;
	356	fTRDMomentum[iPl]= track->GetTRDmomentum(iPl);
	357	for(int isl= 0; isl<= 7;isl++){
	358	fTRDslices[iPl*8+isl]=track->GetTRDslice(iPl,isl);
	359	Double_t sig=track->GetTRDslice(iPl,isl);
	360	if(sig>0){
	361	dEdx+=sig;
	362	}
	363	}
	364	if(dEdx>0){
365	ntracklets++;
366	}
367	}
368	fTRDNtracklets=ntracklets;
369
370	fPDGTRUE=0;
371	if(fMCStack){
372	Int_t label=track->GetLabel();
373	if(label>=0){
374	TParticle *part=fMCStack->Particle(label);
375	if(part){
376	fPDGTRUE=part->GetPdgCode();
377	if(TMath::Abs(fPDGTRUE)==11)fPDGTRUE=-fPDGTRUE; // switch signs (particle code says 11 for e- and -11 for e+)
378	}
379	}
380	}
381
eb10d583	382	AliPID::EParticleType types[]={AliPID::kElectron, AliPID::kMuon, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
eb10d583	383	for(Int_t itrdpid=0; itrdpid<5; itrdpid++){
370f5346	384	fsigmaTRD[itrdpid]= fPIDResponse->NumberOfSigmas(AliPIDResponse::kTRD, track, types[itrdpid]); // (signal-expsig)/(sigma), where sigma = mean*resolution
	385	fdeltaTRD[itrdpid]= fPIDResponse->GetSignalDelta(AliPIDResponse::kTRD, track, types[itrdpid]); // signal/(expsig)
	386	fratioTRD[itrdpid]= fPIDResponse->GetSignalDelta(AliPIDResponse::kTRD, track, types[itrdpid], kTRUE); // difference signal to theor. value divided by resolution
eb10d583	387	}
	388
	389
	390	fTreeTRDPID->Fill();
	391
	392	}
	393
	394	//________________________________________________________________________
	395	Int_t AliTRDPIDTree::CompareFloat(Float_t f1, Float_t f2) const
	396	{
	397	//
	398	//compares if the Float_t f1 is equal with f2 and returns 1 if true and 0 if false
	399	//
	400	Float_t precision = 0.00001;
	401	if (((f1 - precision) < f2) &&
	402	((f1 + precision) > f2))
	403	{
	404	return 1;
	405	}
	406	else
	407	{
	408	return 0;
	409	}
	410	}
	411
	412
	413	//________________________________________________________________________
	414	void AliTRDPIDTree::Terminate(const Option_t *)
	415	{
	416	//
	417	// Terminate function
	418	//
	419
	420
	421	}
	422
	423
	424	//______________________________________________________________________________
	425	void AliTRDPIDTree::FillV0PIDlist(){
	426
	427	//
	428	// Fill the PID object arrays holding the pointers to identified particle tracks
	429	//
	430
	431	// Dynamic cast to ESD events (DO NOTHING for AOD events)
	432	AliESDEvent event = dynamic_cast<AliESDEvent >(InputEvent());
	433	if ( !event ) return;
	434
	435	if(IsPbPb()) {
	436	fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPbPb);
	437	}
	438	else {
	439	fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPP);
	440	}
	441
	442
	443	// V0 selection
	444	// set event
	445	fV0cuts->SetEvent(event);
	446
	447	const Int_t numTracks = event->GetNumberOfTracks();
	448	fV0tags = new Int_t[numTracks];
	449	for (Int_t i = 0; i < numTracks; i++)
	450	fV0tags[i] = 0;
451
452	fNumTagsStored = numTracks;
453
454	// loop over V0 particles
455	for(Int_t iv0=0; iv0<event->GetNumberOfV0s();iv0++){
456
457
458	AliESDv0 v0 = (AliESDv0 ) event->GetV0(iv0);
459
460	if(!v0) continue;
461	if(v0->GetOnFlyStatus()) continue;
462	// Get the particle selection
463	Bool_t foundV0 = kFALSE;
464	Int_t pdgV0, pdgP, pdgN;
465	foundV0 = fV0cuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
466	if(!foundV0) continue;
467	Int_t iTrackP = v0->GetPindex(); // positive track
468	Int_t iTrackN = v0->GetNindex(); // negative track
469
470	// v0 Armenteros plot (QA)
471	Float_t armVar[2] = {0.0,0.0};
472	fV0cuts->Armenteros(v0, armVar);
473
474	if(fListQATRDV0&&fhArmenteros)fhArmenteros->Fill(armVar[0],armVar[1]);
475
476	// fill the Object arrays
477	// positive particles
478	if( pdgP == -11){
479	fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackP));
480	fV0tags[iTrackP] = 11;
481	}
482	else if( pdgP == 211){
483	fV0pions->Add((AliVTrack*)event->GetTrack(iTrackP));
484	fV0tags[iTrackP] = 211;
485	}
486	else if( pdgP == 2212){
487	fV0protons->Add((AliVTrack*)event->GetTrack(iTrackP));
488	fV0tags[iTrackP] = 2212;
489	}
490
491	// negative particles
492	if( pdgN == 11){
493	fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackN));
494	fV0tags[iTrackN] = -11;
495	}
496	else if( pdgN == -211){
497	fV0pions->Add((AliVTrack*)event->GetTrack(iTrackN));
498	fV0tags[iTrackN] = -211;
499	}
500	else if( pdgN == -2212){
501	fV0protons->Add((AliVTrack*)event->GetTrack(iTrackN));
502	fV0tags[iTrackN] = -2212;
503	}
504
505
506	}
507	}
508
509	//______________________________________________________________________________
510	Int_t AliTRDPIDTree::GetV0tag(Int_t trackIndex) const
511	{
512	//
513	// Get the tag for the corresponding trackIndex. Returns -99 in case of invalid index/tag list.
514	//
515
516
517
518	if (trackIndex < 0 \|\| trackIndex >= fNumTagsStored \|\| !fV0tags) return -99;
519	else
520	{
521	return fV0tags[trackIndex];
522	}
523	}
524
525	//______________________________________________________________________________
526	void AliTRDPIDTree::ClearV0PIDlist(){
527
528	//
529	// Clear the PID object arrays
530	//
531
532	fV0electrons->Clear();
533	fV0pions->Clear();
534	fV0protons->Clear();
535
536	delete fV0tags;
537	fV0tags = 0;
538
539	fNumTagsStored = 0;
540	}
541
542
543	//______________________________________________________________________________
544	void AliTRDPIDTree::SetupV0qa()
545	{
546	//
547	// Create the qa objects for V0 Kine cuts
548	//
549
550	fhArmenteros = new TH2F("fhArmenteros","Armenteros plot",200,-1.,1.,200,0.,0.4);
551	fListQATRDV0->Add(fhArmenteros);
552
553	}
554
555	//________________________________________________________________________
556	Double_t AliTRDPIDTree::GetPhi(AliESDtrack *const fTrack,Int_t iPl, Double_t &yposlayer)
557	{
558	//
559	// extrapolate track to TRD radii and convert global phi angle to local coordinate system
560	//
561
562	Double_t phi=-999;
563	yposlayer=-999;
564	// Phi at entrance of TRD
565	Double_t xtrdbeg=AliTRDgeometry::GetXtrdBeg();
566	Double_t xtrdend=AliTRDgeometry::GetXtrdEnd();
567	Double_t x=xtrdbeg+iPl*(xtrdend-xtrdbeg)/6;
568	if(fTrack->GetOuterParam()){
569	AliExternalTrackParam param(*fTrack->GetOuterParam());
570	param.PropagateTo(x,fESDEvent->GetMagneticField());
571	phi=param.Phi()-param.GetAlpha();
572	yposlayer= param.GetY();
573	}
574	if(phi<-TMath::Pi())phi+=2*TMath::Pi();
575	if(phi>TMath::Pi())phi-=2*TMath::Pi();
576
577	return phi;
578	}
579
580
581	//________________________________________________________________________
582	Bool_t AliTRDPIDTree::PassTrackCuts(AliESDtrack *fESDTrack)
583	{
584	//
585	// check if tracks pass minimum quality critieria
586	//
587
588	if(!fESDTrack) return kFALSE;
589
590	if(fhtrackCuts) fhtrackCuts->Fill(0);
591
592
593	// DCA to PV
594	Float_t dca[2];
595	fESDTrack->GetImpactParameters(dca[0],dca[1]);
596	if(dca[0]>5\|\|dca[1]>10){
597	if(fhtrackCuts) fhtrackCuts->Fill(1);
598	return kFALSE;
599	}
600
601	// remove kinks
602	if(fESDTrack->GetKinkIndex(0)>0){
603	if(fhtrackCuts) fhtrackCuts->Fill(2);
604	return kFALSE;
605	}
606
607	// eta cut
608	if((TMath::Abs(fESDTrack->Eta()))>1.0){
609	if(fhtrackCuts) fhtrackCuts->Fill(3);
610	return kFALSE;
611	}
612
613	// Chi2
614	Float_t tpcchi2=99;
615	Int_t tpcnclusF=fESDTrack->GetTPCNclsF();
616	if(tpcnclusF!=0) tpcchi2=(Float_t)fESDTrack->GetTPCchi2()/tpcnclusF;
617	else tpcchi2=1000;
618	if(tpcchi2 > 4){
619	if(fhtrackCuts) fhtrackCuts->Fill(4);
620	return kFALSE;
621	}
622
623	//TRD out
624	if((fESDTrack->GetStatus()&AliVTrack::kTRDout)==0){
625	if(fhtrackCuts) fhtrackCuts->Fill(5);
626	return kFALSE;
627	}
628
629	// QA #TRD PID tracklets (no cut, just for QA)
630	if(fESDTrack->GetTRDntrackletsPID()<4){
631	if(fhtrackCuts) fhtrackCuts->Fill(6);
632	}
633
634	// QA Missing layers (no cut, just for QA)
635	if(HasMissingLayer(fESDTrack)){
636	if(fhtrackCuts) fhtrackCuts->Fill(7);
637	}
638
639	return kTRUE;
640	}
641
642	//_______________________________________________________________________
643	Bool_t AliTRDPIDTree::HasMissingLayer(const AliVTrack *fESDTrack){
644
645	//
646	// returns if slices missing
647	//
648
649	Int_t ntrl=0;
650	for(Int_t jPl=0;jPl<AliVTrack::kTRDnPlanes;jPl++){
651	Double_t signal=0;
652	for(int isl= 0; isl<= 8;isl++){
653	Double_t sigsl=fESDTrack->GetTRDslice(jPl,isl);
654	if(sigsl>0)signal+=sigsl;
655	}
656	// if signal is missing, stop counting
657	if(signal<=0\|\|fESDTrack->GetTRDmomentum(jPl)<=0)break;
658	ntrl++;
659	}
660	//compare with number of layers
661	if(ntrl!=GetNTrackletsPID(fESDTrack))return kTRUE;
662
663	return kFALSE;
664	}
665
666	//_______________________________________________________________________
667	Int_t AliTRDPIDTree::GetNTrackletsPID(const AliVTrack *fESDTrack) const
668	{
669	//
670	// returns the number of PID tracklets, re-calculated from dE/dx stored in slices
671	//
672
673	Int_t ntracklets=0;
674
675	for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
676	Double_t signal=0;
677	for(int isl= 0; isl<= 8;isl++){
678	Double_t sigsl=fESDTrack->GetTRDslice(iPl,isl);
679	if(sigsl>0)signal+=sigsl;
680	}
681	if(signal>0){
682	ntracklets++;
683	}
684	}
685	return ntracklets;
686	}