]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Analysis tasks for photon analysis (by Sudipan)
authorcoppedis <coppedis@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 28 May 2013 12:48:25 +0000 (12:48 +0000)
committercoppedis <coppedis@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 28 May 2013 12:48:25 +0000 (12:48 +0000)
PWGLF/FORWARD/photons/AliAnalysisTaskPMD.cxx [new file with mode: 0644]
PWGLF/FORWARD/photons/AliAnalysisTaskPMD.h [new file with mode: 0644]
PWGLF/FORWARD/photons/AliAnalysisTaskPMDSim.cxx [new file with mode: 0644]
PWGLF/FORWARD/photons/AliAnalysisTaskPMDSim.h [new file with mode: 0644]

diff --git a/PWGLF/FORWARD/photons/AliAnalysisTaskPMD.cxx b/PWGLF/FORWARD/photons/AliAnalysisTaskPMD.cxx
new file mode 100644 (file)
index 0000000..1bd77ff
--- /dev/null
@@ -0,0 +1,201 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDPmdTrack.h"
+#include "AliESDVertex.h"
+#include "AliESDVZERO.h"
+#include "AliAnalysisTaskPMD.h"
+
+// AnalysisTask For PMD
+// Authors: Sudipan De, Subhash Singha
+
+ClassImp(AliAnalysisTaskPMD)
+
+//________________________________________________________________________
+AliAnalysisTaskPMD::AliAnalysisTaskPMD(const char *name) 
+: AliAnalysisTaskSE(name), 
+  fESD(0), 
+  fOutputList(0), 
+  fHistTotEvent(0),
+  fHistTotEventAfterPhySel(0),
+  fHistTotEventAfterVtx(0),
+  fHistVtxZ(0),
+  fHistXYPre(0),
+  fHistEta(0),
+  fHistEta1(0),
+  fHistMultMeas(0),
+  fHistMultMeas1(0)
+{
+  for(Int_t i=0; i<10; i++){
+    fHistMultMeasEtaBinA[i] = 0;
+    fHistMultMeasEtaBinA1[i] = 0;
+  }
+  // Constructor
+  
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #0 id reserved by the base class for AOD
+  // Output slot #1 writes into a TH1 container
+  DefineOutput(1, TList::Class());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPMD::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+  
+  fOutputList = new TList();
+  Int_t kNbinsMultA = 50;  Float_t XminMultA = -0.5; Float_t XmaxMultA = 49.5;
+  Int_t kNbinsXY = 200; Float_t XminXY = -100.0; Float_t XmaxXY  = 100.0;
+  Int_t kNBinsEvent = 10; Float_t XminEvent = 0; Float_t XmaxEvent = 10;
+  Int_t kNBinsEta = 10; Float_t XminEta = 2.1; Float_t XmaxEta = 4.1;
+  fHistTotEvent = new TH1F("TotEvent","TotEvent",
+                          kNBinsEvent,XminEvent,XmaxEvent); 
+  fOutputList->Add(fHistTotEvent);
+  fHistTotEventAfterPhySel = new TH1F("TotEventAfterPhySel","TotEventAfterPhySel",
+                                     kNBinsEvent,XminEvent,XmaxEvent); 
+  fOutputList->Add(fHistTotEventAfterPhySel);
+  fHistTotEventAfterVtx = new TH1F("TotEventAfterVtx","TotEventAfterVtx",
+                                  kNBinsEvent,XminEvent,XmaxEvent); 
+  fOutputList->Add(fHistTotEventAfterVtx);
+  fHistVtxZ = new TH1F("VtxZ","VtxZ",100,-50,50);
+  fOutputList->Add(fHistVtxZ);
+  fHistXYPre = new TH2F("XYPre","XYPre",kNbinsXY,XminXY,XmaxXY,
+                       kNbinsXY,XminXY,XmaxXY);
+  fOutputList->Add(fHistXYPre);
+  fHistEta = new TH1F ("Eta","Eta",kNBinsEta,XminEta,XmaxEta);
+  fOutputList->Add(fHistEta);
+  fHistEta1 = new TH1F ("Eta1","Eta1",kNBinsEta,XminEta,XmaxEta);
+  fOutputList->Add(fHistEta1);
+  fHistMultMeas = new TH1F("MultM","MultM",100,-0.5,99.5);
+  fOutputList->Add(fHistMultMeas);
+  fHistMultMeas1 = new TH1F("MultM1","MultM1",100,-0.5,99.5);
+  fOutputList->Add(fHistMultMeas1);
+  
+  Char_t nameM[256], nameM1[256];
+  for(Int_t i=0; i<10; i++){
+    sprintf(nameM,"MultM_EtaBin%d",i+1);
+    fHistMultMeasEtaBinA[i] =
+      new TH1F(nameM,nameM,kNbinsMultA,
+              XminMultA,XmaxMultA);
+    fOutputList->Add(fHistMultMeasEtaBinA[i]);
+    sprintf(nameM1,"MultM1_EtaBin%d",i+1);
+    fHistMultMeasEtaBinA1[i] =
+      new TH1F(nameM1,nameM1,kNbinsMultA,
+              XminMultA,XmaxMultA);
+    fOutputList->Add(fHistMultMeasEtaBinA1[i]);
+  }//i loop
+  
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPMD::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+  Float_t MipCut1 = 432;//MPV=72=> 6*MPV=432
+  Float_t MipCut2 = 648;//MPV=72=> 9*MPV=648
+  Float_t etacls, theta, rdist;
+  Int_t PhotonCls = 0;//# of photon measured within 2.3 to 3.9 #eta
+  Int_t PhotonCls1 = 0;//# of photon measured within 2.3 to 3.9 #eta
+  Int_t PhotonClsAEtaBin[10] = {0};//# of photon measured with diff Eta bin
+  Int_t PhotonClsAEtaBin1[10] = {0};//# of photon measured with diff Eta bin
+  // Post output data.
+  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+  if (!fESD) {
+    printf("ERROR: fESD not available\n");
+    return;
+  }
+  fHistTotEvent->Fill(5);
+  //event selection
+  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
+  if (! isSelected) return;
+  fHistTotEventAfterPhySel->Fill(5);
+  //printf("There are %d tracks in this event\n", fESD->GetNumberOfPmdTracks());
+  //vertex selection
+  const AliESDVertex *vertex = fESD->GetPrimaryVertex();
+  Float_t Vz = vertex->GetZv();    
+  Bool_t zVerStatus = vertex->GetStatus();
+  if(zVerStatus){
+    fHistVtxZ->Fill(Vz);
+    if(TMath::Abs(Vz)<10){
+      fHistTotEventAfterVtx->Fill(5);
+      Int_t ptracks = fESD->GetNumberOfPmdTracks();
+      //PMDtrack Loop
+      for(Int_t kk=0;kk<ptracks;kk++){
+       AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(kk);
+       Int_t   det   = pmdtr->GetDetector();
+       Float_t clsX  = pmdtr->GetClusterX();
+       Float_t clsY  = pmdtr->GetClusterY();
+       Float_t clsZ  = pmdtr->GetClusterZ();
+       clsZ -= Vz;
+       Float_t ncell = pmdtr->GetClusterCells();
+       Float_t adc   = pmdtr->GetClusterADC();
+       //Int_t smn = pmdtr->GetSmn();
+       //Calculation of #eta
+       rdist = TMath::Sqrt(clsX*clsX + clsY*clsY);
+       if(clsZ!=0) theta = TMath::ATan2(rdist,clsZ);
+       etacls  = -TMath::Log(TMath::Tan(0.5*theta));
+       
+       if(det==0 && adc>72)fHistXYPre->Fill(clsX,clsY);
+       if(det==0 && adc>MipCut1 && ncell>2){
+         if(etacls > 2.1 && etacls<= 2.3)PhotonClsAEtaBin[0]++;
+         if(etacls > 2.3 && etacls<= 2.5)PhotonClsAEtaBin[1]++;
+         if(etacls > 2.5 && etacls<= 2.7)PhotonClsAEtaBin[2]++;
+         if(etacls > 2.7 && etacls<= 2.9)PhotonClsAEtaBin[3]++;
+         if(etacls > 2.9 && etacls<= 3.1)PhotonClsAEtaBin[4]++;
+         if(etacls > 3.1 && etacls<= 3.3)PhotonClsAEtaBin[5]++;
+         if(etacls > 3.3 && etacls<= 3.5)PhotonClsAEtaBin[6]++;
+         if(etacls > 3.5 && etacls<= 3.7)PhotonClsAEtaBin[7]++;
+         if(etacls > 3.7 && etacls<= 3.9)PhotonClsAEtaBin[8]++;
+         if(etacls > 3.9 && etacls<= 4.1)PhotonClsAEtaBin[9]++;
+         if(etacls > 2.3 && etacls<= 3.9)PhotonCls++;
+         fHistEta->Fill(etacls);
+       }//if MipCut1
+       if(det==0 && adc>MipCut2 && ncell>2){
+         if(etacls > 2.1 && etacls<= 2.3)PhotonClsAEtaBin1[0]++;
+         if(etacls > 2.3 && etacls<= 2.5)PhotonClsAEtaBin1[1]++;
+         if(etacls > 2.5 && etacls<= 2.7)PhotonClsAEtaBin1[2]++;
+         if(etacls > 2.7 && etacls<= 2.9)PhotonClsAEtaBin1[3]++;
+         if(etacls > 2.9 && etacls<= 3.1)PhotonClsAEtaBin1[4]++;
+         if(etacls > 3.1 && etacls<= 3.3)PhotonClsAEtaBin1[5]++;
+         if(etacls > 3.3 && etacls<= 3.5)PhotonClsAEtaBin1[6]++;
+         if(etacls > 3.5 && etacls<= 3.7)PhotonClsAEtaBin1[7]++;
+         if(etacls > 3.7 && etacls<= 3.9)PhotonClsAEtaBin1[8]++;
+         if(etacls > 3.9 && etacls<= 4.1)PhotonClsAEtaBin1[9]++;
+         if(etacls > 2.3 && etacls<= 3.9)PhotonCls1++;
+         fHistEta1->Fill(etacls);
+       }//if Mipcut2
+      } //PMDtrack loop 
+      for(Int_t i=0; i<10; i++){
+       fHistMultMeasEtaBinA[i]->Fill(PhotonClsAEtaBin[i]);
+       fHistMultMeasEtaBinA1[i]->Fill(PhotonClsAEtaBin1[i]);
+      }//i loop
+      fHistMultMeas->Fill(PhotonCls);
+      fHistMultMeas1->Fill(PhotonCls1);
+    }// Vz loop 
+  }// Vzcut loop
+  PostData(1, fOutputList);
+}      
+
+//________________________________________________________________________
+void AliAnalysisTaskPMD::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+
+  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutputList) {
+    printf("ERROR: Output list not available\n");
+    return;
+  }
+}
diff --git a/PWGLF/FORWARD/photons/AliAnalysisTaskPMD.h b/PWGLF/FORWARD/photons/AliAnalysisTaskPMD.h
new file mode 100644 (file)
index 0000000..b219d63
--- /dev/null
@@ -0,0 +1,64 @@
+#ifndef AliAnalysisTaskPMD_cxx\r
+#define AliAnalysisTaskPMD_cxx\r
+\r
+// AnalysisTask For PMD \r
+// Authors: Sudipan De, Subhash Singha\r
+\r
+class TH1F;\r
+class TH2F;\r
+class AliESDEvent;\r
+class AliESDPmdTrack;\r
+class AliESDVertex;\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliAnalysisTaskPMD : public AliAnalysisTaskSE {\r
+ public:\r
+ AliAnalysisTaskPMD() : AliAnalysisTaskSE(), \r
+    fESD(0), \r
+    fOutputList(0), \r
+    fHistTotEvent(0),\r
+    fHistTotEventAfterPhySel(0),\r
+    fHistTotEventAfterVtx(0),\r
+    fHistVtxZ(0),\r
+    fHistXYPre(0),\r
+    fHistEta(0),\r
+    fHistEta1(0),\r
+    fHistMultMeas(0),\r
+    fHistMultMeas1(0)\r
+      {\r
+       for(Int_t i=0; i<10; i++){\r
+         fHistMultMeasEtaBinA[i] = 0;\r
+         fHistMultMeasEtaBinA1[i] = 0;\r
+       }\r
+      }\r
+  AliAnalysisTaskPMD(const char *name);\r
+  virtual ~AliAnalysisTaskPMD() {}\r
+  \r
+  virtual void   UserCreateOutputObjects();\r
+  virtual void   UserExec(Option_t *option);\r
+  virtual void   Terminate(Option_t *);\r
+  \r
+ private:\r
+  AliESDEvent *fESD;    //! ESD object\r
+  TList *fOutputList; //! Output list\r
+  TH1F *fHistTotEvent; //total event\r
+  TH1F *fHistTotEventAfterPhySel; //# event after physics selection\r
+  TH1F *fHistTotEventAfterVtx; //# event after vertex cut\r
+  TH1F *fHistVtxZ;//z vertex distribution\r
+  TH2F *fHistXYPre;//2d scatter plot pre\r
+  TH1F *fHistEta; // eta distribution in PMD coverage\r
+  TH1F *fHistEta1; // eta distribution in PMD coverage\r
+  TH1F *fHistMultMeas;//measured multiplicity (2.3-3.9)\r
+  TH1F *fHistMultMeas1;//measured multiplicity (2.3-3.9)\r
+  \r
+  TH1F *fHistMultMeasEtaBinA[10];//meas. mult. dist. for diff. eta bins\r
+  TH1F *fHistMultMeasEtaBinA1[10];//meas. mult. dist. for diff. eta bins\r
+  \r
+  AliAnalysisTaskPMD(const AliAnalysisTaskPMD&); // not implemented\r
+  AliAnalysisTaskPMD& operator=(const AliAnalysisTaskPMD&); // not implemented\r
+  \r
+  ClassDef(AliAnalysisTaskPMD, 1); // example of analysis\r
+};\r
+\r
+#endif\r
diff --git a/PWGLF/FORWARD/photons/AliAnalysisTaskPMDSim.cxx b/PWGLF/FORWARD/photons/AliAnalysisTaskPMDSim.cxx
new file mode 100644 (file)
index 0000000..04e9cc0
--- /dev/null
@@ -0,0 +1,295 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDPmdTrack.h"
+#include "AliESDVertex.h"
+#include "AliAnalysisTaskPMDSim.h"
+#include <TParticle.h>
+#include <AliLog.h>
+#include <AliStack.h>
+#include <AliHeader.h>
+#include <AliGenEventHeader.h>
+#include <AliHeader.h>
+#include <AliGenPythiaEventHeader.h>
+#include <AliMCEventHandler.h>
+#include <AliMCEvent.h>
+#include "AliTriggerAnalysis.h"
+#include "AliPhysicsSelection.h"
+
+// AnalysisTask For PMD
+// Authors: Sudipan De, Subhash Singha
+
+ClassImp(AliAnalysisTaskPMDSim)
+
+//________________________________________________________________________
+AliAnalysisTaskPMDSim::AliAnalysisTaskPMDSim(const char *name) 
+  : AliAnalysisTaskSE(name), 
+  fESD(0), 
+  fOutputList(0), 
+  fHistTotEvent(0),
+  fHistTotEventAfterPhySel(0),
+  fHistTotEventAfterVtx(0),
+  fVtxZ(0),
+  fHistXYPre(0),
+  fHistEtaPhM(0),
+  fHistEtaPhM1(0),
+  fHistEtaT(0),
+  fMultMeasured(0),
+  fMultMeasured1(0),
+  fMultTrue(0),
+  fMultCorr(0),
+  fMultCorr1(0)
+{
+  for(Int_t i=0; i<10; i++){
+    fHistMultMeasEtaBinA[i] = 0;
+    fHistMultMeasEtaBinA1[i] = 0;
+    fHistMultTrueEtaBinA[i] = 0;
+    fHistMultCorrEtaBinA[i] = 0;
+    fHistMultCorrEtaBinA1[i] = 0;
+  }
+  
+  // Constructor
+  
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #0 id reserved by the base class for AOD
+  // Output slot #1 writes into a TH1 container
+  DefineOutput(1, TList::Class());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPMDSim::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+  
+  fOutputList = new TList();
+  
+  Int_t kNbinsMult = 100;  Float_t XminMult = -0.5; 
+  Float_t XmaxMult = 99.5;//mult
+  Int_t kNbinsMultA = 50;  Float_t XminMultA = -0.5; 
+  Float_t XmaxMultA = 49.5;//mult
+  Int_t kNbinEtaPMDCov = 10; Float_t XminEtaPMDCov = 2.1; Float_t XmaxEtaPMDCov = 4.1;//etaPMD
+  Int_t kNBinsEvent = 10; Float_t XminEvent = 0; Float_t XmaxEvent = 10;
+  Int_t kNbinsXY = 200; Float_t XminXY = -100.0; Float_t XmaxXY  = 100.0;
+  fHistTotEvent = new TH1F("TotEvent","TotEvent",
+                          kNBinsEvent,XminEvent,XmaxEvent); 
+  fOutputList->Add(fHistTotEvent);
+  fHistTotEventAfterPhySel = new TH1F("TotEventAfterPhySel","TotEventAfterPhySel",
+                                     kNBinsEvent,XminEvent,XmaxEvent); 
+  fOutputList->Add(fHistTotEventAfterPhySel);
+  fHistTotEventAfterVtx = new TH1F("TotEventAfterVtx","TotEventAfterVtx",
+                                 kNBinsEvent,XminEvent,XmaxEvent); 
+  fOutputList->Add(fHistTotEventAfterVtx);
+  fVtxZ     = new TH1F("VtxZ","VtxZ",200,-50.0,50.0);
+  fOutputList->Add(fVtxZ);
+  fHistXYPre = new TH2F("XYPre","XYPre",kNbinsXY,XminXY,XmaxXY,
+                        kNbinsXY,XminXY,XmaxXY);
+  fOutputList->Add(fHistXYPre);
+  fHistEtaPhM = new TH1F("fHistEtaPhM","fHistEtaPhM",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
+  fOutputList->Add(fHistEtaPhM);
+  fHistEtaPhM1 = new TH1F("fHistEtaPhM1","fHistEtaPhM1",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
+  fOutputList->Add(fHistEtaPhM1);
+  fHistEtaT = new TH1F("fHistEtaT","fHistEtaT",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
+  fOutputList->Add(fHistEtaT);
+  fMultMeasured = new TH1F("MultM","MultM",kNbinsMult,XminMult,XmaxMult);
+  fOutputList->Add(fMultMeasured);
+  fMultMeasured1 = new TH1F("MultM1","MultM1",kNbinsMult,XminMult,XmaxMult);
+  fOutputList->Add(fMultMeasured1);
+  fMultTrue = new TH1F("MultT","MultT",kNbinsMult,XminMult,XmaxMult);
+  fOutputList->Add(fMultTrue);
+  fMultCorr = new TH2F("MultCorr","MultCorr",kNbinsMult,XminMult,XmaxMult,
+                      kNbinsMult,XminMult,XmaxMult);
+  fOutputList->Add(fMultCorr);
+  fMultCorr1 = new TH2F("MultCorr1","MultCorr1",kNbinsMult,XminMult,XmaxMult,
+                      kNbinsMult,XminMult,XmaxMult);
+  fOutputList->Add(fMultCorr1);
+   
+  
+  Char_t nameT[256], nameM[256], nameCorr[256],nameM1[256], nameCorr1[256];
+  for(Int_t i=0; i<10; i++){
+    sprintf(nameM,"MultM_EtaBin%d",i+1);
+    sprintf(nameM1,"MultM1_EtaBin%d",i+1);
+    sprintf(nameT,"MultT_EtaBin%d",i+1);
+    sprintf(nameCorr,"MultCorr_EtaBin%d",i+1);
+    sprintf(nameCorr1,"MultCorr1_EtaBin%d",i+1);
+    fHistMultMeasEtaBinA[i] = 
+      new TH1F(nameM,nameM,kNbinsMultA,XminMultA,XmaxMultA);
+    fHistMultMeasEtaBinA1[i] = 
+      new TH1F(nameM1,nameM1,kNbinsMultA,XminMultA,XmaxMultA);
+    fHistMultTrueEtaBinA[i] = 
+      new TH1F(nameT,nameT,kNbinsMultA,XminMultA,XmaxMultA);
+    fHistMultCorrEtaBinA[i] = 
+      new TH2F(nameCorr,nameCorr,kNbinsMultA,XminMultA,XmaxMultA,
+              kNbinsMultA,XminMultA,XmaxMultA);
+    fHistMultCorrEtaBinA1[i] = 
+      new TH2F(nameCorr1,nameCorr1,kNbinsMultA,XminMultA,XmaxMultA,
+              kNbinsMultA,XminMultA,XmaxMultA);
+    fOutputList->Add(fHistMultMeasEtaBinA[i]);
+    fOutputList->Add(fHistMultMeasEtaBinA1[i]);
+    fOutputList->Add(fHistMultTrueEtaBinA[i]);
+    fOutputList->Add(fHistMultCorrEtaBinA[i]);
+    fOutputList->Add(fHistMultCorrEtaBinA1[i]);
+  }//i loop
+  //fOutputList->Add(fHistPt);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskPMDSim::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+  Float_t MipCut1 = 432;//MPV=432=> 6*MPV=432
+  Float_t MipCut2 = 648;//MPV=432=> 6*MPV=648 
+  Float_t etacls, theta, rdist;
+  Int_t PhotonCls = 0;
+  Int_t PhotonCls1 = 0;
+  Int_t PhotonClsAEtaBin[10] = {0};//#of photon measured with diff Eta bin
+  Int_t PhotonClsAEtaBin1[10] = {0};//#of photon measured with diff Eta bin
+  Int_t PhotonTrueAEtaBin[10] = {0};//#of photon Incident with diff Eta bin
+  // Post output data.
+  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+  if (!fESD) {
+    printf("ERROR: fESD not available\n");
+    return;
+  }
+  fHistTotEvent->Fill(5);
+  //event selection
+  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
+  if (! isSelected) return;
+  fHistTotEventAfterPhySel->Fill(5);
+  //Vertex selection
+  const AliESDVertex *vertex = fESD->GetPrimaryVertex();
+  Float_t Vz = vertex->GetZv();    
+  // Float_t Vx = vertex->GetXv();    
+  //Float_t Vy = vertex->GetYv();
+  Bool_t zVerStatus = vertex->GetStatus();
+  if(zVerStatus){
+    fVtxZ->Fill(Vz);
+    if(TMath::Abs(Vz)<10){   
+      fHistTotEventAfterVtx->Fill(5);   
+      Int_t ptracks = fESD->GetNumberOfPmdTracks();
+      for(Int_t kk=0;kk<ptracks;kk++){
+       AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(kk);
+       Int_t   det   = pmdtr->GetDetector();
+       Float_t clsX  = pmdtr->GetClusterX();
+       Float_t clsY  = pmdtr->GetClusterY();
+       Float_t clsZ  = pmdtr->GetClusterZ();
+       clsZ -= Vz;
+       Float_t ncell = pmdtr->GetClusterCells();
+       Float_t adc   = pmdtr->GetClusterADC();
+       //Float_t pid   = pmdtr->GetClusterPID();
+       //Float_t isotag = pmdtr->GetClusterSigmaX();
+       //Int_t trno = pmdtr->GetClusterTrackNo();
+       //calculation of #eta
+       rdist = TMath::Sqrt(clsX*clsX + clsY*clsY);
+       if(clsZ!=0) theta = TMath::ATan2(rdist,clsZ);
+       etacls   = -TMath::Log(TMath::Tan(0.5*theta));
+       
+       if(det==0 && adc>0)fHistXYPre->Fill(clsX,clsY);
+       if(det==0 && adc>MipCut1 && ncell>2){
+         fHistEtaPhM->Fill(etacls);
+         if(etacls > 2.3 && etacls <= 3.9)  PhotonCls++;
+         if(etacls > 2.1 && etacls <= 2.3)  PhotonClsAEtaBin[0]++;
+         if(etacls > 2.3 && etacls <= 2.5)  PhotonClsAEtaBin[1]++;
+         if(etacls > 2.5 && etacls <= 2.7)  PhotonClsAEtaBin[2]++;
+         if(etacls > 2.7 && etacls <= 2.9)  PhotonClsAEtaBin[3]++;
+         if(etacls > 2.9 && etacls <= 3.1)  PhotonClsAEtaBin[4]++;
+         if(etacls > 3.1 && etacls <= 3.3)  PhotonClsAEtaBin[5]++;
+         if(etacls > 3.3 && etacls <= 3.5)  PhotonClsAEtaBin[6]++;
+         if(etacls > 3.5 && etacls <= 3.7)  PhotonClsAEtaBin[7]++;
+         if(etacls > 3.7 && etacls <= 3.9)  PhotonClsAEtaBin[8]++;
+         if(etacls > 3.9 && etacls <= 4.1)  PhotonClsAEtaBin[9]++;
+       }//if MipCut1
+       if( det==0 && adc>MipCut2 && ncell>2){
+         fHistEtaPhM1->Fill(etacls);
+         if(etacls > 2.3 && etacls <= 3.9)  PhotonCls1++;
+         if(etacls > 2.1 && etacls <= 2.3)  PhotonClsAEtaBin1[0]++;
+         if(etacls > 2.3 && etacls <= 2.5)  PhotonClsAEtaBin1[1]++;
+         if(etacls > 2.5 && etacls <= 2.7)  PhotonClsAEtaBin1[2]++;
+         if(etacls > 2.7 && etacls <= 2.9)  PhotonClsAEtaBin1[3]++;
+         if(etacls > 2.9 && etacls <= 3.1)  PhotonClsAEtaBin1[4]++;
+         if(etacls > 3.1 && etacls <= 3.3)  PhotonClsAEtaBin1[5]++;
+         if(etacls > 3.3 && etacls <= 3.5)  PhotonClsAEtaBin1[6]++;
+         if(etacls > 3.5 && etacls <= 3.7)  PhotonClsAEtaBin1[7]++;
+         if(etacls > 3.7 && etacls <= 3.9)  PhotonClsAEtaBin1[8]++;
+         if(etacls > 3.9 && etacls <= 4.1)  PhotonClsAEtaBin1[9]++;
+       }//if MipCut2
+      } //track loop 
+      //reading MC info.
+      AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*>
+       (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+      if (!eventHandler) {
+       Printf("ERROR: Could not retrieve MC event handler");
+       return;
+      }
+      AliMCEvent* mcEvent = eventHandler->MCEvent();
+      if (!mcEvent) {
+       Printf("ERROR: Could not retrieve MC event");
+       return;
+      }
+      AliStack* stack = mcEvent->Stack();
+      if (!stack)
+       {
+         AliDebug(AliLog::kError, "Stack not available");
+         return;
+       }
+      if(stack){
+       Int_t nPrim  = stack->GetNprimary();
+       Int_t PhotonTrue = 0;
+       for (Int_t iMc = 0; iMc < nPrim; ++iMc)
+         {
+           TParticle *MPart = stack->Particle(iMc);
+           Int_t mpart  = MPart->GetPdgCode();
+           Float_t eta   = MPart->Eta();
+           if(mpart == 22){  
+             fHistEtaT->Fill(eta);
+             if(eta > 2.3 && eta <= 3.9)PhotonTrue++;
+             if(eta > 2.1 && eta <= 2.3)PhotonTrueAEtaBin[0]++;
+             if(eta > 2.3 && eta <= 2.5)PhotonTrueAEtaBin[1]++;
+             if(eta > 2.5 && eta <= 2.7)PhotonTrueAEtaBin[2]++;
+             if(eta > 2.7 && eta <= 2.9)PhotonTrueAEtaBin[3]++;
+             if(eta > 2.9 && eta <= 3.1)PhotonTrueAEtaBin[4]++;
+             if(eta > 3.1 && eta <= 3.3)PhotonTrueAEtaBin[5]++;
+             if(eta > 3.3 && eta <= 3.5)PhotonTrueAEtaBin[6]++;
+             if(eta > 3.5 && eta <= 3.7)PhotonTrueAEtaBin[7]++;
+             if(eta > 3.7 && eta <= 3.9)PhotonTrueAEtaBin[8]++;
+             if(eta > 3.9 && eta <= 4.1)PhotonTrueAEtaBin[9]++;
+           }//mpart
+         }//track loop
+       for(Int_t i=0; i<10; i++){
+         fHistMultMeasEtaBinA[i]->Fill(PhotonClsAEtaBin[i]);
+         fHistMultMeasEtaBinA1[i]->Fill(PhotonClsAEtaBin1[i]);
+         fHistMultTrueEtaBinA[i]->Fill(PhotonTrueAEtaBin[i]);
+         fHistMultCorrEtaBinA[i]->Fill(PhotonTrueAEtaBin[i],PhotonClsAEtaBin[i]);
+         fHistMultCorrEtaBinA1[i]->Fill(PhotonTrueAEtaBin[i],PhotonClsAEtaBin1[i]);
+       }//i loop
+       fMultMeasured->Fill(PhotonCls);
+       fMultMeasured1->Fill(PhotonCls1);
+       fMultTrue->Fill(PhotonTrue);
+       fMultCorr->Fill(PhotonTrue,PhotonCls);
+       fMultCorr1->Fill(PhotonTrue,PhotonCls1);
+      }//if stack
+    }//if Vz!= 0
+  }//Vz loop
+  PostData(1, fOutputList);
+}//UserExec()     
+
+//_______________________________________________________________________
+void AliAnalysisTaskPMDSim::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutputList) {
+    printf("ERROR: Output list not available\n");
+    return;
+  }
+}//Terminate
diff --git a/PWGLF/FORWARD/photons/AliAnalysisTaskPMDSim.h b/PWGLF/FORWARD/photons/AliAnalysisTaskPMDSim.h
new file mode 100644 (file)
index 0000000..12e748e
--- /dev/null
@@ -0,0 +1,81 @@
+#ifndef AliAnalysisTaskPMDSim_cxx\r
+#define AliAnalysisTaskPMDSim_cxx\r
+\r
+// AnalysisTask For PMD\r
+// Authors: Sudipan De, Subhash Singha\r
+\r
+class TH1F;\r
+class TH2F;\r
+class AliESDEvent;\r
+class AliESDPmdTrack;\r
+class AliESDVertex;\r
+class AliStack;\r
+class AliHeader;\r
+class AliGenEventHeader;\r
+class TParticle;\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliAnalysisTaskPMDSim : public AliAnalysisTaskSE {\r
+ public:\r
+ AliAnalysisTaskPMDSim() : AliAnalysisTaskSE(), \r
+    fESD(0), \r
+    fOutputList(0), \r
+    fHistTotEvent(0),\r
+    fHistTotEventAfterPhySel(0),\r
+    fHistTotEventAfterVtx(0),\r
+    fVtxZ(0),\r
+    fHistXYPre(0),\r
+    fHistEtaPhM(0),\r
+    fHistEtaPhM1(0),\r
+    fHistEtaT(0),\r
+    fMultMeasured(0),\r
+    fMultMeasured1(0),\r
+    fMultTrue(0),\r
+    fMultCorr(0),\r
+    fMultCorr1(0) {\r
+    for(Int_t i=0; i<10; i++){\r
+      fHistMultMeasEtaBinA[i] = 0;\r
+      fHistMultMeasEtaBinA1[i] = 0;\r
+      fHistMultTrueEtaBinA[i] = 0;\r
+      fHistMultCorrEtaBinA[i] = 0;\r
+      fHistMultCorrEtaBinA1[i] = 0;\r
+    }\r
+  }\r
+  AliAnalysisTaskPMDSim(const char *name);\r
+  virtual ~AliAnalysisTaskPMDSim() {}\r
+  \r
+  virtual void   UserCreateOutputObjects();\r
+  virtual void   UserExec(Option_t *option);\r
+  virtual void   Terminate(Option_t *);\r
\r
+ private:\r
+  AliESDEvent *fESD;    //! ESD object\r
+  TList *fOutputList; //! Output list\r
+  TH1F *fHistTotEvent; //total event\r
+  TH1F *fHistTotEventAfterPhySel; //total event after physel\r
+  TH1F *fHistTotEventAfterVtx; //# event after vertex cut\r
+  TH1F *fVtxZ;//Vertex Z\r
+  TH2F *fHistXYPre;//2d scatter plot pre\r
+  TH1F *fHistEtaPhM;\r
+  TH1F *fHistEtaPhM1;\r
+  TH1F *fHistEtaT;\r
+  TH1F *fMultMeasured;\r
+  TH1F *fMultMeasured1;\r
+  TH1F *fMultTrue;\r
+  TH2F *fMultCorr;\r
+  TH2F *fMultCorr1;\r
+  \r
+  TH2F *fHistMultCorrEtaBinA[10];//mult. corr. for diff. eta bin\r
+  TH2F *fHistMultCorrEtaBinA1[10];//mult. corr. for diff. eta bin\r
+  TH1F *fHistMultTrueEtaBinA[10];//multTrue\r
+  TH1F *fHistMultMeasEtaBinA[10];//meas. mult. dist. for diff. eta bins\r
+  TH1F *fHistMultMeasEtaBinA1[10];//meas. mult. dist. for diff. eta bins\r
+  \r
+  AliAnalysisTaskPMDSim(const AliAnalysisTaskPMDSim&); // not implemented\r
+  AliAnalysisTaskPMDSim& operator=(const AliAnalysisTaskPMDSim&); // not implemented\r
+  \r
+  ClassDef(AliAnalysisTaskPMDSim, 1); // example of analysis\r
+};\r
+\r
+#endif\r