Updates on single muon HF analysis (Shuang)
authorlaphecet <laurent.aphecetche@subatech.in2p3.fr>
Fri, 16 Jan 2015 09:26:09 +0000 (10:26 +0100)
committerlaphecet <laurent.aphecetche@subatech.in2p3.fr>
Fri, 16 Jan 2015 09:26:09 +0000 (10:26 +0100)
PWG/muon/AliAnalysisTaskSEMuonsHF.cxx
PWG/muon/AliAnalysisTaskSEMuonsHF.h
PWG/muon/AliMuonInfoStoreMC.cxx
PWG/muon/AliMuonInfoStoreMC.h
PWG/muon/AliMuonInfoStoreRD.cxx
PWG/muon/AliMuonInfoStoreRD.h
PWG/muon/AliMuonsHFHeader.cxx
PWG/muon/AliMuonsHFHeader.h

index 4e23217..aad5be9 100644 (file)
@@ -33,6 +33,7 @@
 
 #include "AliESDEvent.h"
 #include "AliAODEvent.h"
+#include "AliMCEvent.h"
 #include "AliESDHeader.h"
 #include "AliAODHeader.h"
 #include "AliESDMuonTrack.h"
@@ -57,6 +58,7 @@ AliAnalysisTaskSE(),
 fAnaMode(0),
 fIsOutputTree(kFALSE),
 fIsMC(kFALSE),
+fIsFull(kFALSE),
 fCutsMuon(0x0),
 fCutsDimu(0x0),
 fHeader(0),
@@ -75,6 +77,7 @@ AliAnalysisTaskSE(name),
 fAnaMode(0),
 fIsOutputTree(kFALSE),
 fIsMC(kFALSE),
+fIsFull(kFALSE),
 fCutsMuon(new AliMuonTrackCuts(cutsMuon)),
 fCutsDimu(new AliMuonPairCuts(cutsDimu)),
 fHeader(0),
@@ -108,12 +111,11 @@ void AliAnalysisTaskSEMuonsHF::Init()
   // Initialization
   // Setting and initializing the running mode and status
 
-  AliMuonsHFHeader::SetAnaMode(fAnaMode);
-  AliMuonsHFHeader::SetIsMC(fIsMC);
   if (!fHeader) {
     fHeader = new AliMuonsHFHeader();
     fHeader->SetName(AliMuonsHFHeader::StdBranchName());
-  }
+  } fHeader->SetAnaMode(fAnaMode);
+    fHeader->SetIsMC(fIsMC);
 
   if (!fMuonClArr) {
     if (fIsMC) { 
@@ -154,6 +156,7 @@ void AliAnalysisTaskSEMuonsHF::UserCreateOutputObjects()
 
   fCutsMuon->Print("mask");
   fCutsDimu->Print("mask");
+
   return;
 }
 
@@ -164,11 +167,11 @@ void AliAnalysisTaskSEMuonsHF::UserExec(Option_t *)
   // muon event header & (di)muon info store
 
   if (fIsMC) {
-    if (MCEvent()) {
-      if (MCEvent()->GetNumberOfTracks()<=0)
+    if (fInputHandler->MCEvent()) {
+      if (fInputHandler->MCEvent()->GetNumberOfTracks()<=0)
            { AliError("MC event not found. Nothing done!"); return; }
     } else { AliError("MC event not found. Nothing done!"); return; }
-  }
+  } 
 
   if ( !fCutsMuon) { AliError("AliMuonTrackCuts should be loaded!"); return;}
   if ((!fCutsDimu) && (fAnaMode!=1)) { AliError("AliMuonPairCuts should be loaded!"); return; }
@@ -188,7 +191,7 @@ void AliAnalysisTaskSEMuonsHF::UserExec(Option_t *)
     ntrks = esd->GetNumberOfMuonTracks();
   } if (fIsOutputTree) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
 
-  fHeader->SetEventInfo(fInputHandler, MCEvent());
+  fHeader->SetEventInfo(fInputHandler);
   fHeader->FillHistosEvnH(fListOutput);
 
   fMuonClArr->Delete();
@@ -204,13 +207,13 @@ void AliAnalysisTaskSEMuonsHF::UserExec(Option_t *)
     if (aod) {
       trkAOD = (AliAODTrack*)aod->GetTrack(itrk);
       if (!trkAOD->IsMuonTrack())        { trkAOD=0; continue; }
-      if (fIsMC) muonMC = new AliMuonInfoStoreMC(trkAOD,MCEvent(),fCutsMuon->GetSelectionMask(trkAOD));
+      if (fIsMC) muonMC = new AliMuonInfoStoreMC(trkAOD,fInputHandler->MCEvent(),fCutsMuon->GetSelectionMask(trkAOD),fIsFull);
       else muonRD = new AliMuonInfoStoreRD(trkAOD,fCutsMuon->GetSelectionMask(trkAOD));
       trkAOD = 0;
     } else {
       trkESD = (AliESDMuonTrack*)esd->GetMuonTrack(itrk);
       if (!trkESD->ContainTrackerData()) { trkESD=0; continue; }
-      if (fIsMC) muonMC = new AliMuonInfoStoreMC(trkESD,MCEvent(),fCutsMuon->GetSelectionMask(trkESD));
+      if (fIsMC) muonMC = new AliMuonInfoStoreMC(trkESD,fInputHandler->MCEvent(),fCutsMuon->GetSelectionMask(trkESD),fIsFull);
       else muonRD = new AliMuonInfoStoreRD(trkESD,fCutsMuon->GetSelectionMask(trkESD));
       trkESD = 0;
     } if (muonRD) {
@@ -252,13 +255,14 @@ void AliAnalysisTaskSEMuonsHF::UserExec(Option_t *)
         if (fIsOutputTree) new(dimuRef[countN++]) AliDimuInfoStoreMC(*dimuMC);
       }
 
-      if (dimuRD) { delete dimuRD; dimuRD=0; }
-      if (dimuMC) { delete dimuMC; dimuMC=0; }
+      if (dimuRD) { delete dimuRD; dimuRD = 0; }
+      if (dimuMC) { delete dimuMC; dimuMC = 0; }
     }  // end 2nd loop of muon tracks
   }  // end 1st loop of muon tracks
 
   aod = 0; esd = 0;
   PostData(1, fListOutput);
+
   return;
 }
 
@@ -274,9 +278,10 @@ void AliAnalysisTaskSEMuonsHF::Terminate(Option_t *)
 
 void AliAnalysisTaskSEMuonsHF::NotifyRun()
 {
-  // Notify of the current run number
+  // Notify of the input handler 
 
   if (fCutsMuon) fCutsMuon->SetRun(fInputHandler);
   if (fCutsDimu) fCutsDimu->SetRun(fInputHandler);
+
   return;
 }
index f942ad3..d345566 100644 (file)
@@ -36,32 +36,34 @@ class AliAnalysisTaskSEMuonsHF : public AliAnalysisTaskSE {
   virtual void Terminate(Option_t *opt);
   virtual void NotifyRun();
 
-  void SetAnaMode(Int_t mode)      { fAnaMode      = ((mode>=0 && mode<3) ? mode : 0); }
-  void SetIsOutputTree(Bool_t ist) { fIsOutputTree = ist;                              }
-  void SetUseMC(Bool_t isMC)       { fIsMC         = isMC;                             }
+  void SetAnaMode(Int_t mode)      { fAnaMode      = ((mode>=0 && mode<=2) ? mode : 0); }
+  void SetIsOutputTree(Bool_t ist) { fIsOutputTree = ist;                               }
+  void SetUseMC(Bool_t isMC)       { fIsMC         = isMC;                              }
+  void SetIsFull(Bool_t isFull)    { fIsFull       = isFull;                            }
 
-  void SetEvsHCuts(Double_t cuts[5])  const { AliMuonsHFHeader::SetSelectionCuts(cuts);   }
+  void SetEvsHCuts(Double_t cuts[5])  const { AliMuonsHFHeader::SetSelectionCuts(cuts); }
 
  private:
 
   AliAnalysisTaskSEMuonsHF(const AliAnalysisTaskSEMuonsHF&);
   AliAnalysisTaskSEMuonsHF& operator=(const AliAnalysisTaskSEMuonsHF&);
 
-  Int_t fAnaMode;        // = 0, ana both single muon and dimuon
-                         // = 1, ana single muon
-                         // = 2, ana dimuon
-  Bool_t fIsOutputTree;  // flag used to switch on/off tree output
-  Bool_t fIsMC;          // flag of whether the input is MC
+  Int_t fAnaMode;       // = 0, ana both single muon and dimuon
+                        // = 1, ana single muon
+                        // = 2, ana dimuon
+  Bool_t fIsOutputTree; // flag used to switch on/off tree output
+  Bool_t fIsMC;         // flag to use MC
+  Bool_t fIsFull;       // flag to save the parton info in MC (PYTHIA)
 
-  AliMuonTrackCuts *fCutsMuon;  // single muon selection cuts
-  AliMuonPairCuts  *fCutsDimu;  // dimuon selection cuts
+  AliMuonTrackCuts *fCutsMuon; // single muon selection cuts
+  AliMuonPairCuts  *fCutsDimu; // dimuon selection cuts
 
-  AliMuonsHFHeader *fHeader;  // output for info at ev level
-  TClonesArray  *fMuonClArr;  // output clones array for single mu
-  TClonesArray  *fDimuClArr;  // output clones array for dimu
-  TList *fListOutput;         // output list of histos
+  AliMuonsHFHeader *fHeader;     // output for info at ev level
+  TClonesArray     *fMuonClArr;  // output clones array for single mu
+  TClonesArray     *fDimuClArr;  // output clones array for dimu
+  TList            *fListOutput; // output list of histos
 
-  ClassDef(AliAnalysisTaskSEMuonsHF, 7);
+  ClassDef(AliAnalysisTaskSEMuonsHF, 8);
 };
 
 #endif
index 30e5bee..e47d774 100644 (file)
@@ -31,7 +31,6 @@
 #include "AliAODTrack.h"
 #include "AliMuonInfoStoreRD.h"
 #include "AliMuonInfoStoreMC.h"
-#include "AliMCEvent.h"
 
 ClassImp(AliMuonInfoStoreMC)
 
@@ -48,7 +47,9 @@ fTrackPDGCode(0),
 fSource(-1),
 fNParents(0),
 fOscillation(kFALSE),
-fWeight(0.)
+fWeight(0.),
+fIsPhyPrim(kFALSE),
+fGeneratorIndex(-1)
 {
   //
   // default constructor
@@ -67,7 +68,9 @@ fTrackPDGCode(0),
 fSource(-1),
 fNParents(0),
 fOscillation(kFALSE),
-fWeight(0.)
+fWeight(0.),
+fIsPhyPrim(kFALSE),
+fGeneratorIndex(-1)
 {
   //
   // default constructor
@@ -88,7 +91,9 @@ fTrackPDGCode(0),
 fSource(-1),
 fNParents(0),
 fOscillation(kFALSE),
-fWeight(0.)
+fWeight(0.),
+fIsPhyPrim(kFALSE),
+fGeneratorIndex(-1)
 {
   //
   // default constructor
@@ -109,7 +114,9 @@ fTrackPDGCode(src.fTrackPDGCode),
 fSource(src.fSource),
 fNParents(src.fNParents),
 fOscillation(src.fOscillation),
-fWeight(src.fWeight)
+fWeight(src.fWeight),
+fIsPhyPrim(src.fIsPhyPrim),
+fGeneratorIndex(src.fGeneratorIndex)
 {
   //
   // copy constructor
@@ -130,17 +137,19 @@ AliMuonInfoStoreMC& AliMuonInfoStoreMC::operator=(const AliMuonInfoStoreMC &src)
   //
   // assignment constructor
   //
-  if(&src==this) return *this;
+  if (&src==this) return *this;
   AliMuonInfoStoreRD::operator=(src);
 
-  fIsFull       = src.fIsFull;
-  fLorentzP     = src.fLorentzP;
-  fTrackIndex   = src.fTrackIndex;
-  fTrackPDGCode = src.fTrackPDGCode;
-  fSource       = src.fSource;
-  fNParents     = src.fNParents;
-  fOscillation  = src.fOscillation;
-  fWeight       = src.fWeight;
+  fIsFull         = src.fIsFull;
+  fLorentzP       = src.fLorentzP;
+  fTrackIndex     = src.fTrackIndex;
+  fTrackPDGCode   = src.fTrackPDGCode;
+  fSource         = src.fSource;
+  fNParents       = src.fNParents;
+  fOscillation    = src.fOscillation;
+  fWeight         = src.fWeight;
+  fIsPhyPrim      = src.fIsPhyPrim;
+  fGeneratorIndex = src.fGeneratorIndex;
   for (Int_t i=5; i--;) {
     fParentIndex[i]   = src.fParentIndex[i];
     fParentPDGCode[i] = src.fParentPDGCode[i];
@@ -166,32 +175,36 @@ void AliMuonInfoStoreMC::SetMCInfoAOD(AliMCEvent *mcEvent, Int_t label)
 {
   // fill track MC info with AOD base
   fTrackIndex = label;
-  if (fTrackIndex<0) { fSource=5; return; }
+  if (fTrackIndex<0) { fSource = 5; return; }
 
   AliAODMCParticle *pMC = (AliAODMCParticle*)mcEvent->GetTrack(fTrackIndex);
   fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
 
   fTrackPDGCode = pMC->PdgCode();
-  if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; } 
+  if (TMath::Abs(fTrackPDGCode)!=13) { fSource = 4; return; } 
+
+  fGeneratorIndex = pMC->GetGeneratorIndex();
 
   Int_t lineM = pMC->GetMother();
-  if (lineM<0) { fSource=2; return; }
+  Bool_t isPhysicalPrimary = ((AliAODMCParticle*)mcEvent->GetTrack(lineM))->IsPhysicalPrimary();
+  if (isPhysicalPrimary) fIsPhyPrim = kTRUE;
+  if (lineM<0) { fSource = 2; return; }
 
   Bool_t isPrimary = ((AliAODMCParticle*)mcEvent->GetTrack(lineM))->IsPrimary();
-  if (!isPrimary) { fSource=3; return; }
+  if (!isPrimary) { fSource = 3; return; }
 
-  Int_t countP=-1, pdg=0;
+  Int_t countP = -1, pdg = 0;
   Int_t parents[10], parLine[10];
   AliAODMCParticle *mother = 0;
-  while(lineM>=0){
+  while (lineM>=0) {
     mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
     pdg = mother->GetPdgCode();
-    if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
+    if (pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
     parents[++countP] = pdg;
     parLine[countP] = lineM;
     lineM = mother->GetMother();
   }
-  for(Int_t i=0; i<=countP; i++){
+  for (Int_t i=0; i<=countP; i++) {
     fParentIndex[i] = parLine[countP-i];
     fParentPDGCode[i] = parents[countP-i];
   }
@@ -200,6 +213,7 @@ void AliMuonInfoStoreMC::SetMCInfoAOD(AliMCEvent *mcEvent, Int_t label)
   if (fIsFull && lineM>=0) this->FillHistoryQuarksAOD(mcEvent, lineM);
 
   fSource = this->SelectHFMuon();
+
   return;
 }
 
@@ -208,31 +222,35 @@ void AliMuonInfoStoreMC::SetMCInfoESD(AliMCEvent *mcEvent, Int_t label)
 {
   // fill track MC info with ESD base
   fTrackIndex = label;
-  if (fTrackIndex<0) { fSource=5; return; }
+  if (fTrackIndex<0) { fSource = 5; return; }
 
   TParticle *pMC = ((AliMCParticle*)mcEvent->GetTrack(fTrackIndex))->Particle();
   fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
 
   fTrackPDGCode = pMC->GetPdgCode();
-  if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
+  if (TMath::Abs(fTrackPDGCode)!=13) { fSource = 4; return; }
+
+  fGeneratorIndex = ((AliMCParticle*)mcEvent->GetTrack(fTrackIndex))->GetGeneratorIndex();
 
   Int_t lineM = pMC->GetFirstMother();
-  if (lineM<0) { fSource=2; return; }
+  Bool_t isPhysicalPrimary = mcEvent->Stack()->IsPhysicalPrimary(lineM);
+  if (isPhysicalPrimary) fIsPhyPrim = kTRUE;
+  if (lineM<0) { fSource = 2; return; } 
 
-  if (lineM>=mcEvent->Stack()->GetNprimary()) { fSource=3; return; }
+  if (lineM>=mcEvent->Stack()->GetNprimary()) { fSource = 3; return; }
 
-  Int_t countP=-1, pdg=0;
+  Int_t countP = -1, pdg = 0;
   Int_t parents[10], parLine[10];
   TParticle *mother = 0;
-  while(lineM>=0){
+  while (lineM>=0) {
     mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
     pdg = mother->GetPdgCode();
-    if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
+    if (pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
     parents[++countP] = pdg;
     parLine[countP] = lineM;
     lineM = mother->GetFirstMother();
   }
-  for(Int_t i=0; i<=countP; i++){
+  for (Int_t i=0; i<=countP; i++) {
     fParentIndex[i] = parLine[countP-i];
     fParentPDGCode[i] = parents[countP-i];
   }
@@ -241,6 +259,7 @@ void AliMuonInfoStoreMC::SetMCInfoESD(AliMCEvent *mcEvent, Int_t label)
   if (fIsFull && lineM>=0) this->FillHistoryQuarksESD(mcEvent, lineM);
 
   fSource = this->SelectHFMuon();
+
   return;
 }
 
@@ -250,9 +269,9 @@ void AliMuonInfoStoreMC::FillHistoryQuarksAOD(AliMCEvent* const mcEvent, Int_t l
   // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx 
 
   if (lineM<0) return;
-  Int_t countP=-1, pdg=0;
+  Int_t countP = -1, pdg = 0;
   AliAODMCParticle *mother = 0;
-  while(lineM>=0){
+  while (lineM>=0) {
     mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
     pdg = mother->GetPdgCode();
     fQuarkIndex[++countP] = lineM;
@@ -262,21 +281,21 @@ void AliMuonInfoStoreMC::FillHistoryQuarksAOD(AliMCEvent* const mcEvent, Int_t l
 
   // for PYTHIA checking
   countP = 1;
-  for(Int_t par=0; par<4; par++) {
-    if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
+  for (Int_t par=0; par<4; par++) {
+    if (TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
   }
-  if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
-    if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
+  if (this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
+    if (this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
                  this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
 
       pdg = this->QuarkPDGCode(countP);
       Int_t line = this->QuarkIndex(countP);
       this->ResetQuarkInfo();
-      while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
+      while (TMath::Abs(pdg)!=this->ParentFlavour(0)) {
         pdg = ((AliAODMCParticle*)mcEvent->GetTrack(++line))->GetPdgCode();
       }
-      while(line>=0){
+      while (line>=0) {
         mother = (AliAODMCParticle*)mcEvent->GetTrack(line);
         pdg = mother->GetPdgCode();
         fQuarkIndex[countP] = line;
@@ -285,6 +304,7 @@ void AliMuonInfoStoreMC::FillHistoryQuarksAOD(AliMCEvent* const mcEvent, Int_t l
       }
     }
   }
+
   return;
 }
 
@@ -294,9 +314,9 @@ void AliMuonInfoStoreMC::FillHistoryQuarksESD(AliMCEvent* const mcEvent, Int_t l
   // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx 
 
   if (lineM<0) return;
-  Int_t countP=-1, pdg=0;
+  Int_t countP = -1, pdg = 0;
   TParticle *mother = 0;
-  while(lineM>=0){
+  while (lineM>=0) {
     mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
     pdg = mother->GetPdgCode();
     fQuarkIndex[++countP] = lineM;
@@ -306,22 +326,22 @@ void AliMuonInfoStoreMC::FillHistoryQuarksESD(AliMCEvent* const mcEvent, Int_t l
 
   // for PYTHIA checking
   countP = 1;
-  for(Int_t par=0; par<4; par++) {
-    if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
+  for (Int_t par=0; par<4; par++) {
+    if (TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
   }
-  if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
-    if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
+  if (this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
+    if (this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
                  this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
 
       pdg = this->QuarkPDGCode(countP);
       Int_t line = this->QuarkIndex(countP);
       this->ResetQuarkInfo();
-      while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
-        pdg = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle()->GetPdgCode();
+      while (TMath::Abs(pdg)!=this->ParentFlavour(0)) {
+        pdg = ((AliMCParticle*)mcEvent->GetTrack(++line))->Particle()->GetPdgCode();
       }
       while(line>=0){
-        mother = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle();
+        mother = ((AliMCParticle*)mcEvent->GetTrack(++line))->Particle();
         pdg = mother->GetPdgCode();
         fQuarkIndex[countP] = line;
         fQuarkPDGCode[countP++] = pdg;
@@ -329,6 +349,7 @@ void AliMuonInfoStoreMC::FillHistoryQuarksESD(AliMCEvent* const mcEvent, Int_t l
       }
     }
   }
+
   return;
 }
 
@@ -343,7 +364,7 @@ Int_t AliMuonInfoStoreMC::SelectHFMuon()
   Bool_t isRes = kFALSE;
   Int_t i=0, nparents=this->ParentsN();
   while (i<nparents && !isRes) {
-    isRes = IsMotherAResonance(i++);
+    isRes = IsMotherAResonance(i++); // quarkonium mother has no anti-particle
   }
 
   if (isRes) return 2;
@@ -368,6 +389,7 @@ void AliMuonInfoStoreMC::ResetQuarkInfo()
     fQuarkIndex[pos] = -1;
     fQuarkPDGCode[pos] = 0;
   }
+
   return;
 }
 
@@ -378,6 +400,7 @@ Int_t AliMuonInfoStoreMC::ParentFlavour(Int_t i) const
   Int_t pdg = ParentPDGCode(i);
   pdg = TMath::Abs(pdg/100);
   if(pdg>9) pdg /= 10;
+
   return pdg;
 }
 
@@ -387,5 +410,6 @@ Bool_t AliMuonInfoStoreMC::IsMotherAResonance(Int_t i) const
   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
   Int_t pdg = ParentPDGCode(i);
   Int_t id=pdg%100000;
+
   return (!((id-id%10)%110));
 }
index 2a2192f..a5ad99e 100644 (file)
@@ -41,16 +41,18 @@ class AliMuonInfoStoreMC : public AliMuonInfoStoreRD {
   Bool_t IsMotherAResonance(Int_t i) const;
 
   TLorentzVector LorentzP()         const { return fLorentzP; }
-  Int_t    Source()                 const { return fSource; }
-  Int_t    TrackIndex()             const { return fTrackIndex; }
+  Int_t    Source()                 const { return fSource;   }
+  Int_t    TrackIndex()             const { return fTrackIndex;   }
   Int_t    TrackPDGCode()           const { return fTrackPDGCode; }
   Int_t    ParentsN()               const { return fNParents; }
-  Int_t    ParentIndex(Int_t i=0)   const { return (i<fNParents ? fParentIndex[i] : -1); }
+  Int_t    ParentIndex(Int_t i=0)   const { return (i<fNParents ? fParentIndex[i] : -1);  }
   Int_t    ParentPDGCode(Int_t i=0) const { return (i<fNParents ? fParentPDGCode[i] : 0); }
-  Int_t    QuarkIndex(Int_t i=0)    const { return (i<4 ? fQuarkIndex[i] : -1); }
+  Int_t    QuarkIndex(Int_t i=0)    const { return (i<4 ? fQuarkIndex[i] : -1);  }
   Int_t    QuarkPDGCode(Int_t i=0)  const { return (i<4 ? fQuarkPDGCode[i] : 0); }
   Bool_t   IsOscillation()          const { return fOscillation; }
-  Double_t Weight()                 const { return fWeight; }
+  Double_t Weight()                 const { return fWeight;    }
+  Bool_t   IsPhyPrimKPi()           const { return fIsPhyPrim; }
+  Short_t  GetGeneratorIndex()      const { return fGeneratorIndex; }
 
   static const char* StdBranchName() { return fgkStdBranchName.Data(); }
   static Int_t SourcesN()            { return fgkSourcesN;             }
@@ -66,31 +68,33 @@ class AliMuonInfoStoreMC : public AliMuonInfoStoreRD {
   Bool_t IsDiquark(Int_t pdg);
   void ResetQuarkInfo();
 
-  static const TString fgkStdBranchName;  // Standard branch name
-  static const Int_t   fgkSourcesN;       // num. of muon sources
-
-  Bool_t fIsFull;            // whether to use full mode (Pb-Pb)
-  TLorentzVector fLorentzP;  // lorentz momentum of particle
-  Int_t fTrackIndex;         // index of the MC particle
-  Int_t fTrackPDGCode;       // PDG code of the MC particle
-  Int_t fSource;  // = 0, mu<-b 
-                  // = 1, mu<-c 
-                  // = 2, primary mu
-                  // = 3, secondary mu
-                  // = 4, not mu
-                  // = 5, unidentified track
-
-  Int_t fParentIndex[5];    // index of parents
-  Int_t fParentPDGCode[5];  // PDG code of parents
-  Int_t fNParents;          // num. of parents
-  Bool_t fOscillation;      // flag of oscillation
-
-  Int_t fQuarkIndex[4];    // index of quarks
-  Int_t fQuarkPDGCode[4];  // PDG code of quarks
-
-  Double_t fWeight;  // for PbPb collisoions
-
-  ClassDef(AliMuonInfoStoreMC, 5);
+  static const TString fgkStdBranchName; // Standard branch name
+  static const Int_t   fgkSourcesN;      // num. of muon sources
+
+  Bool_t fIsFull;           // whether to use full mode (Pb-Pb)
+  TLorentzVector fLorentzP; // lorentz momentum of particle
+  Int_t fTrackIndex;        // index of the MC particle
+  Int_t fTrackPDGCode;      // PDG code of the MC particle
+  Int_t fSource; // = 0, mu<-b 
+                 // = 1, mu<-c 
+                 // = 2, primary mu
+                 // = 3, secondary mu
+                 // = 4, not mu
+                 // = 5, unidentified track
+
+  Int_t fParentIndex[5];   // index of parents
+  Int_t fParentPDGCode[5]; // PDG code of parents
+  Int_t fNParents;         // num. of parents
+  Bool_t fOscillation;     // flag of oscillation
+
+  Int_t fQuarkIndex[4];   // index of quarks
+  Int_t fQuarkPDGCode[4]; // PDG code of quarks
+
+  Double_t fWeight;         // for PbPb collisoions
+  Bool_t   fIsPhyPrim;      // select physical primary decay mu
+  Short_t  fGeneratorIndex; // index for the employed generator
+
+  ClassDef(AliMuonInfoStoreMC, 6);
 };
 
 #endif
index 7ae50d8..6634120 100644 (file)
@@ -156,6 +156,7 @@ void AliMuonInfoStoreRD::FillMuonInfo(AliAODTrack *trk)
   this->SetChi2FitMomentum(trk->Chi2perNDF());
   this->SetChi2MatchTrigger(trk->GetChi2MatchTrigger());
   this->SetRabsEnd(trk->GetRAtAbsorberEnd());
+
   return;
 }
 
@@ -175,5 +176,6 @@ void AliMuonInfoStoreRD::FillMuonInfo(AliESDMuonTrack *trk)
   this->SetChi2FitMomentum(trk->GetChi2()/(2.*trk->GetNHit()-5.));
   this->SetChi2MatchTrigger(trk->GetChi2MatchTrigger());
   this->SetRabsEnd(trk->GetRAtAbsorberEnd());
+
   return;
 }
index 0958b36..1940e09 100644 (file)
@@ -52,31 +52,31 @@ class AliMuonInfoStoreRD : public TObject {
   void FillMuonInfo(AliAODTrack *trk);
   void FillMuonInfo(AliESDMuonTrack *trk);
 
-  void SetMomentumAtVtx(Double_t p[3]) { fMomentumAtVtx.SetXYZ(p[0],p[1],p[2]); }
-  void SetMomentumAtDCA(Double_t p[3]) { fMomentumAtDCA.SetXYZ(p[0],p[1],p[2]); }
-  void SetMomentumUncor(Double_t p[3]) { fMomentumUncor.SetXYZ(p[0],p[1],p[2]); }
-  void SetDCA(Double_t dca[3])         { for (Int_t i=3; i--;) fDCA[i]=dca[i];  }
-  void SetCharge(Short_t charge)           { fCharge           = charge;  }
-  void SetMatchTrigger(Int_t trigger)      { fMatchTrigger     = trigger; }
-  void SetChi2FitMomentum(Double_t chi2)   { fChi2FitMomentum  = chi2;    }
-  void SetChi2MatchTrigger(Double_t chi2)  { fChi2MatchTrigger = chi2;    }
-  void SetRabsEnd(Double_t rAbsEnd)        { fRabsEnd          = rAbsEnd; }
-
-  static const TString fgkStdBranchName;  // Standard branch name
-
-  TVector3 fMomentumAtVtx;  // momentum corrected w vtx
-  TVector3 fMomentumAtDCA;  // momentum at DCA in vtx plane
-  TVector3 fMomentumUncor;  // momentum at first station
-
-  Double_t fDCA[3];            // distance of closet approach
-  Short_t  fCharge;            // track charge
-  Int_t    fMatchTrigger;      // type of match trigger
-  Double_t fChi2FitMomentum;   // chi2/NDF of momentum fit
-  Double_t fChi2MatchTrigger;  // chi2 of trigger matching
-  Double_t fRabsEnd;  // position at the end of front absorber
-  UInt_t   fSelMask;  // mask of single muon selection
-
-  ClassDef(AliMuonInfoStoreRD, 7);
+  void SetMomentumAtVtx(Double_t p[3])    { fMomentumAtVtx.SetXYZ(p[0],p[1],p[2]); }
+  void SetMomentumAtDCA(Double_t p[3])    { fMomentumAtDCA.SetXYZ(p[0],p[1],p[2]); }
+  void SetMomentumUncor(Double_t p[3])    { fMomentumUncor.SetXYZ(p[0],p[1],p[2]); }
+  void SetDCA(Double_t dca[3])            { for (Int_t i=3; i--;) fDCA[i]=dca[i];  }
+  void SetCharge(Short_t charge)          { fCharge           = charge;  }
+  void SetMatchTrigger(Int_t trigger)     { fMatchTrigger     = trigger; }
+  void SetChi2FitMomentum(Double_t chi2)  { fChi2FitMomentum  = chi2;    }
+  void SetChi2MatchTrigger(Double_t chi2) { fChi2MatchTrigger = chi2;    }
+  void SetRabsEnd(Double_t rAbsEnd)       { fRabsEnd          = rAbsEnd; }
+
+  static const TString fgkStdBranchName; // Standard branch name
+
+  TVector3 fMomentumAtVtx; // momentum corrected w vtx
+  TVector3 fMomentumAtDCA; // momentum at DCA in vtx plane
+  TVector3 fMomentumUncor; // momentum at first station
+
+  Double_t fDCA[3];           // distance of closet approach
+  Short_t  fCharge;           // track charge
+  Int_t    fMatchTrigger;     // type of match trigger
+  Double_t fChi2FitMomentum;  // chi2/NDF of momentum fit
+  Double_t fChi2MatchTrigger; // chi2 of trigger matching
+  Double_t fRabsEnd;          // position at the end of front absorber
+  UInt_t   fSelMask;          // mask of single muon selection
+
+  ClassDef(AliMuonInfoStoreRD, 8);
 };
 
 #endif
index 396a4cf..48b11f2 100644 (file)
@@ -36,7 +36,8 @@
 #include "AliVVertex.h"
 #include "AliCentrality.h"
 #include "AliEventplane.h"
-
+#include "AliAnalysisUtils.h"
+#include "AliMultiplicity.h"
 #include "AliMuonTrackCuts.h"
 #include "AliMuonInfoStoreRD.h"
 #include "AliMuonInfoStoreMC.h"
 #include "AliDimuInfoStoreMC.h"
 #include "AliMuonsHFHeader.h"
 
+#include "AliMCEvent.h"
+#include "AliGenEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
+
 ClassImp(AliMuonsHFHeader)
 
 const TString AliMuonsHFHeader::fgkStdBranchName("MuEvsH");
-Int_t         AliMuonsHFHeader::fgAnaMode = 0;
-Bool_t        AliMuonsHFHeader::fgIsMC    = kFALSE;
 Double_t      AliMuonsHFHeader::fgCuts[5] = { -999999., 999999., 999999., -999999., 999999. };
 
 //_____________________________________________________________________________
 AliMuonsHFHeader::AliMuonsHFHeader() :
 TNamed(),
+fAnaMode(0),
+fIsMC(kFALSE),
 fSelMask(AliVEvent::kAny),
 fIsMB(kFALSE),
 fIsMU(kFALSE),
-fIsPileupSPD(kFALSE),
 fVtxContrsN(0),
 fFiredTriggerClass(),
-fCentrality(-1.),
 fCentQA(-1),
-fEventPlane(0.)
+fEventPlane(0.),
+fIsEvtInChunk(kFALSE),
+fIsVtxSeled2013pA(kFALSE),
+fNumOfTrklets(-1),
+fTrgInpts(0),
+fImpParam(-1.),
+fCentralityV0M(-1.),
+fCentralityV0A(-1.),
+fCentralityV0C(-1.),
+fCentralityCL1(-1.),
+fCentralityZNA(-1.),
+fCentralityZNC(-1.),
+fPUMask(0)
 {
   //
   // default constructor
@@ -74,15 +89,27 @@ fEventPlane(0.)
 //_____________________________________________________________________________
 AliMuonsHFHeader::AliMuonsHFHeader(const AliMuonsHFHeader &src) :
 TNamed(),
+fAnaMode(src.fAnaMode),
+fIsMC(src.fIsMC),
 fSelMask(src.fSelMask),
 fIsMB(src.fIsMB),
 fIsMU(src.fIsMU),
-fIsPileupSPD(src.fIsPileupSPD),
 fVtxContrsN(src.fVtxContrsN),
 fFiredTriggerClass(src.fFiredTriggerClass),
-fCentrality(src.fCentrality),
 fCentQA(src.fCentQA),
-fEventPlane(src.fEventPlane)
+fEventPlane(src.fEventPlane),
+fIsEvtInChunk(src.fIsEvtInChunk),
+fIsVtxSeled2013pA(src.fIsVtxSeled2013pA),
+fNumOfTrklets(src.fNumOfTrklets),
+fTrgInpts(src.fTrgInpts),
+fImpParam(src.fImpParam),
+fCentralityV0M(src.fCentralityV0M),
+fCentralityV0A(src.fCentralityV0A),
+fCentralityV0C(src.fCentralityV0C),
+fCentralityCL1(src.fCentralityCL1),
+fCentralityZNA(src.fCentralityZNA),
+fCentralityZNC(src.fCentralityZNC),
+fPUMask(src.fPUMask)
 {
   //
   // copy constructor
@@ -100,15 +127,27 @@ AliMuonsHFHeader& AliMuonsHFHeader::operator=(const AliMuonsHFHeader &src)
 
   if(&src==this) return *this;
 
+  fAnaMode           = src.fAnaMode;
+  fIsMC              = src.fIsMC;
   fSelMask           = src.fSelMask;
   fIsMB              = src.fIsMB;
   fIsMU              = src.fIsMU;
-  fIsPileupSPD       = src.fIsPileupSPD;
   fVtxContrsN        = src.fVtxContrsN;
   fFiredTriggerClass = src.fFiredTriggerClass;
-  fCentrality        = src.fCentrality;
   fCentQA            = src.fCentQA;
   fEventPlane        = src.fEventPlane;
+  fIsEvtInChunk      = src.fIsEvtInChunk;
+  fIsVtxSeled2013pA  = src.fIsVtxSeled2013pA;
+  fNumOfTrklets      = src.fNumOfTrklets;
+  fTrgInpts          = src.fTrgInpts;
+  fImpParam          = src.fImpParam;
+  fCentralityV0M     = src.fCentralityV0M;
+  fCentralityV0A     = src.fCentralityV0A;
+  fCentralityV0C     = src.fCentralityV0C; 
+  fCentralityCL1     = src.fCentralityCL1; 
+  fCentralityZNA     = src.fCentralityZNA;
+  fCentralityZNC     = src.fCentralityZNC; 
+  fPUMask            = src.fPUMask;
   for (Int_t i=3; i--;) fVtx[i] = src.fVtx[i];
   for (Int_t i=3; i--;) fVMC[i] = src.fVMC[i];
 
@@ -124,39 +163,66 @@ AliMuonsHFHeader::~AliMuonsHFHeader()
 }
 
 //_____________________________________________________________________________
-void AliMuonsHFHeader::SetEventInfo(AliInputEventHandler* const handler, AliMCEvent* const eventMC)
+void AliMuonsHFHeader::SetEventInfo(AliInputEventHandler* const handler)
 {
   // fill info at event level
 
-  AliVEvent *event = handler->GetEvent();
-  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(event);
-  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
+  AliVEvent   *event = handler->GetEvent();
+  AliAODEvent *aod   = dynamic_cast<AliAODEvent*>(event);
+  AliESDEvent *esd   = dynamic_cast<AliESDEvent*>(event);
 
   fSelMask = handler->IsEventSelected();
-  if (aod) fFiredTriggerClass = aod->GetFiredTriggerClasses();
-  if (esd) fFiredTriggerClass = esd->GetFiredTriggerClasses();
+  if (aod) { fFiredTriggerClass = aod->GetFiredTriggerClasses(); fTrgInpts = aod->GetHeader()->GetL0TriggerInputs(); }
+  if (esd) { fFiredTriggerClass = esd->GetFiredTriggerClasses(); fTrgInpts = esd->GetHeader()->GetL0TriggerInputs(); }
   fIsMB = fSelMask & AliVEvent::kMB;
   fIsMU = fSelMask & AliVEvent::kMUON;
 
   const AliVVertex *vertex = event->GetPrimaryVertex();
   vertex->GetXYZ(fVtx);
   fVtxContrsN = vertex->GetNContributors();
-  if (fgIsMC) {
-    if (esd)   eventMC->GetPrimaryVertex()->GetXYZ(fVMC);
+  if (fIsMC) {
+    AliMCEvent *mcEvent = handler->MCEvent();
+    AliGenEventHeader *genHeader = mcEvent->GenEventHeader();
+    if (!genHeader) {  AliError("Header not found. Nothing done!"); return; }
+    AliGenDPMjetEventHeader *dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
+    if (dpmHeader) fImpParam = dpmHeader->ImpactParameter();
+
+    if (esd)   mcEvent->GetPrimaryVertex()->GetXYZ(fVMC);
     if (aod) ((AliAODMCHeader*)aod->FindListObject(AliAODMCHeader::StdBranchName()))->GetVertex(fVMC);
   } this->SetTitle(vertex->GetTitle());
-  fIsPileupSPD = (aod && !aod->GetTracklets()) ? event->IsPileupFromSPD(3,0.8,3.,2.,5.) : event->IsPileupFromSPDInMultBins();
-
-  AliCentrality *cent = event->GetCentrality();
-  if (cent) {
-    fCentrality = cent->GetCentralityPercentileUnchecked("V0M");
-    fCentQA     = cent->GetQuality();
+  //(aod && !aod->GetTracklets()) ? event->IsPileupFromSPD(3,0.8,3.,2.,5.) : event->IsPileupFromSPDInMultBins();
+
+  AliAnalysisUtils *anaUtils = new AliAnalysisUtils();
+  if (esd) {
+    fIsEvtInChunk     = anaUtils->IsFirstEventInChunk(esd);
+    fIsVtxSeled2013pA = anaUtils->IsVertexSelected2013pA(esd);
+    fNumOfTrklets     = esd->GetMultiplicity()->GetNumberOfTracklets();
+  } 
+  if (aod) {
+    fIsEvtInChunk     = anaUtils->IsFirstEventInChunk(aod);
+    fIsVtxSeled2013pA = anaUtils->IsVertexSelected2013pA(aod);
+    fNumOfTrklets     = aod->GetTracklets()->GetNumberOfTracklets();
   }
 
+  fPUMask = this->CollectPUMask(event);
+
   AliEventplane *evnP = event->GetEventplane();
   if (evnP) fEventPlane = evnP->GetEventplane("Q");
 //if (evnP) fEventPlane = evnP->GetEventplane("V0A");
 
+   AliCentrality *cent = event->GetCentrality(); if (cent) {
+    fCentQA        = cent->GetQuality();
+    fCentralityV0M = cent->GetCentralityPercentileUnchecked("V0M");
+    fCentralityV0A = cent->GetCentralityPercentileUnchecked("V0A");
+    fCentralityV0C = cent->GetCentralityPercentileUnchecked("V0C");
+    fCentralityCL1 = cent->GetCentralityPercentileUnchecked("CL1");
+    fCentralityZNA = cent->GetCentralityPercentileUnchecked("ZNA");
+    fCentralityZNC = cent->GetCentralityPercentileUnchecked("ZNC");
+  }
+
+  aod = 0; esd = 0;
+  delete anaUtils; anaUtils = 0;
+  
   return;
 }
 
@@ -164,13 +230,25 @@ void AliMuonsHFHeader::SetEventInfo(AliInputEventHandler* const handler, AliMCEv
 Bool_t AliMuonsHFHeader::IsSelected()
 {
   // select event according to the event selection cuts
-  if (this->VtxContrsN()<fgCuts[0])       return kFALSE;
-  if (TMath::Abs(this->Vz())>fgCuts[1])   return kFALSE;
-  if (this->Vt()>fgCuts[2])               return kFALSE;
+  if (this->VtxContrsN()<fgCuts[0])     return kFALSE;
+  if (TMath::Abs(this->Vz())>fgCuts[1]) return kFALSE;
+  if (this->Vt()>fgCuts[2])             return kFALSE;
 
   // centrality selection
-  Float_t centr = this->Centrality();
-  if (centr<fgCuts[3] || centr>fgCuts[4]) return kFALSE;
+  Float_t centrV0M = this->Centrality(kV0M);
+  Float_t centrV0A = this->Centrality(kV0A);
+  Float_t centrV0C = this->Centrality(kV0C); 
+  Float_t centrCL1 = this->Centrality(kCL1); 
+  Float_t centrZNA = this->Centrality(kZNA);
+  Float_t centrZNC = this->Centrality(kZNC); 
+
+  if (centrV0M<fgCuts[3] || centrV0M>fgCuts[4]) return kFALSE;
+  if (centrV0A<fgCuts[3] || centrV0A>fgCuts[4]) return kFALSE;
+  if (centrV0C<fgCuts[3] || centrV0C>fgCuts[4]) return kFALSE;
+  if (centrCL1<fgCuts[3] || centrCL1>fgCuts[4]) return kFALSE;
+  if (centrZNA<fgCuts[3] || centrZNA>fgCuts[4]) return kFALSE;
+  if (centrZNC<fgCuts[3] || centrZNC>fgCuts[4]) return kFALSE;
+
   return kTRUE;
 }
 
@@ -179,13 +257,13 @@ void AliMuonsHFHeader::CreateHistograms(TList *list)
 {
   // create output histos of muon analysis according to the analysis mode & MC flag
 
-  if (fgIsMC) {
+  if (fIsMC) {
     this->CreateHistosEvnH(list);
-    if (fgAnaMode!=2) {
+    if (fAnaMode!=2) {
       TString sName[7] = { "Unidentified", "Hadron", "SecondaryMu", "PrimaryMu", "CharmMu", "BottomMu", "" };
       for (Int_t i=7; i--;) this->CreateHistosMuon(list, sName[i]);
     }
-    if (fgAnaMode!=1) {
+    if (fAnaMode!=1) {
       TString sName[7] = { "Uncorr", "Resonance", "DDsame", "DDdiff", "BBsame", "BBdiff", "" };
       for (Int_t i=7; i--;) this->CreateHistosDimu(list, sName[i]);
     }
@@ -193,8 +271,9 @@ void AliMuonsHFHeader::CreateHistograms(TList *list)
   }
 
   this->CreateHistosEvnH(list,"MB"); this->CreateHistosEvnH(list,"MU");
-  if (fgAnaMode!=2) { this->CreateHistosMuon(list,"MB"); this->CreateHistosMuon(list,"MU"); }
-  if (fgAnaMode!=1) { this->CreateHistosDimu(list,"MB"); this->CreateHistosDimu(list,"MU"); }
+  if (fAnaMode!=2) { this->CreateHistosMuon(list,"MB"); this->CreateHistosMuon(list,"MU"); }
+  if (fAnaMode!=1) { this->CreateHistosDimu(list,"MB"); this->CreateHistosDimu(list,"MU"); }
+
   return;
 }
 
@@ -222,6 +301,7 @@ void AliMuonsHFHeader::CreateHistosEvnH(TList *list, TString sName)
   }
 
   TH1::AddDirectory(oldStatus);
+
   return;
 }
 
@@ -249,6 +329,7 @@ void AliMuonsHFHeader::CreateHistosMuon(TList *list, TString sName)
   }
 
   TH1::AddDirectory(oldStatus);
+
   return;
 }
 
@@ -278,6 +359,7 @@ void AliMuonsHFHeader::CreateHistosDimu(TList *list, TString sName)
   }
 
   TH1::AddDirectory(oldStatus);
+
   return;
 }
 
@@ -292,12 +374,13 @@ void AliMuonsHFHeader::FillHistosEvnH(TList *list)
   const Int_t nhs    = 3;
   TString tName[nhs] = {       "Vz",       "Vt",        "VtxNcontr" };
   Double_t dist[nhs] = { this->Vz(), this->Vt(), static_cast<Double_t>(this->VtxContrsN()) };
-  if (fgIsMC && (fSelMask & AliVEvent::kAny)) {
+  if (fIsMC && (fSelMask & AliVEvent::kAny)) {
     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h_%s",tName[i].Data())))->Fill(dist[i]);
   } else {
     if (fIsMB && (fSelMask & AliVEvent::kMB))   { for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s","MB",tName[i].Data())))->Fill(dist[i]); }
     if (fIsMU && (fSelMask & AliVEvent::kMUON)) { for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s","MU",tName[i].Data())))->Fill(dist[i]); }
   }
+
   return;
 }
 
@@ -320,7 +403,7 @@ void AliMuonsHFHeader::FillHistosMuon(TList *list, AliMuonInfoStoreRD* const inf
                          static_cast<Double_t>(infoStore->Charge()),
                          infoStore->RabsEnd() };
 
-  if (fgIsMC && (fSelMask & AliVEvent::kAny)) {
+  if (fIsMC && (fSelMask & AliVEvent::kAny)) {
     TString sName[7] = { "BottomMu", "CharmMu", "PrimaryMu", "SecondaryMu", "Hadron", "Unidentified", "" };
     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s",sName[6].Data(),tName[i].Data())))->Fill(dist[i]);
     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s",sName[s].Data(),tName[i].Data())))->Fill(dist[i]);
@@ -351,7 +434,7 @@ void AliMuonsHFHeader::FillHistosDimu(TList *list, AliDimuInfoStoreRD* const inf
                          infoStore->Momentum().Pt(),
                          infoStore->InvM() };
 
-  if (fgIsMC && (fSelMask & AliVEvent::kAny)) {
+  if (fIsMC && (fSelMask & AliVEvent::kAny)) {
     TString sName[7] = { "BBdiff", "BBsame", "DDdiff", "DDsame", "Resonance", "Uncorr", "" };
     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s_%s",sName[6].Data(),dimuName.Data(),tName[i].Data())))->Fill(dist[i]);
     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s_%s",sName[s].Data(),dimuName.Data(),tName[i].Data())))->Fill(dist[i]);
@@ -366,3 +449,70 @@ void AliMuonsHFHeader::FillHistosDimu(TList *list, AliDimuInfoStoreRD* const inf
 
   return;
 }
+
+//_____________________________________________________________________________
+Float_t AliMuonsHFHeader::Centrality(Int_t centrality)
+{
+  // obtain centrality via selected estimators
+
+  if (centrality==kV0M)      return fCentralityV0M;
+  else if (centrality==kV0A) return fCentralityV0A; 
+  else if (centrality==kV0C) return fCentralityV0C;
+  else if (centrality==kCL1) return fCentralityCL1;
+  else if (centrality==kZNA) return fCentralityZNA;
+  else if (centrality==kZNC) return fCentralityZNC;
+  else return -1.;
+}
+
+//_____________________________________________________________________________
+UInt_t AliMuonsHFHeader::CollectPUMask(AliVEvent *event)
+{
+  // collect the mask for different combination of the parameters;
+  // used to tag pile-up events (IsPileupFromSPD, no IsPileupFromSPDInMultBins at the moment) 
+
+  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(event);
+  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
+
+  UInt_t collectMask = 0;
+  Bool_t ISc1z1 = kFALSE, ISc1z2 = kFALSE, ISc1z3 = kFALSE, ISc1z4 = kFALSE,
+         ISc2z1 = kFALSE, ISc2z2 = kFALSE, ISc2z3 = kFALSE, ISc2z4 = kFALSE,
+         ISc3z1 = kFALSE, ISc3z2 = kFALSE, ISc3z3 = kFALSE, ISc3z4 = kFALSE,
+         ISc4z1 = kFALSE, ISc4z2 = kFALSE, ISc4z3 = kFALSE, ISc4z4 = kFALSE;
+
+  ISc1z1 = (aod) ? aod->IsPileupFromSPD(3,0.5,3.,2.,5.) : esd->IsPileupFromSPD(3,0.5,3.,2.,5.);
+  ISc1z2 = (aod) ? aod->IsPileupFromSPD(3,0.6,3.,2.,5.) : esd->IsPileupFromSPD(3,0.6,3.,2.,5.);
+  ISc1z3 = (aod) ? aod->IsPileupFromSPD(3,0.8,3.,2.,5.) : esd->IsPileupFromSPD(3,0.8,3.,2.,5.);
+  ISc1z4 = (aod) ? aod->IsPileupFromSPD(3,0.9,3.,2.,5.) : esd->IsPileupFromSPD(3,0.9,3.,2.,5.);
+  ISc2z1 = (aod) ? aod->IsPileupFromSPD(4,0.5,3.,2.,5.) : esd->IsPileupFromSPD(4,0.5,3.,2.,5.);
+  ISc2z2 = (aod) ? aod->IsPileupFromSPD(4,0.6,3.,2.,5.) : esd->IsPileupFromSPD(4,0.6,3.,2.,5.);
+  ISc2z3 = (aod) ? aod->IsPileupFromSPD(4,0.8,3.,2.,5.) : esd->IsPileupFromSPD(4,0.8,3.,2.,5.);
+  ISc2z4 = (aod) ? aod->IsPileupFromSPD(4,0.9,3.,2.,5.) : esd->IsPileupFromSPD(4,0.9,3.,2.,5.);
+  ISc3z1 = (aod) ? aod->IsPileupFromSPD(5,0.5,3.,2.,5.) : esd->IsPileupFromSPD(5,0.5,3.,2.,5.);
+  ISc3z2 = (aod) ? aod->IsPileupFromSPD(5,0.6,3.,2.,5.) : esd->IsPileupFromSPD(5,0.6,3.,2.,5.);
+  ISc3z3 = (aod) ? aod->IsPileupFromSPD(5,0.8,3.,2.,5.) : esd->IsPileupFromSPD(5,0.8,3.,2.,5.);
+  ISc3z4 = (aod) ? aod->IsPileupFromSPD(5,0.9,3.,2.,5.) : esd->IsPileupFromSPD(5,0.9,3.,2.,5.);
+  ISc4z1 = (aod) ? aod->IsPileupFromSPD(6,0.5,3.,2.,5.) : esd->IsPileupFromSPD(6,0.5,3.,2.,5.);
+  ISc4z2 = (aod) ? aod->IsPileupFromSPD(6,0.6,3.,2.,5.) : esd->IsPileupFromSPD(6,0.6,3.,2.,5.);
+  ISc4z3 = (aod) ? aod->IsPileupFromSPD(6,0.8,3.,2.,5.) : esd->IsPileupFromSPD(6,0.8,3.,2.,5.);
+  ISc4z4 = (aod) ? aod->IsPileupFromSPD(6,0.9,3.,2.,5.) : esd->IsPileupFromSPD(6,0.9,3.,2.,5.);
+  if (ISc1z1) collectMask |= kPUc1z1;
+  if (ISc1z2) collectMask |= kPUc1z2;
+  if (ISc1z3) collectMask |= kPUc1z3;
+  if (ISc1z4) collectMask |= kPUc1z4;
+  if (ISc2z1) collectMask |= kPUc2z1;
+  if (ISc2z2) collectMask |= kPUc2z2;
+  if (ISc2z3) collectMask |= kPUc2z3;
+  if (ISc2z4) collectMask |= kPUc2z4;
+  if (ISc3z1) collectMask |= kPUc3z1;
+  if (ISc3z2) collectMask |= kPUc3z2;
+  if (ISc3z3) collectMask |= kPUc3z3;
+  if (ISc3z4) collectMask |= kPUc3z4;
+  if (ISc4z1) collectMask |= kPUc4z1;
+  if (ISc4z2) collectMask |= kPUc4z2;
+  if (ISc4z3) collectMask |= kPUc4z3;
+  if (ISc4z4) collectMask |= kPUc4z4;
+  aod = 0; esd = 0;
+
+  return collectMask;
+}
index d365096..ff411e4 100644 (file)
@@ -26,31 +26,29 @@ class AliMuonInfoStoreMC;
 class AliDimuInfoStoreMC;
 
 class AliMuonsHFHeader : public TNamed {
- public :
+ public:
 
   AliMuonsHFHeader();
   AliMuonsHFHeader(const AliMuonsHFHeader &src);
   AliMuonsHFHeader& operator=(const AliMuonsHFHeader &src);
   ~AliMuonsHFHeader();
 
-  void GetVMC(Double_t *vtx)  const { for (Int_t i=3; i--;) vtx[i]=fVMC[i]; }
-  void GetXYZ(Double_t *vtx)  const { for (Int_t i=3; i--;) vtx[i]=fVtx[i]; }
-  Double_t Vx()               const { return fVtx[0]; }
-  Double_t Vy()               const { return fVtx[1]; }
-  Double_t Vz()               const { return fVtx[2]; }
-  Double_t Vt()               const { return TMath::Sqrt(fVtx[0]*fVtx[0] + fVtx[1]*fVtx[1]); }
-  Int_t VtxContrsN()          const { return fVtxContrsN; }
-  TString FiredTriggerClass() const { return fFiredTriggerClass; }
-  UInt_t SelectionMask()      const { return fSelMask; }
-  Bool_t IsMB()               const { return fIsMB; }
-  Bool_t IsMU()               const { return fIsMU; }
-  Bool_t IsPileupSPD()        const { return fIsPileupSPD; }
-  Float_t    Centrality()     const { return fCentrality; }
-  Int_t      CentQA()         const { return fCentQA;}
-  Double32_t EventPlane()     const { return fEventPlane; }
-  Bool_t IsSelected();
-
-  void SetEventInfo(AliInputEventHandler* const handler, AliMCEvent* const eventMC);
+  void       GetVMC(Double_t *vtx) const { for (Int_t i=3; i--;) vtx[i]=fVMC[i]; }
+  void       GetXYZ(Double_t *vtx) const { for (Int_t i=3; i--;) vtx[i]=fVtx[i]; }
+  Double_t   Vx()                  const { return fVtx[0]; }
+  Double_t   Vy()                  const { return fVtx[1]; }
+  Double_t   Vz()                  const { return fVtx[2]; }
+  Double_t   Vt()                  const { return TMath::Sqrt(fVtx[0]*fVtx[0] + fVtx[1]*fVtx[1]); }
+  Int_t      VtxContrsN()          const { return fVtxContrsN; }
+  TString    FiredTriggerClass()   const { return fFiredTriggerClass; }
+  UInt_t     SelectionMask()       const { return fSelMask; }
+  Bool_t     IsMB()                const { return fIsMB; }
+  Bool_t     IsMU()                const { return fIsMU; }
+  Int_t      CentQA()              const { return fCentQA;}
+  Double32_t EventPlane()          const { return fEventPlane; }
+  Bool_t     IsSelected();
+
+  void SetEventInfo(AliInputEventHandler* const handler);
 
   void CreateHistograms(TList *list);
   void FillHistosEvnH(TList *list);
@@ -58,40 +56,84 @@ class AliMuonsHFHeader : public TNamed {
   void FillHistosDimu(TList *list, AliDimuInfoStoreRD* infoStore, Int_t src=0);
 
   static const char* StdBranchName()             { return fgkStdBranchName.Data();          }
-  static void SetAnaMode(Int_t anaMode=0)        { fgAnaMode=anaMode;                       }
-  static void SetIsMC(Int_t isMC=kFALSE)         { fgIsMC   =isMC;                          }
   static void SetSelectionCuts(Double_t cuts[5]) { for (Int_t i=5; i--;) fgCuts[i]=cuts[i]; }
 
+  void    SetAnaMode(Int_t anaMode)       { fAnaMode = anaMode;       }
+  void    SetIsMC(Int_t isMC)             { fIsMC    = isMC;          }
+  Int_t   GetAnaMode()              const { return fAnaMode;          }
+  Int_t   GetNumOfTracklets()       const { return fNumOfTrklets;     }
+  UInt_t  GetTrgInpts()             const { return fTrgInpts;         }
+  Float_t GetImpParam()             const { return fImpParam;         }
+
+  Bool_t  IsMC()                    const { return fIsMC;             }
+  Bool_t  IsEvtInChunk()            const { return fIsEvtInChunk;     }
+  Bool_t  IsVtxSeled2013pA()        const { return fIsVtxSeled2013pA; }
+
+  enum { kV0M=0, kV0A, kV0C, kCL1, kZNA, kZNC };
+  Float_t Centrality(Int_t centEstor=0);
+
+  enum {
+    kPUc1z1 = BIT(0),  // (minContributor=3, minZdis=0.5)
+    kPUc1z2 = BIT(1),  // (3, 0.6)
+    kPUc1z3 = BIT(2),  // (3, 0.8)
+    kPUc1z4 = BIT(3),  // (3, 0.9)
+    kPUc2z1 = BIT(4),  // (minContributor=4, minZdis=0.5)
+    kPUc2z2 = BIT(5),  // (4, 0.6)
+    kPUc2z3 = BIT(6),  // (4, 0.8)
+    kPUc2z4 = BIT(7),  // (4, 0.9)
+    kPUc3z1 = BIT(8),  // (minContributor=5, minZdis=0.5)
+    kPUc3z2 = BIT(9),  // (5, 0.6)
+    kPUc3z3 = BIT(10), // (5, 0.8)
+    kPUc3z4 = BIT(11), // (5, 0.9)
+    kPUc4z1 = BIT(12), // (minContributor=6, minZdis=0.5)
+    kPUc4z2 = BIT(13), // (6, 0.6)
+    kPUc4z3 = BIT(14), // (6, 0.8)
+    kPUc4z4 = BIT(15)  // (6, 0.9)
+  };
+  UInt_t GetPUMask() const { return fPUMask; }
+
  private :
 
   void CreateHistosEvnH(TList *list, TString sName="");
   void CreateHistosMuon(TList *list, TString sName="");
   void CreateHistosDimu(TList *list, TString sName="");
 
-  static const TString fgkStdBranchName;  // Standard branch name
-  static Int_t  fgAnaMode;                // analysis mode
-  static Bool_t fgIsMC;                   // flag to use MC
-  static Double_t fgCuts[5];  // 0, low limit of num. of vtx contributors
-                              // 1, up limit of vz
-                              // 2, up limit of vt
-                              // 3, centrality max
-                              // 4, centrality min
-
-  UInt_t fSelMask;     // mask of physics selection
-  Bool_t fIsMB;        // is min. bias triggered event (for real data)
-  Bool_t fIsMU;        // is MUON triggered event (for real data)
-  Bool_t fIsPileupSPD; // is pileup from SPD
-  Double_t fVtx[3];    // position of vtx
-  Double_t fVMC[3];    // position of vtx in MC
-  Int_t fVtxContrsN;   // num. of contributors of vtx rec
-
-  TString fFiredTriggerClass; // trigger class
-
-  Float_t fCentrality;  // event centrality class
-  Int_t   fCentQA;      // quality of centrality determination
-  Double32_t fEventPlane; // event plane angle
-
-  ClassDef(AliMuonsHFHeader, 7)
+  UInt_t CollectPUMask(AliVEvent *event); // collect the mask for the plie-up identification 
+
+  static const TString fgkStdBranchName; // Standard branch name
+  static Double_t fgCuts[5]; // 0, low limit of num. of vtx contributors
+                             // 1, up limit of vz
+                             // 2, up limit of vt
+                             // 3, centrality max
+                             // 4, centrality min
+
+  Int_t      fAnaMode;           // analysis mode
+  Bool_t     fIsMC;              // flag to use MC
+  UInt_t     fSelMask;           // mask of physics selection
+  Bool_t     fIsMB;              // is min. bias triggered event (for real data)
+  Bool_t     fIsMU;              // is MUON triggered event (for real data)
+  Double_t   fVtx[3];            // position of vtx
+  Double_t   fVMC[3];            // position of vtx in MC
+  Int_t      fVtxContrsN;        // num. of contributors of vtx rec
+  TString    fFiredTriggerClass; // trigger class
+  Int_t      fCentQA;            // quality of centrality determination
+  Double32_t fEventPlane;        // event plane angle
+  Bool_t     fIsEvtInChunk;      // flag denotes whether the event is in the chunk
+  Bool_t     fIsVtxSeled2013pA;  // vertex selection based on 2013pA criteria
+  Int_t      fNumOfTrklets;      // num. of tracklets
+  UInt_t     fTrgInpts;          // trigger inputs (ID)
+  Float_t    fImpParam;          // impact parameter
+  
+  Float_t fCentralityV0M; //centrality determinated via employed estimators
+  Float_t fCentralityV0A;
+  Float_t fCentralityV0C;
+  Float_t fCentralityCL1;
+  Float_t fCentralityZNA;
+  Float_t fCentralityZNC;
+  
+  UInt_t fPUMask; // mask for tagging pile-up event
+
+  ClassDef(AliMuonsHFHeader, 8)
 };
 
 #endif