[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEdebugTreeTask.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.                  *
**************************************************************************/
//
// Debug tree task
// 
// Authors:
//   Markus Fasel <M.Fasel@gsi.de>
//
#include <TBits.h>
#include <TString.h>

#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "AliESDVertex.h"
#include "AliHFEcuts.h"
#include "AliHFEsignalCuts.h"
#include "AliInputEventHandler.h"
#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliMCParticle.h"
#include "AliPIDResponse.h"
#include "AliVEvent.h"
#include "TTreeStream.h"

#include "AliHFEdebugTreeTask.h"

ClassImp(AliHFEdebugTreeTask)

AliHFEdebugTreeTask::AliHFEdebugTreeTask():
  AliAnalysisTaskSE(),
  fTrackCuts(NULL),
  fSignalCuts(NULL),
  fNclustersTPC(70),
  fNclustersTPCPID(0),
  fNclustersITS(2),
  fFilename("HFEtree.root"),
  fDebugTree(NULL)
{

}

AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name):
  AliAnalysisTaskSE(name),
  fTrackCuts(NULL),
  fSignalCuts(NULL),
  fNclustersTPC(70),
  fNclustersTPCPID(0),
  fNclustersITS(2),
  fFilename("HFEtree.root"),
  fDebugTree(NULL)
{

}

AliHFEdebugTreeTask::~AliHFEdebugTreeTask(){
  if(fDebugTree) delete fDebugTree;
}

void AliHFEdebugTreeTask::UserCreateOutputObjects(){
  //
  // Create debug tree, signal cuts and track cuts
  //
  fDebugTree = new TTreeSRedirector(fFilename.Data());

  fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");

  fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
  fTrackCuts->CreateStandardCuts();
  // Track cuts
  fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
  fTrackCuts->SetMinRatioTPCclusters(0);
  fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable); 
  fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
  fTrackCuts->SetMinNClustersITS(fNclustersITS);
  // Event cuts
  fTrackCuts->SetUseMixedVertex(kTRUE);
  fTrackCuts->SetVertexRange(10.);
  fTrackCuts->Initialize();
}

void AliHFEdebugTreeTask::UserExec(Option_t *){
  //
  // User Exec: Fill debug Tree
  // 

  // Get PID response
  AliPIDResponse *pid = NULL;
  AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  if(handler){
    pid = handler->GetPIDResponse();
  } else {
    AliError("No Handler");
  }
  if(!pid){
    AliError("No PID response\n");
    return;
  }
  
  fTrackCuts->SetRecEvent(fInputEvent);

  // Cut event
  if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
   AliDebug(1, "Event flagged as pileup\n");
    return;
  }
  if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
    AliDebug(1, "Event rejected by the event cuts\n");
    return;
  }

  // Get run number
  Int_t run = fInputEvent->GetRunNumber();

  // Derive trigger 
  UInt_t trigger = fInputHandler->IsEventSelected();
  Bool_t isMBTrigger = trigger & AliVEvent::kMB;
  Bool_t isCentralTrigger = trigger & AliVEvent::kCentral;
  Bool_t isSemicentralTrigger = trigger & AliVEvent::kSemiCentral;
  Bool_t isEMCALTrigger = trigger & AliVEvent::kEMCEJE;

  // Check if MC information is available
  Bool_t mcthere = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) != NULL;
  if(mcthere){ 
    fTrackCuts->SetMCEvent(fMCEvent);
    fSignalCuts->SetMCEvent(fMCEvent);
  }

  // Get Primary Vertex
  const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();

  // Get centrality
  Float_t centrality = -1.;
  TString beamtype = (dynamic_cast<AliESDEvent *>(fInputEvent))->GetESDRun()->GetBeamType();
  if(!beamtype.CompareTo("Pb-Pb") || !beamtype.CompareTo("A-A")){
    // Heavy ion run
    AliDebug(1, "Heavy-Ion event\n");
    AliCentrality *hicent = fInputEvent->GetCentrality();
    centrality = hicent->GetCentralityPercentile("V0M");
  }

  // Common variables
  Double_t charge, eta, phi, momentum, transversemomentum;
  Int_t source;
      
  // Monte-Carlo loop
  if(mcthere){
    AliMCParticle * mcpart;
    for(Int_t itrk = 0; itrk < fMCEvent->GetNumberOfTracks(); itrk++){
      mcpart = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrk));
      if(!mcpart) continue;
      if(!fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mcpart)) continue;
      // Kinematics
      charge = mcpart->Charge() > 0. ? 1. : -1.;
      momentum = mcpart->P() * charge;
      transversemomentum = mcpart->Pt() * charge;
      eta = mcpart->Eta();
      phi = mcpart->Phi();
      Int_t pdg = mcpart->Particle()->GetPdgCode();
      
      // Get Production Vertex in radial direction
      Double_t vx = mcpart->Particle()->Vx(), 
               vy = mcpart->Particle()->Vy(); 
      Double_t productionVertex = TMath::Sqrt(vx*vx+vy*vy);
      
      // Get Mother PDG code of the particle
      Int_t motherPdg = 0;
      Int_t motherlabel = mcpart->Particle()->GetFirstMother();
      if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){
        AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherlabel));
        if(mother) motherPdg = mother->Particle()->GetPdgCode();
      }
      
      // derive source
      source = 5;
      if(fSignalCuts->IsCharmElectron(mcpart)) source = 0;
      else if(fSignalCuts->IsBeautyElectron(mcpart)) source = 1;
      else if(fSignalCuts->IsGammaElectron(mcpart)) source = 2;
      else if(fSignalCuts->IsNonHFElectron(mcpart)) source = 3;
      else if(TMath::Abs(pdg) == 11) source = 4;
      else source = 5;
      
      (*fDebugTree) << "MCDebug" 
                    << "centrality="          << centrality
                    << "MBtrigger="           << isMBTrigger 
                    << "CentralTrigger="      << isCentralTrigger
                    << "SemicentralTrigger="  << isSemicentralTrigger
                    << "EMCALtrigger="        << isEMCALTrigger
                    << "run="                 << run
                    << "p="                   << momentum
                    << "pt="                  << transversemomentum
                    << "eta="                 << eta
                    << "phi="                 << phi
                    << "pdg="                 << pdg
                    << "ProductionVertex="    << productionVertex
                    << "motherPdg="           << motherPdg
                    << "source="              << source
                    << "\n";
    }
  }
  
  AliESDtrack *track;
  for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
    // fill the tree
    track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrack));
    if(!track) continue;
    // Cut track (Only basic track cuts)
    if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
    // Debug streaming of PID-related quantities
    Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
    Double_t nSigmaTPC = pid->NumberOfSigmasTPC(track, AliPID::kElectron);
    if(TMath::Abs(nSigmaTOF) > 5) continue;
    // we are not interested in tracks which are more than 5 sigma away from the electron hypothesis in either TOF or TPC
    Double_t tPCdEdx = track->GetTPCsignal();
    // Signal, source and MCPID
    Bool_t signal = kTRUE;
    source = 5;
    if(mcthere){      
      // Signal
      AliMCParticle *mctrack;
      if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))){
        if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE; 
      }
  
      // Source
      if(fSignalCuts->IsCharmElectron(track)) source = 0;
      else if(fSignalCuts->IsBeautyElectron(track)) source = 1;
      else if(fSignalCuts->IsGammaElectron(track)) source = 2;
      else if(fSignalCuts->IsNonHFElectron(track)) source = 3;
      else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)) source = 4;
      else source = 5;
    }
    // Get V0 tag (if available)
    Int_t v0pid = -1;
    if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron;
    else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion;
    else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton;
    // Kinematics
    charge = track->Charge() > 0 ? 1. : -1.;
    eta = track->Eta();
    phi = track->Phi();
    momentum = track->P() * charge;
    transversemomentum = track->Pt() * charge;
    // ITS number of clusters
    UChar_t nclustersITS = track->GetITSclusters(NULL);
    Double_t chi2matching =  track->GetChi2TPCConstrainedVsGlobal(dynamic_cast<const AliESDVertex *>(vertex));
    Double_t chi2PerClusterITS = 0.0;
    if (nclustersITS != 0) chi2PerClusterITS = track->GetITSchi2() / Float_t(nclustersITS);
    // TPC number of clusters (different definitions)
    UChar_t nclustersTPC = track->GetTPCncls();
    UChar_t nclustersTPCPID = track->GetTPCsignalN();
    UChar_t nfindableTPC =  track->GetTPCNclsF();
    Double_t clusterRatioTPC = 0.0;
    if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPC = static_cast<Double_t>(nclustersTPC)/static_cast<Double_t>(nfindableTPC);
    UChar_t nclustersTPCshared = 0;
    Float_t ncrossedRowsTPC = track->GetTPCCrossedRows();
    const TBits &sharedTPC = track->GetTPCSharedMap();
    for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
    // TRD clusters and tracklets
    UChar_t nclustersTRD = track->GetTRDncls();
    UChar_t ntracklets = track->GetTRDntrackletsPID();
    // ITS and TRD acceptance maps
    UChar_t hasClusterITS[6], hasTrackletTRD[6];
    UChar_t itsPixel = track->GetITSClusterMap();
    for(Int_t icl = 0; icl < 6; icl++) hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
    for(Int_t itl = 0; itl < 6; itl++){
      Int_t nSliceNonZero = 0;
      for(Int_t islice = 0; islice < 8; islice++){
        if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++;
      }
      hasTrackletTRD[itl] = nSliceNonZero ? 1 : 0;
    }
    // TRD PID
    Double_t pidprobs[5];
    track->GetTRDpid(pidprobs);
    Double_t likeEleTRD = pidprobs[0];
    Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]);
    // DCA
    Float_t b[2] = {0.,0.};
    Float_t bCov[3] = {0.,0.,0.};
    track->GetImpactParameters(b,bCov);
    Double_t dca = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); // impact parameter space
    Double_t dcaSR=0, dcaSZ=0;
    if(bCov[0]>0) dcaSR = b[0]/TMath::Sqrt(bCov[0]); // normalised impact parameter xy
    if(bCov[2]>0) dcaSZ = b[1]/TMath::Sqrt(bCov[2]); // normalised impact parameter z
    Double_t dcaS = AliESDtrackCuts::GetSigmaToVertex(track); // n_sigma
    // Vertex
    Double_t vtx[3];
    vertex->GetXYZ(vtx);
    Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
    // Fill Tree
    (*fDebugTree) << "PIDdebug"
                  << "centrality="          << centrality
                  << "MBtrigger="           << isMBTrigger 
                  << "CentralTrigger="      << isCentralTrigger
                  << "SemicentralTrigger="  << isSemicentralTrigger
                  << "EMCALtrigger="        << isEMCALTrigger
                  << "signal="              << signal
                  << "source="              << source
                  << "v0pid="               << v0pid
                  << "run="                 << run
                  << "p="                   << momentum
                  << "pt="                  << transversemomentum
                  << "eta="                 << eta
                  << "phi="                 << phi
                  << "ntracklets="          << ntracklets
                  << "nclustersTPC="        << nclustersTPC
                  << "nclustersTPCPID="     << nclustersTPCPID
                  << "nclustersTPCshared="  << nclustersTPCshared
                  << "ncrossedRowsTPC="     << ncrossedRowsTPC
                  << "clusterRatioTPC="     << clusterRatioTPC
                  << "nclustersITS="        << nclustersITS
                  << "nclusters="           << nclustersTRD
                  << "chi2matching="        << chi2matching
		  << "chi2PerClusterITS="   << chi2PerClusterITS
                  << "its0="                << hasClusterITS[0]
                  << "its1="                << hasClusterITS[1]
                  << "its2="                << hasClusterITS[2]
                  << "its3="                << hasClusterITS[3]
                  << "its4="                << hasClusterITS[4]
                  << "its5="                << hasClusterITS[5]
                  << "trd0="                << hasTrackletTRD[0]
                  << "trd1="                << hasTrackletTRD[1]
                  << "trd2="                << hasTrackletTRD[2]
                  << "trd3="                << hasTrackletTRD[3]
                  << "trd4="                << hasTrackletTRD[4]
                  << "trd5="                << hasTrackletTRD[5]
                  << "TOFsigmaEl="          << nSigmaTOF
                  << "TPCsigmaEl="          << nSigmaTPC
                  << "TPCdEdx="             << tPCdEdx
                  << "TRDlikeEl="           << likeEleTRD
                  << "TRDlikeEln="          << likeEleTRDn
                  << "dcaR="                << b[0]
                  << "dcaZ="                << b[1]
                  << "dca="                 << dca
                  << "dcaSR="               << dcaSR
                  << "dcaSZ="               << dcaSZ
                  << "dcaS="                << dcaS
                  << "vx="                  << vtx[0]
                  << "vy="                  << vtx[1]
                  << "vz="                  << vtx[2] 
                  << "ncontrib="            << ncontrib
                  << "\n";
  }
}


void AliHFEdebugTreeTask::SetFileName(const char *filename){ fFilename = filename; }
Commit	Line	Data
fda75e2d	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	// Debug tree task
	17	//
	18	// Authors:
	19	// Markus Fasel <M.Fasel@gsi.de>
	20	//
	21	#include <TBits.h>
	22	#include <TString.h>
	23
	24	#include "AliAnalysisManager.h"
	25	#include "AliCentrality.h"
	26	#include "AliESDEvent.h"
	27	#include "AliESDtrack.h"
	28	#include "AliESDtrackCuts.h"
	29	#include "AliESDVertex.h"
	30	#include "AliHFEcuts.h"
	31	#include "AliHFEsignalCuts.h"
	32	#include "AliInputEventHandler.h"
	33	#include "AliLog.h"
	34	#include "AliMCEvent.h"
	35	#include "AliMCEventHandler.h"
	36	#include "AliMCParticle.h"
	37	#include "AliPIDResponse.h"
	38	#include "AliVEvent.h"
	39	#include "TTreeStream.h"
	40
	41	#include "AliHFEdebugTreeTask.h"
	42
	43	ClassImp(AliHFEdebugTreeTask)
	44
	45	AliHFEdebugTreeTask::AliHFEdebugTreeTask():
	46	AliAnalysisTaskSE(),
	47	fTrackCuts(NULL),
	48	fSignalCuts(NULL),
	49	fNclustersTPC(70),
	50	fNclustersTPCPID(0),
	51	fNclustersITS(2),
	52	fFilename("HFEtree.root"),
	53	fDebugTree(NULL)
	54	{
	55
	56	}
	57
	58	AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name):
	59	AliAnalysisTaskSE(name),
	60	fTrackCuts(NULL),
	61	fSignalCuts(NULL),
	62	fNclustersTPC(70),
	63	fNclustersTPCPID(0),
	64	fNclustersITS(2),
65	fFilename("HFEtree.root"),
66	fDebugTree(NULL)
67	{
68
69	}
70
71	AliHFEdebugTreeTask::~AliHFEdebugTreeTask(){
72	if(fDebugTree) delete fDebugTree;
73	}
74
75	void AliHFEdebugTreeTask::UserCreateOutputObjects(){
76	//
77	// Create debug tree, signal cuts and track cuts
78	//
79	fDebugTree = new TTreeSRedirector(fFilename.Data());
80
81	fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
82
83	fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
84	fTrackCuts->CreateStandardCuts();
85	// Track cuts
86	fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
87	fTrackCuts->SetMinRatioTPCclusters(0);
88	fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
89	fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
90	fTrackCuts->SetMinNClustersITS(fNclustersITS);
91	// Event cuts
92	fTrackCuts->SetUseMixedVertex(kTRUE);
93	fTrackCuts->SetVertexRange(10.);
94	fTrackCuts->Initialize();
95	}
96
97	void AliHFEdebugTreeTask::UserExec(Option_t *){
98	//
99	// User Exec: Fill debug Tree
100	//
101
102	// Get PID response
103	AliPIDResponse *pid = NULL;
104	AliInputEventHandler handler = dynamic_cast<AliInputEventHandler >(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
105	if(handler){
106	pid = handler->GetPIDResponse();
107	} else {
108	AliError("No Handler");
109	}
110	if(!pid){
111	AliError("No PID response\n");
112	return;
113	}
114
115	fTrackCuts->SetRecEvent(fInputEvent);
116
117	// Cut event
118	if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
119	AliDebug(1, "Event flagged as pileup\n");
120	return;
121	}
122	if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
123	AliDebug(1, "Event rejected by the event cuts\n");
124	return;
125	}
126
127	// Get run number
128	Int_t run = fInputEvent->GetRunNumber();
129
130	// Derive trigger
131	UInt_t trigger = fInputHandler->IsEventSelected();
132	Bool_t isMBTrigger = trigger & AliVEvent::kMB;
133	Bool_t isCentralTrigger = trigger & AliVEvent::kCentral;
134	Bool_t isSemicentralTrigger = trigger & AliVEvent::kSemiCentral;
135	Bool_t isEMCALTrigger = trigger & AliVEvent::kEMCEJE;
136
137	// Check if MC information is available
138	Bool_t mcthere = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) != NULL;
139	if(mcthere){
140	fTrackCuts->SetMCEvent(fMCEvent);
141	fSignalCuts->SetMCEvent(fMCEvent);
142	}
143
144	// Get Primary Vertex
145	const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
146
147	// Get centrality
148	Float_t centrality = -1.;
149	TString beamtype = (dynamic_cast<AliESDEvent *>(fInputEvent))->GetESDRun()->GetBeamType();
150	if(!beamtype.CompareTo("Pb-Pb") \|\| !beamtype.CompareTo("A-A")){
151	// Heavy ion run
152	AliDebug(1, "Heavy-Ion event\n");
153	AliCentrality *hicent = fInputEvent->GetCentrality();
154	centrality = hicent->GetCentralityPercentile("V0M");
155	}
156
157	// Common variables
158	Double_t charge, eta, phi, momentum, transversemomentum;
159	Int_t source;
160
161	// Monte-Carlo loop
162	if(mcthere){
163	AliMCParticle * mcpart;
164	for(Int_t itrk = 0; itrk < fMCEvent->GetNumberOfTracks(); itrk++){
165	mcpart = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrk));
166	if(!mcpart) continue;
167	if(!fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mcpart)) continue;
168	// Kinematics
169	charge = mcpart->Charge() > 0. ? 1. : -1.;
170	momentum = mcpart->P() * charge;
171	transversemomentum = mcpart->Pt() * charge;
172	eta = mcpart->Eta();
173	phi = mcpart->Phi();
174	Int_t pdg = mcpart->Particle()->GetPdgCode();
175
176	// Get Production Vertex in radial direction
177	Double_t vx = mcpart->Particle()->Vx(),
178	vy = mcpart->Particle()->Vy();
179	Double_t productionVertex = TMath::Sqrt(vxvx+vyvy);
180
181	// Get Mother PDG code of the particle
182	Int_t motherPdg = 0;
183	Int_t motherlabel = mcpart->Particle()->GetFirstMother();
184	if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){
185	AliMCParticle mother = dynamic_cast<AliMCParticle >(fMCEvent->GetTrack(motherlabel));
186	if(mother) motherPdg = mother->Particle()->GetPdgCode();
187	}
188
189	// derive source
190	source = 5;
191	if(fSignalCuts->IsCharmElectron(mcpart)) source = 0;
192	else if(fSignalCuts->IsBeautyElectron(mcpart)) source = 1;
193	else if(fSignalCuts->IsGammaElectron(mcpart)) source = 2;
194	else if(fSignalCuts->IsNonHFElectron(mcpart)) source = 3;
195	else if(TMath::Abs(pdg) == 11) source = 4;
196	else source = 5;
197
198	(*fDebugTree) << "MCDebug"
199	<< "centrality=" << centrality
200	<< "MBtrigger=" << isMBTrigger
201	<< "CentralTrigger=" << isCentralTrigger
202	<< "SemicentralTrigger=" << isSemicentralTrigger
203	<< "EMCALtrigger=" << isEMCALTrigger
204	<< "run=" << run
205	<< "p=" << momentum
206	<< "pt=" << transversemomentum
207	<< "eta=" << eta
208	<< "phi=" << phi
209	<< "pdg=" << pdg
210	<< "ProductionVertex=" << productionVertex
211	<< "motherPdg=" << motherPdg
212	<< "source=" << source
213	<< "\n";
214	}
215	}
216
217	AliESDtrack *track;
218	for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
219	// fill the tree
220	track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrack));
221	if(!track) continue;
222	// Cut track (Only basic track cuts)
223	if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
224	// Debug streaming of PID-related quantities
225	Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
226	Double_t nSigmaTPC = pid->NumberOfSigmasTPC(track, AliPID::kElectron);
227	if(TMath::Abs(nSigmaTOF) > 5) continue;
228	// we are not interested in tracks which are more than 5 sigma away from the electron hypothesis in either TOF or TPC
229	Double_t tPCdEdx = track->GetTPCsignal();
230	// Signal, source and MCPID
231	Bool_t signal = kTRUE;
232	source = 5;
233	if(mcthere){
234	// Signal
235	AliMCParticle *mctrack;
236	if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))){
237	if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
238	}
239
240	// Source
241	if(fSignalCuts->IsCharmElectron(track)) source = 0;
242	else if(fSignalCuts->IsBeautyElectron(track)) source = 1;
243	else if(fSignalCuts->IsGammaElectron(track)) source = 2;
244	else if(fSignalCuts->IsNonHFElectron(track)) source = 3;
245	else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)) source = 4;
246	else source = 5;
247	}
248	// Get V0 tag (if available)
249	Int_t v0pid = -1;
250	if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron;
251	else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion;
252	else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton;
253	// Kinematics
254	charge = track->Charge() > 0 ? 1. : -1.;
255	eta = track->Eta();
256	phi = track->Phi();
257	momentum = track->P() * charge;
258	transversemomentum = track->Pt() * charge;
259	// ITS number of clusters
260	UChar_t nclustersITS = track->GetITSclusters(NULL);
261	Double_t chi2matching = track->GetChi2TPCConstrainedVsGlobal(dynamic_cast<const AliESDVertex *>(vertex));
262	Double_t chi2PerClusterITS = 0.0;
263	if (nclustersITS != 0) chi2PerClusterITS = track->GetITSchi2() / Float_t(nclustersITS);
264	// TPC number of clusters (different definitions)
265	UChar_t nclustersTPC = track->GetTPCncls();
266	UChar_t nclustersTPCPID = track->GetTPCsignalN();
267	UChar_t nfindableTPC = track->GetTPCNclsF();
268	Double_t clusterRatioTPC = 0.0;
269	if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPC = static_cast<Double_t>(nclustersTPC)/static_cast<Double_t>(nfindableTPC);
270	UChar_t nclustersTPCshared = 0;
271	Float_t ncrossedRowsTPC = track->GetTPCCrossedRows();
272	const TBits &sharedTPC = track->GetTPCSharedMap();
273	for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
274	// TRD clusters and tracklets
275	UChar_t nclustersTRD = track->GetTRDncls();
276	UChar_t ntracklets = track->GetTRDntrackletsPID();
277	// ITS and TRD acceptance maps
278	UChar_t hasClusterITS[6], hasTrackletTRD[6];
279	UChar_t itsPixel = track->GetITSClusterMap();
280	for(Int_t icl = 0; icl < 6; icl++) hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
281	for(Int_t itl = 0; itl < 6; itl++){
282	Int_t nSliceNonZero = 0;
283	for(Int_t islice = 0; islice < 8; islice++){
284	if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++;
285	}
286	hasTrackletTRD[itl] = nSliceNonZero ? 1 : 0;
287	}
288	// TRD PID
289	Double_t pidprobs[5];
290	track->GetTRDpid(pidprobs);
291	Double_t likeEleTRD = pidprobs[0];
292	Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]);
293	// DCA
294	Float_t b[2] = {0.,0.};
295	Float_t bCov[3] = {0.,0.,0.};
296	track->GetImpactParameters(b,bCov);
297	Double_t dca = TMath::Sqrt(b[0]b[0]+b[1]b[1]); // impact parameter space
298	Double_t dcaSR=0, dcaSZ=0;
299	if(bCov[0]>0) dcaSR = b[0]/TMath::Sqrt(bCov[0]); // normalised impact parameter xy
300	if(bCov[2]>0) dcaSZ = b[1]/TMath::Sqrt(bCov[2]); // normalised impact parameter z
301	Double_t dcaS = AliESDtrackCuts::GetSigmaToVertex(track); // n_sigma
302	// Vertex
303	Double_t vtx[3];
304	vertex->GetXYZ(vtx);
305	Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
306	// Fill Tree
307	(*fDebugTree) << "PIDdebug"
308	<< "centrality=" << centrality
309	<< "MBtrigger=" << isMBTrigger
310	<< "CentralTrigger=" << isCentralTrigger
311	<< "SemicentralTrigger=" << isSemicentralTrigger
312	<< "EMCALtrigger=" << isEMCALTrigger
313	<< "signal=" << signal
314	<< "source=" << source
315	<< "v0pid=" << v0pid
316	<< "run=" << run
317	<< "p=" << momentum
318	<< "pt=" << transversemomentum
319	<< "eta=" << eta
320	<< "phi=" << phi
321	<< "ntracklets=" << ntracklets
322	<< "nclustersTPC=" << nclustersTPC
323	<< "nclustersTPCPID=" << nclustersTPCPID
324	<< "nclustersTPCshared=" << nclustersTPCshared
325	<< "ncrossedRowsTPC=" << ncrossedRowsTPC
326	<< "clusterRatioTPC=" << clusterRatioTPC
327	<< "nclustersITS=" << nclustersITS
328	<< "nclusters=" << nclustersTRD
329	<< "chi2matching=" << chi2matching
330	<< "chi2PerClusterITS=" << chi2PerClusterITS
331	<< "its0=" << hasClusterITS[0]
332	<< "its1=" << hasClusterITS[1]
333	<< "its2=" << hasClusterITS[2]
334	<< "its3=" << hasClusterITS[3]
335	<< "its4=" << hasClusterITS[4]
336	<< "its5=" << hasClusterITS[5]
337	<< "trd0=" << hasTrackletTRD[0]
338	<< "trd1=" << hasTrackletTRD[1]
339	<< "trd2=" << hasTrackletTRD[2]
340	<< "trd3=" << hasTrackletTRD[3]
341	<< "trd4=" << hasTrackletTRD[4]
342	<< "trd5=" << hasTrackletTRD[5]
343	<< "TOFsigmaEl=" << nSigmaTOF
344	<< "TPCsigmaEl=" << nSigmaTPC
345	<< "TPCdEdx=" << tPCdEdx
346	<< "TRDlikeEl=" << likeEleTRD
347	<< "TRDlikeEln=" << likeEleTRDn
348	<< "dcaR=" << b[0]
349	<< "dcaZ=" << b[1]
350	<< "dca=" << dca
351	<< "dcaSR=" << dcaSR
352	<< "dcaSZ=" << dcaSZ
353	<< "dcaS=" << dcaS
354	<< "vx=" << vtx[0]
355	<< "vy=" << vtx[1]
356	<< "vz=" << vtx[2]
357	<< "ncontrib=" << ncontrib
358	<< "\n";
359	}
360	}
361
362
363	void AliHFEdebugTreeTask::SetFileName(const char *filename){ fFilename = filename; }