#include "AliCentrality.h"
#include "AliAODTrack.h"
#include "AliAODEvent.h"
+#include "AliAODPid.h"
#include "AliHFEcuts.h"
#include "AliHFEextraCuts.h"
#include "AliInputEventHandler.h"
//printf("test\n");
fDebugTree = new TTreeSRedirector(fFilename.Data());
- //printf("testa\n");
+ // printf("testa\n");
fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
//printf("testb\n");
AliPIDResponse *pid = NULL;
AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
if(handler){
- //printf("testb\n");
+// printf("testb\n");
pid = handler->GetPIDResponse();
} else {
AliError("No Handler");
}
if(!pid){
- printf("testc\n");
+ // printf("testc\n");
AliError("No PID response");
return;
}
return;
}
+ AliAODTrack copyTrack;
+
// MC info
Bool_t mcthere = kTRUE;
AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
- if(!aodE){
+ if(!aodE){
+ // printf("testd\n");
AliError("No AOD Event");
return;
}
fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
if(!fAODMCHeader){
- mcthere = kFALSE;
- return;
+ mcthere = kFALSE;
+ // return;
}
fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
if(!fAODArrayMCInfo){
- mcthere = kFALSE;
- return;
+ mcthere = kFALSE;
+ // return;
}
else {
fSignalCuts->SetMCAODInfo(fAODArrayMCInfo);
}
fExtraCuts->SetRecEventInfo(fInputEvent);
+
+
+
// Get run number
Int_t run = fInputEvent->GetRunNumber();
AliCentrality *hicent = fInputEvent->GetCentrality();
centrality = hicent->GetCentralityPercentile("V0M");
+ // Store event selection variables
+ (*fDebugTree) << "EventDebug"
+ << "Centrality=" << centrality
+ << "VertexZ=" << vtx[2]
+ << "NumberOfContributors=" << ncontrib
+ << "run=" << run
+ << "\n";
+
+
+
// Look for kink mother
Int_t numberofvertices = aodE->GetNumberOfVertices();
Double_t listofmotherkink[numberofvertices];
}
}
//printf("Number of kink mother in the events %d\n",numberofmotherkink);
+
+ //
+ // Loop on MC tracks only
+ //
+ AliAODMCParticle *mctrack = NULL;
+ // Monte-Carlo info
+ Double_t eR,vx,vy,vz;
+ Double_t chargemc, etamc, phimc, momentummc, transversemomentummc;
+ Int_t source,pdg,signal;
+ if(mcthere) {
+ for(Int_t itrack = 0; itrack < fAODArrayMCInfo->GetEntriesFast(); itrack++) {
+ mctrack = (AliAODMCParticle *)(fAODArrayMCInfo->At(itrack));
+ if(!mctrack) continue;
+ signal = 0;
+ if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) signal = 1;
+ // Kinematics
+ chargemc = mctrack->Charge() > 0. ? 1. : -1.;
+ momentummc = mctrack->P() * chargemc;
+ transversemomentummc = mctrack->Pt() * chargemc;
+ etamc = mctrack->Eta();
+ phimc = mctrack->Phi();
+ pdg = mctrack->GetPdgCode();
+
+ // Get Production Vertex in radial direction
+ vx = mctrack->Xv();
+ vy = mctrack->Yv();
+ vz = mctrack->Zv();
+ eR = TMath::Sqrt(vx*vx+vy*vy);
+
+ // Get Mother PDG code of the particle
+ Int_t motherPdg = 0;
+ Int_t motherlabel = TMath::Abs(mctrack->GetMother());
+ if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
+ AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
+ if(mother) motherPdg = mother->GetPdgCode();
+ }
+
+ // derive source
+ source = 5;
+ if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
+ else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
+ else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
+ else if(fSignalCuts->IsNonHFElectron(mctrack)) 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=" << momentummc
+ << "pt=" << transversemomentummc
+ << "eta=" << etamc
+ << "phi=" << phimc
+ << "pdg=" << pdg
+ << "px=" << vx
+ << "py=" << vy
+ << "pz=" << vz
+ << "ProductionVertex=" << eR
+ << "motherPdg=" << motherPdg
+ << "source=" << source
+ << "\n";
+ }
+ }
// Common variables
Double_t charge, eta, phi, momentum, momentumTPC, transversemomentum;
//
TArrayI *arraytrack = new TArrayI(fInputEvent->GetNumberOfTracks());
+ Int_t counterdc=0;
AliAODTrack *track = 0x0;
- AliAODMCParticle *mctrack = NULL;
for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
// fill the tree
track = dynamic_cast<AliAODTrack *>(fInputEvent->GetTrack(itrack));
if(!track) continue;
// Cut track (Only basic track cuts)
- //printf("testv\n");
+ // printf("testv\n");
if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
//
//printf("testu\n");
Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
Double_t nSigmaTPC = pid->NumberOfSigmasTPC(track, AliPID::kElectron);
Double_t tPCdEdx = track->GetDetPid() ? track->GetDetPid()->GetTPCsignal() : 0.;
+ // Eta correction
+ copyTrack.~AliAODTrack();
+ new(©Track) AliAODTrack(*track);
+ if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(©Track, static_cast<Double_t>(ncontrib),AliHFEpidObject::kAODanalysis);
+ if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(©Track, AliHFEpidObject::kAODanalysis);
+ Double_t nSigmaTPCcorr = pid->NumberOfSigmasTPC(©Track, AliPID::kElectron);
+ Double_t tPCdEdxcorr = copyTrack.GetDetPid() ? copyTrack.GetDetPid()->GetTPCsignal() : 0.;
// Kinematics
charge = track->Charge() > 0 ? 1. : -1.;
if((status & AliESDtrack::kITSrefit) == AliESDtrack::kITSrefit) itsrefit = 1;
Int_t tpcrefit=0;
if((status & AliESDtrack::kTPCrefit) == AliESDtrack::kTPCrefit) tpcrefit = 1;
+ Int_t tofpid=0;
+ if((status & AliESDtrack::kTOFpid) == AliESDtrack::kTOFpid) tofpid = 1;
// ITS number of clusters
UChar_t nclustersITS = track->GetITSNcls();
// TRD clusters and tracklets
UChar_t nclustersTRD = track->GetTRDncls();
UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID();
- Int_t nslicesTRD = track->GetNumberOfTRDslices();
+ Int_t nslicesTRD = track->GetDetPid()->GetTRDnSlices();
Int_t chi2TRD = track->GetTRDchi2();
+ AliAODPid* aodpid= track->GetDetPid();
+ Double_t* arraytrdsignals;
+ arraytrdsignals=aodpid->GetTRDslices();
+ Int_t nslicetemp=0;
+ Int_t trdlayer[6];
+ for(Int_t iplane = 0; iplane < 6; iplane++){
+ trdlayer[iplane]=0;
+ }
+
+ for(Int_t iplane = 0; iplane < 6; iplane++){
+ nslicetemp=0;
+ for(Int_t b=(iplane*8);b<((iplane*8)+8);b++)
+ {
+ if(ntrackletsTRDPID>0)
+ {
+ if(arraytrdsignals[b]>0.001) nslicetemp++;
+ }
+ }
+ if(nslicetemp > 0) trdlayer[iplane]=1;
+ }
+
// ITS and TRD acceptance maps
UChar_t itsPixel = track->GetITSClusterMap();
Bool_t statusL0 = kFALSE;
// Double counted
Int_t doublec = 0;
- for(Int_t l=0; l < itrack; l++){
+ for(Int_t l=0; l < counterdc; l++){
Int_t iTrack2 = arraytrack->At(l);
if(iTrack2==id) doublec=1;
}
//printf("Doublec %d\n",doublec);
// Add the id at this place
- arraytrack->AddAt(id,itrack);
+ arraytrack->AddAt(id,counterdc);
+ counterdc++;
// Filter
Int_t filter[20];
//printf("track\n");
- // Monte-Carlo info
- Double_t eR,vx,vy,vz;
- Double_t chargemc, etamc, phimc, momentummc, transversemomentummc;
- Int_t source,pdg,signal;
signal = 0;
if(mcthere){
Int_t label = TMath::Abs(track->GetLabel());
if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) signal = 1;
// Kinematics
chargemc = mctrack->Charge() > 0. ? 1. : -1.;
- momentummc = mctrack->P() * charge;
- transversemomentummc = mctrack->Pt() * charge;
+ momentummc = mctrack->P() * chargemc;
+ transversemomentummc = mctrack->Pt() * chargemc;
etamc = mctrack->Eta();
phimc = mctrack->Phi();
pdg = mctrack->GetPdgCode();
<< "phi=" << phi
<< "itsrefit=" << itsrefit
<< "tpcrefit=" << tpcrefit
+ << "tofpid=" << tofpid
<< "nclustersTPC=" << nclustersTPCfit
<< "nclustersTPCall=" << nclustersTPCall
<< "nclustersTPCPID=" << nclustersTPCPID
<< "nclustersTRD=" << nclustersTRD
<< "ntrackletsTRD=" << ntrackletsTRDPID
<< "nslicesTRD=" << nslicesTRD
+ << "trd0=" << trdlayer[0]
+ << "trd1=" << trdlayer[1]
+ << "trd2=" << trdlayer[2]
+ << "trd3=" << trdlayer[3]
+ << "trd4=" << trdlayer[4]
+ << "trd5=" << trdlayer[5]
<< "chi2TRD=" << chi2TRD
<< "statusITS0=" << statusL0
<< "statusITS1=" << statusL1
<< "TOFsigmaEl=" << nSigmaTOF
<< "TPCsigmaEl=" << nSigmaTPC
+ << "TPCsigmaElcorr=" << nSigmaTPCcorr
<< "TPCdEdx=" << tPCdEdx
+ << "TPCdEdxcorr=" << tPCdEdxcorr
<< "dcaR=" << dcaxy
<< "dcaZ=" << dcaz
<< "kinkdaughter=" << kink