fix of mc particles indexing
authormcosenti <mcosenti@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 3 May 2013 18:55:13 +0000 (18:55 +0000)
committermcosenti <mcosenti@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 3 May 2013 18:55:13 +0000 (18:55 +0000)
including daughter information and status code in mc branch
adding isolation at mc level
fixing trigger bit selection

PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx
PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhoton.cxx
PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhoton.h

index f4d1e90..b573cb5 100644 (file)
@@ -234,6 +234,9 @@ void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
   const Int_t ndims =   fNDimensions;
   Double_t xmin[] = { 0.,   0.,  -10.,   -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,-1.5};
   Double_t xmax[] = { 100., 4., 190., 190., 190.,  190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,99.5};
+  if(fPeriod.Contains("11h")){
+    xmax[12]=3999.5;
+  }
   fHnOutput =  new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, CeIsoNoUESub, AllIsoNoUESub, d#phi_{trk},d#eta_{trk},#eta_{clus},#phi_{clus},T_{max},mult,mc-p_{T}^{#gamma}", ndims, bins, xmin, xmax);
   fOutputList->Add(fHnOutput);
 
@@ -264,6 +267,12 @@ void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
       else
        isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
   }
+  if(fPeriod.Contains("11h")){
+    if(fTrigBit.Contains("kAny"))
+      isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);
+    if(fTrigBit.Contains("kAnyINT"))
+      isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAnyINT);
+  }
   if(fIsMc)
     isSelected = kTRUE;
   if(fDebug)
index ebaa5dc..662a95c 100644 (file)
@@ -325,8 +325,10 @@ void AliAnalysisTaskEMCALPhoton::UserExec(Option_t *)
   fHeader->fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5); 
 
   fMCEvent = MCEvent();
-  if(fMCEvent)
+  if(fMCEvent){
     fStack = (AliStack*)fMCEvent->Stack();
+    fHeader->fNMcParts = fStack->GetNtrack();
+  }
 
   
   FindConversions();
@@ -679,45 +681,59 @@ void  AliAnalysisTaskEMCALPhoton::GetMcParts()
   const AliVVertex *evtVtx = fMCEvent->GetPrimaryVertex();
   if (!evtVtx)
     return;
-  Int_t ipart = 0;
   Int_t nTracks = fStack->GetNtrack();
+  AliPhotonMcPartObj *mcp = 0;
   for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
     TParticle *mcP = static_cast<TParticle*>(fStack->Particle(iTrack));
-    if (!mcP)
+    if (!mcP){
+      mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
       continue;
+    }
     // primary particle
     Double_t dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) + 
                               (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
                               (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
-    if(dR > 0.5)
+    if(dR > 0.5){
+      mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
       continue;
+    }
     
     // kinematic cuts
     Double_t pt = mcP->Pt() ;
-    if (pt<0.5)
+    if (pt<0.5){
+      mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
       continue;
+    }
     Double_t eta = mcP->Eta();
-    if (TMath::Abs(eta)>0.7)
+    if (TMath::Abs(eta)>0.7){
+      mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
       continue;
+    }
     Double_t phi  = mcP->Phi();
-    if (phi<1.0||phi>3.3)
+    if (phi<1.0||phi>3.3){
+      mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
       continue;
+    }
     // pion or eta meson or direct photon
     if(mcP->GetPdgCode() == 111) {
     } else if(mcP->GetPdgCode() == 221) {
     } else if(mcP->GetPdgCode() == 22 ) {
-    } else
+    } else {
+      mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
       continue;
+    }
+    
     bool checkIfAlreadySaved = false;
     for(Int_t imy=0;imy<fMyMcParts->GetEntries();imy++){
       AliPhotonMcPartObj *mymc = static_cast<AliPhotonMcPartObj*>(fMyMcParts->At(imy));
       if(!mymc)
        continue;
-      if(mymc->fLabel == iTrack)
+      if(imy == iTrack)
        checkIfAlreadySaved = true;
     }
-    if(!checkIfAlreadySaved)
-      FillMcPart(mcP, ipart++, iTrack);
+    if(!checkIfAlreadySaved){
+      FillMcPart(mcP,  iTrack);
+    }
     for(Int_t id=mcP->GetFirstDaughter(); id <= mcP->GetLastDaughter(); id++){
       if(id<=mcP->GetMother(0))
        continue;
@@ -726,20 +742,20 @@ void  AliAnalysisTaskEMCALPhoton::GetMcParts()
       TParticle *mcD = static_cast<TParticle*>(fStack->Particle(id));
       if(!mcD)
        continue;
-      FillMcPart(mcD, ipart++, id);
-    }
+      FillMcPart(mcD, id);
+      }
   }
 }
 
 //________________________________________________________________________
-void  AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP, Int_t ipart, Int_t itrack)
+void  AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP,  Int_t itrack)
 {
   // Fill MC particles.
 
   if(!mcP)
     return;
   TVector3 vmcv(mcP->Vx(),mcP->Vy(), mcP->Vz());
-  AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(ipart));
+  AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(itrack));
   mcp->fLabel = itrack ;
   mcp->fPdg = mcP->GetPdgCode() ;
   mcp->fPt = mcP->Pt() ;
@@ -751,6 +767,46 @@ void  AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP, Int_t ipart, Int_t
     mcp->fVPhi = vmcv.Phi() ;
   }
   mcp->fMother = mcP->GetMother(0) ;
+  mcp->fFirstD = mcP->GetFirstDaughter() ;
+  mcp->fLastD = mcP->GetLastDaughter() ;
+  mcp->fStatus = mcP->GetStatusCode();
+  mcp->fIso = AliAnalysisTaskEMCALPhoton::GetMcIsolation(mcP, itrack, 0.4 , 0.2);
+}
+//________________________________________________________________________                                                                                                                                   
+Double_t AliAnalysisTaskEMCALPhoton::GetMcIsolation(TParticle *mcP, Int_t itrack, Double_t radius, Double_t pt) const
+{
+  if (!fStack)
+    return -1;
+  if(!mcP)
+    return -1;
+  if(itrack<6 || itrack>8)
+    return -1;
+  if(mcP->GetPdgCode()!=22)
+    return -1;
+  Double_t sumpt=0;
+  Float_t eta = mcP->Eta();
+  Float_t phi = mcP->Phi();
+  Float_t dR;
+  Int_t nparts =  fStack->GetNtrack();
+  for(Int_t ip = 0; ip<nparts; ip++){
+    TParticle *mcisop = static_cast<TParticle*>(fStack->Particle(ip));
+    if(!mcisop)
+      continue;
+    if(ip==itrack)
+      continue;
+    if(mcisop->GetStatusCode()!=1)
+      continue;
+    if(mcisop->Pt()<pt)
+      continue;
+    TVector3 vmcv(mcisop->Vx(),mcisop->Vy(), mcisop->Vz());  
+    if(vmcv.Perp()>1)
+      continue;
+    dR = TMath::Sqrt((phi-mcisop->Phi())*(phi-mcisop->Phi())+(eta-mcisop->Eta())*(eta-mcisop->Eta()));
+    if(dR>radius)
+      continue;
+    sumpt += mcisop->Pt();
+  }
+  return sumpt;
 }
 
 //________________________________________________________________________
index 5a1a939..ce308db 100644 (file)
@@ -50,8 +50,9 @@ class AliAnalysisTaskEMCALPhoton : public AliAnalysisTaskSE {
   void         FillMyClusters();\r
   void         FillMyAltClusters();\r
   void         FillIsoTracks();\r
-  void         FillMcPart(TParticle *mcP, Int_t ipart, Int_t itrack);\r
+  void         FillMcPart(TParticle *mcP,  Int_t itrack);\r
   void         GetMcParts();\r
+  Double_t     GetMcIsolation(TParticle *mcP, Int_t itrack, Double_t radius, Double_t pt)                 const;\r
   Double_t     GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius=0.2, Double_t pt=0.)       const;\r
   Double_t     GetPhiBandEt(Double_t cEta, Double_t cPhi, Double_t radius=0.2, Double_t pt=0.)            const;\r
  // Double_t     GetPhiBandEt(Double_t cEta, Double_t cPhi, Double_t radius=0.2, Double_t pt=0.)            const;\r
@@ -121,7 +122,7 @@ class AliPhotonHeaderObj : public TObject
 {\r
   public: AliPhotonHeaderObj() :\r
   TObject(), fInputFileName(""), fTrClassMask(0), fTrCluster(0), fV0Cent(0), fV0(0), fCl1Cent(0), \r
-         fCl1(0), fTrCent(0), fTr(0), fNClus(0), fNCells(0), fTrackMult(0)  {;}\r
+    fCl1(0), fTrCent(0), fTr(0), fNClus(0), fNCells(0), fTrackMult(0), fNMcParts(0)  {;}\r
   public:\r
   TString       fInputFileName;  // used for normalization purposes in MC productions\r
   ULong64_t     fTrClassMask;    //         trigger class mask\r
@@ -135,8 +136,9 @@ class AliPhotonHeaderObj : public TObject
   Int_t         fNClus;\r
   Int_t         fNCells;\r
   Int_t         fTrackMult;\r
+  Int_t         fNMcParts;\r
 \r
-  ClassDef(AliPhotonHeaderObj,4)\r
+  ClassDef(AliPhotonHeaderObj,5)\r
 };\r
 \r
 class AliPhotonConvObj : public TObject\r
@@ -251,8 +253,9 @@ class AliPhotonTrackObj : public TObject
 class AliPhotonMcPartObj : public TObject\r
 {\r
   public: AliPhotonMcPartObj() :\r
-        TObject(), fLabel(-1), fPdg(0), fPt(0), fEta(0), fPhi(0), \r
-        fVR(0), fVEta(0), fVPhi(0), fMother(-1)  {;}\r
+  TObject(), fLabel(-1), fPdg(0), fPt(0), fEta(0), fPhi(0), \r
+    fVR(0), fVEta(0), fVPhi(0), fMother(-1), fFirstD(-1),\r
+    fLastD(-1), fStatus(-1), fIso(-1) {;}\r
   public:\r
   Short_t    fLabel;\r
   Short_t    fPdg;\r
@@ -263,8 +266,12 @@ class AliPhotonMcPartObj : public TObject
   Double32_t fVEta;\r
   Double32_t fVPhi;\r
   Short_t    fMother;\r
+  Short_t    fFirstD;\r
+  Short_t    fLastD;\r
+  Short_t    fStatus;\r
+  Double32_t fIso;\r
 \r
-  ClassDef(AliPhotonMcPartObj,1)\r
+  ClassDef(AliPhotonMcPartObj,2)\r
 };\r
 \r
 #endif\r