#include <TParticle.h>
#include <AliLog.h>
-#include <AliStack.h>
+#include <AliMCEvent.h>
#include <AliGenEventHeader.h>
#include <AliAODMCParticle.h>
//_______________________________________________________________________________________________
AliHFEmcQA::AliHFEmcQA() :
- fStack(NULL)
+ fMCEvent(NULL)
,fMCHeader(NULL)
,fMCArray(NULL)
,fQAhistos(NULL)
//_______________________________________________________________________________________________
AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
TObject(p)
- ,fStack(NULL)
+ ,fMCEvent(NULL)
,fMCHeader(NULL)
,fMCArray(NULL)
,fQAhistos(p.fQAhistos)
{
// create histograms
- if (kquark != kCharm && kquark != kBeauty) {
+ if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
AliDebug(1, "This task is only for heavy quark QA, return\n");
return;
}
kqTypeLabel[kDeHadron]="bDeHadron";
kqTypeLabel[kElectron]="be";
kqTypeLabel[kElectron2nd]="bce";
+ } else if (kquark == kOthers){
+ kqTypeLabel[kGamma-4]="gammae";
+ kqTypeLabel[kPi0-4]="pi0e";
+ kqTypeLabel[kElse-4]="elsee";
+ kqTypeLabel[kMisID-4]="miside";
}
-
-
+/*
+ const Double_t kPtbound[2] = {0.001, 50.};
+ Int_t iBin[1];
+ iBin[0] = 100; // bins in pt
+ Double_t* binEdges[1];
+ binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
+ */
+
TString hname;
+ if(kquark == kOthers){
+ for (Int_t iqType = 0; iqType < 4; iqType++ ){
+ hname = hnopt+"Pt_"+kqTypeLabel[iqType];
+ //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
+ fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",500,0,50);
+ hname = hnopt+"Y_"+kqTypeLabel[iqType];
+ fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
+ hname = hnopt+"Eta_"+kqTypeLabel[iqType];
+ fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
+ // Fill List
+ if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
+ }
+ return;
+ }
for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
hname = hnopt+"Pt_"+kqTypeLabel[iqType];
- fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",250,0,50);
+ //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
+ fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",500,0,50);
hname = hnopt+"Y_"+kqTypeLabel[iqType];
fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
hname = hnopt+"Eta_"+kqTypeLabel[iqType];
}
Int_t iq = kquark - kCharm;
- if (iTrack < 0) {
- AliDebug(1, "Stack label is negative, return\n");
+ if (iTrack < 0 || !part) {
+ AliDebug(1, "Stack label is negative or no mcparticle, return\n");
return;
}
- //TParticle *part = fStack->Particle(iTrack);
+ AliMCParticle *mctrack = NULL;
Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
// select heavy hadron or not fragmented heavy quark
iLabel = iTrack;
} else{ // in case of heavy hadron, start to search for mother heavy parton
iLabel = part->GetFirstMother();
- if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
+ if (iLabel>-1) { partMother = mctrack->Particle(); }
else {
AliDebug(1, "Stack label is negative, return\n");
return;
Bool_t isSameString = kTRUE;
for (Int_t i=1; i<fgkMaxIter; i++){
iLabel = iLabel - 1;
- if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
+ if (iLabel>-1) { partMother = mctrack->Particle(); }
else {
AliDebug(1, "Stack label is negative, return\n");
return;
}
Int_t iq = kquark - kCharm;
- //TParticle* mcpart = fStack->Particle(iTrack);
+ if(!mcpart){
+ AliDebug(1, "no mc particle, return\n");
+ return;
+ }
Int_t iLabel = mcpart->GetFirstMother();
if (iLabel<0){
Int_t pdgcode = mcpart->GetPdgCode();
Int_t pdgcodeCopy = pdgcode;
+ AliMCParticle *mctrack = NULL;
+
// if the mother is charmed hadron
Bool_t isDirectCharm = kFALSE;
if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
return;
}
// if there is an ancester
- TParticle* mother = fStack->Particle(jLabel);
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
+ TParticle* mother = mctrack->Particle();
Int_t motherPDG = mother->GetPdgCode();
for (Int_t j=0; j<fNparents; j++){
{
// decay electron kinematics
- if (kquark != kCharm && kquark != kBeauty) {
+ if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
AliDebug(1, "This task is only for heavy quark QA, return\n");
return;
}
Int_t iq = kquark - kCharm;
Bool_t isFinalOpenCharm = kFALSE;
-/*
- if (iTrack < 0) {
- AliDebug(1, "Stack label is negative, return\n");
- return;
+ if(!mcpart){
+ AliDebug(1, "no mcparticle, return\n");
+ return;
}
- */
- //TParticle* mcpart = fStack->Particle(iTrack);
+ if(kquark==kOthers){
+ Int_t esource = -1;
+ if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
+ else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
+ if(esource==0|| esource==1 || esource==2 || esource==3){
+ fHist[iq][esource][icut].fPt->Fill(mcpart->Pt());
+ fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+ fHist[iq][esource][icut].fEta->Fill(mcpart->Eta());
+ return;
+ }
+ else {
+ AliDebug(1, "e source is out of defined ranges, return\n");
+ return;
+ }
+ }
if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
return;
}
- TParticle *partMother = fStack->Particle(iLabel);
+ AliMCParticle *mctrack = NULL;
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
+ TParticle *partMother = mctrack->Particle();
TParticle *partMotherCopy = partMother;
Int_t maPdgcode = partMother->GetPdgCode();
Int_t maPdgcodeCopy = maPdgcode;
// get mc primary vertex
+ /*
TArrayF mcPrimVtx;
if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
// calculated production vertex to primary vertex (in xy plane)
Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
+ */
+ Float_t decayLxy = 0;
// if the mother is charmed hadron
Bool_t isMotherDirectCharm = kFALSE;
}
// if there is an ancester
- TParticle* grandMa = fStack->Particle(jLabel);
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
+ TParticle* grandMa = mctrack->Particle();
Int_t grandMaPDG = grandMa->GetPdgCode();
for (Int_t j=0; j<fNparents; j++){
Int_t iq = kquark - kCharm;
Bool_t isFinalOpenCharm = kFALSE;
+ if(!mcpart){
+ AliDebug(1, "no mcparticle, return\n");
+ return;
+ }
+
if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
// mother
AliDebug(1, "Stack label is negative, return\n");
return;
}
- TParticle *heavysMother = fStack->Particle(motherlabel);
+ AliMCParticle *mctrack = NULL;
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
+ TParticle *heavysMother = mctrack->Particle();
motherpdg = heavysMother->GetPdgCode();
grandmotherlabel = heavysMother->GetFirstMother();
AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
{
// mothertype -1 means this heavy quark coming from hard vertex
- TParticle *afterinitialrad1 = fStack->Particle(4);
- TParticle *afterinitialrad2 = fStack->Particle(5);
+ AliMCParticle *mctrack1 = NULL;
+ AliMCParticle *mctrack2 = NULL;
+ if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
+ if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
+ TParticle *afterinitialrad1 = mctrack1->Particle();
+ TParticle *afterinitialrad2 = mctrack2->Particle();
motherlabel = -1;
{
// mothertype -2 means this heavy quark coming from initial state
+ AliMCParticle *mctrack = NULL;
if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
- TParticle *heavysMother = fStack->Particle(inputmotherlabel);
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
+ TParticle *heavysMother = mctrack->Particle();
motherID = heavysMother->GetPdgCode();
mothertype = -2; // there is mother before initial state radiation
motherlabel = inputmotherlabel;
{
// mothertype 2 means this heavy quark coming from final state
+ AliMCParticle *mctrack = NULL;
if (inputmotherlabel > 5){ // mother exist after hard scattering
- TParticle *heavysMother = fStack->Particle(inputmotherlabel);
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
+ TParticle *heavysMother = mctrack->Particle();
motherID = heavysMother->GetPdgCode();
mothertype = 2; //
motherlabel = inputmotherlabel;
Int_t origin = -1;
Bool_t isFinalOpenCharm = kFALSE;
+ if(!mcpart){
+ AliDebug(1, "Stack label is negative or no mcparticle, return\n");
+ return -1;
+ }
+
// mother
Int_t iLabel = mcpart->GetMother();
if (iLabel<0){
Int_t origin = -1;
Bool_t isFinalOpenCharm = kFALSE;
+ if(!mcpart){
+ AliDebug(1, "no mcparticle, return\n");
+ return -1;
+ }
+
Int_t iLabel = mcpart->GetFirstMother();
if (iLabel<0){
AliDebug(1, "Stack label is negative, return\n");
return -1;
}
- TParticle *partMother = fStack->Particle(iLabel);
+ AliMCParticle *mctrack = NULL;
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
+ TParticle *partMother = mctrack->Particle();
Int_t maPdgcode = partMother->GetPdgCode();
// if the mother is charmed hadron
}
// if there is an ancester
- TParticle* grandMa = fStack->Particle(jLabel);
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
+ TParticle* grandMa = mctrack->Particle();
Int_t grandMaPDG = grandMa->GetPdgCode();
for (Int_t j=0; j<fNparents; j++){
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff