[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 "AliHFEpidTRD.h"
#include "TTreeStream.h"

#include "AliHFEdebugTreeTask.h"

ClassImp(AliHFEdebugTreeTask)

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

}

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

}

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

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();

  fTRDpid = new AliHFEpidTRD("QAtrdPID");

}

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");
    return;
  }
  if(!fInputEvent) {
    AliError("No Input event");
    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.;
  AliESDEvent *event = (dynamic_cast<AliESDEvent *>(fInputEvent));
  if(!event) return;
  TString beamtype = event->GetBeamType();
  //printf("Beam type %s\n",(const char*)beamtype);
  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 ntrackletsTRDPID = 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;
    Double_t trddEdxSum[6];
    for(Int_t a=0;a<6;a++) { trddEdxSum[a]= 0.;}
    for(Int_t itl = 0; itl < 6; itl++){
	Int_t nSliceNonZero = 0;
        trddEdxSum[itl] = track->GetTRDslice(itl, 0); // in new reconstruction slice 0 contains the total charge
      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]);
    Double_t trdtruncmean1 = fTRDpid->GetTRDSignalV1(track, 0.6);
    Double_t trdtruncmean2 = fTRDpid->GetTRDSignalV2(track, 0.6);

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