Gamma conversion analysis moved to its own subdirectory from PartCorr
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 Nov 2008 17:53:32 +0000 (17:53 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 Nov 2008 17:53:32 +0000 (17:53 +0000)
PWG4/CMake_libPWG4GammaConv.txt [new file with mode: 0644]
PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx [new file with mode: 0644]
PWG4/GammaConv/AliAnalysisTaskGammaConversion.h [new file with mode: 0644]
PWG4/GammaConv/AliGammaConversionHistograms.cxx [new file with mode: 0644]
PWG4/GammaConv/AliGammaConversionHistograms.h [new file with mode: 0644]
PWG4/GammaConv/AliV0Reader.cxx [new file with mode: 0644]
PWG4/GammaConv/AliV0Reader.h [new file with mode: 0644]
PWG4/PWG4GammaConvLinkDef.h [new file with mode: 0644]
PWG4/libPWG4GammaConv.pkg [new file with mode: 0644]

diff --git a/PWG4/CMake_libPWG4GammaConv.txt b/PWG4/CMake_libPWG4GammaConv.txt
new file mode 100644 (file)
index 0000000..f5ad0f2
--- /dev/null
@@ -0,0 +1,12 @@
+# -*- mode: cmake -*-
+
+set(SRCS
+       GammaConv/AliV0Reader.cxx GammaConv/AliAnalysisTaskGammaConversion.cxx GammaConv/AliGammaConversionHistograms.cxx
+)
+
+# fill list of header files from list of source files
+# by exchanging the file extension
+String(REPLACE ".cxx" ".h" HDRS "${SRCS}")
+
+AddLibrary(PWG4PartCorr "${SRCS}" "${HDRS}")
+
diff --git a/PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx b/PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx
new file mode 100644 (file)
index 0000000..f6edf1f
--- /dev/null
@@ -0,0 +1,618 @@
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt                        *\r
+ * Version 1.0                                                            *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+////////////////////////////////////////////////\r
+//--------------------------------------------- \r
+// Class used to do analysis on conversion pairs\r
+//---------------------------------------------\r
+////////////////////////////////////////////////\r
+\r
+// root\r
+#include <TChain.h>\r
+\r
+// analysis\r
+#include "AliAnalysisTaskGammaConversion.h"\r
+#include "AliStack.h"\r
+#include "AliLog.h"\r
+#include <vector>\r
+\r
+class AliKFVertex;\r
+class AliAODHandler;\r
+class AliAODEvent;\r
+class ALiESDEvent;\r
+class AliMCEvent;\r
+class AliMCEventHandler;\r
+class AliESDInputHandler;\r
+class AliAnalysisManager;\r
+class Riostream;\r
+class TFile;\r
+class TInterpreter;\r
+class TSystem;\r
+class TROOT;\r
+\r
+ClassImp(AliAnalysisTaskGammaConversion)\r
+\r
+\r
+AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():\r
+  AliAnalysisTaskSE(),\r
+  fV0Reader(NULL),\r
+  fStack(NULL),\r
+  fOutputContainer(NULL),\r
+  fHistograms(NULL),\r
+  fDoMCTruth(kFALSE),\r
+  fMCAllGammas(),\r
+  fMCPi0s(),\r
+  fMCEtas(),\r
+  fMCGammaChic(),\r
+  fKFReconstructedGammas(),\r
+  fElectronMass(-1),\r
+  fGammaMass(-1),\r
+  fPi0Mass(-1),\r
+  fEtaMass(-1),\r
+  fGammaWidth(-1),\r
+  fPi0Width(-1),\r
+  fEtaWidth(-1),\r
+  fCalculateBackground(kFALSE)\r
+{\r
+  // Default constructor\r
+  // Common I/O in slot 0\r
+  DefineInput (0, TChain::Class());\r
+  DefineOutput(0, TTree::Class());\r
+\r
+  // Your private output\r
+  DefineOutput(1, TList::Class());\r
+}\r
+\r
+AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
+  AliAnalysisTaskSE(name),\r
+  fV0Reader(NULL),\r
+  fStack(NULL),\r
+  fOutputContainer(0x0),\r
+  fHistograms(NULL),\r
+  fDoMCTruth(kFALSE),\r
+  fMCAllGammas(),\r
+  fMCPi0s(),\r
+  fMCEtas(),\r
+  fMCGammaChic(),\r
+  fKFReconstructedGammas(),\r
+  fElectronMass(-1),\r
+  fGammaMass(-1),\r
+  fPi0Mass(-1),\r
+  fEtaMass(-1),\r
+  fGammaWidth(-1),\r
+  fPi0Width(-1),\r
+  fEtaWidth(-1),\r
+  fCalculateBackground(kFALSE)\r
+{\r
+  // Common I/O in slot 0\r
+  DefineInput (0, TChain::Class());\r
+  DefineOutput(0, TTree::Class());\r
+  \r
+  // Your private output\r
+  DefineOutput(1, TList::Class());\r
+}\r
+\r
+AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
+{\r
+  // Remove all pointers\r
\r
+  if(fOutputContainer){\r
+    fOutputContainer->Clear() ; \r
+    delete fOutputContainer ;\r
+  }\r
+  if(fHistograms){\r
+    delete fHistograms;\r
+  }\r
+  if(fV0Reader){\r
+    delete fV0Reader;\r
+  }\r
+}\r
+\r
+\r
+void AliAnalysisTaskGammaConversion::Init()\r
+{\r
+  // Initialization\r
+  AliLog::SetGlobalLogLevel(AliLog::kError);\r
+}\r
+\r
+\r
+void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
+{\r
+  // Execute analysis for current event\r
+  \r
+  ConnectInputData("");\r
+  \r
+  //clear vectors\r
+  fMCAllGammas.clear();\r
+  fMCPi0s.clear();\r
+  fMCEtas.clear();\r
+  fMCGammaChic.clear();\r
+\r
+  fKFReconstructedGammas.clear();\r
+\r
+  //Clear the data in the v0Reader\r
+  fV0Reader->UpdateEventByEventData();\r
+\r
+  // Process the MC information\r
+  if(fDoMCTruth){\r
+    ProcessMCData();\r
+  }\r
+\r
+  // Process the v0 information\r
+  ProcessV0s();\r
+\r
+  //calculate background if flag is set\r
+  if(fCalculateBackground){\r
+    CalculateBackground();\r
+  }\r
+\r
+  // Process reconstructed gammas\r
+  ProcessGammasForNeutralMesonAnalysis();\r
+\r
+  PostData(1, fOutputContainer);\r
+  \r
+}\r
+\r
+void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
+  // see header file for documentation\r
+\r
+  if(fV0Reader == NULL){\r
+    // Write warning here cuts and so on are default if this ever happens\r
+  }\r
+  fV0Reader->Initialize();\r
+}\r
+\r
+void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
+  // see header file for documentation\r
+  \r
+  fStack = fV0Reader->GetMCStack();\r
+\r
+  for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
+    TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
+\r
+    if (!particle) {\r
+      //print warning here\r
+      continue;\r
+    }\r
+    \r
+    if(particle->Pt()<fV0Reader->GetPtCut()){\r
+      continue;\r
+    }\r
+    \r
+    if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut()){\r
+      continue;\r
+    }\r
+\r
+    if(particle->R()>fV0Reader->GetMaxRCut()){ // cuts on distance from collision point\r
+      continue;\r
+    }\r
+\r
+    Double_t tmpPhi=particle->Phi();\r
+    if(particle->Phi()> TMath::Pi()){\r
+      tmpPhi = particle->Phi()-(2*TMath::Pi());\r
+    }\r
+\r
+    \r
+    //process the gammas\r
+    if (particle->GetPdgCode()== 22){\r
+      fMCAllGammas.push_back(particle);\r
+      if(particle->GetMother(0)>-1){ //Means we have a mother\r
+       if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 ){//Checks for a non gamma mother.\r
+         fHistograms->FillHistogram("MC_Gamma_Energy", particle->Energy());\r
+         fHistograms->FillHistogram("MC_Gamma_Pt", particle->Pt());\r
+         fHistograms->FillHistogram("MC_Gamma_Eta", particle->Eta());\r
+         \r
+         fHistograms->FillHistogram("MC_Gamma_Phi", tmpPhi);\r
+\r
+         //adding the conversion points from all gammas with e+e- daughters\r
+         if(particle->GetNDaughters() >= 2){\r
+           TParticle* daughter0 = NULL;\r
+           TParticle* daughter1 = NULL;\r
+           \r
+           for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
+             TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
+             if(tmpDaughter->GetUniqueID() == 5){\r
+               if(tmpDaughter->GetPdgCode() == 11){\r
+                  daughter0 = tmpDaughter;\r
+               }\r
+               else if(tmpDaughter->GetPdgCode() == -11){\r
+                 daughter1 = tmpDaughter;\r
+               }\r
+             }\r
+           }\r
+\r
+           if(daughter0 == NULL || daughter1 == NULL){ // means we do not have two daughters from pair production\r
+             continue;\r
+           }\r
+\r
+           if(daughter0->R()>fV0Reader->GetMaxRCut() || daughter1->R()>fV0Reader->GetMaxRCut()){\r
+             continue;\r
+           }\r
+           \r
+           if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode()) == 11 ||\r
+              (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
+\r
+             // begin Mapping \r
+             Int_t rBin    = fHistograms->GetRBin(daughter0->R());\r
+             Int_t phiBin  = fHistograms->GetPhiBin(daughter0->Phi());\r
+             \r
+             TString nameMCMappingPhiR="";\r
+             nameMCMappingPhiR.Form("MC_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
+             fHistograms->FillHistogram(nameMCMappingPhiR, daughter0->Vz(), particle->Eta());\r
+             \r
+             TString nameMCMappingPhi="";\r
+             nameMCMappingPhi.Form("MC_EP_Mapping-Phi%02d",phiBin);\r
+             fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
+             \r
+             TString nameMCMappingR="";\r
+             nameMCMappingR.Form("MC_EP_Mapping-R%02d",rBin);\r
+             fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
+             \r
+             TString nameMCMappingPhiInR="";\r
+             nameMCMappingPhiInR.Form("MC_EP_Mapping_Phi_R-%02d",rBin);\r
+             fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
+             //end mapping\r
+\r
+             fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
+             fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
+             fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
+             fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
+           }// end if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode()) == 11 ||....... approx 20 lines above\r
+         }// end if(particle->GetNDaughters() >= 2){\r
+       } // end if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 )\r
+       if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
+           fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
+           fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2\r
+           ){ \r
+         fMCGammaChic.push_back(particle);\r
+       }\r
+      }// end if(particle->GetMother(0)>-1)\r
+      else{//means we have a primary particle\r
+       fHistograms->FillHistogram("MC_DirectGamma_Energy",particle->Energy());\r
+       fHistograms->FillHistogram("MC_DirectGamma_Pt", particle->Pt());\r
+       fHistograms->FillHistogram("MC_DirectGamma_Eta", particle->Eta());\r
+       fHistograms->FillHistogram("MCDirectGammaPhi", tmpPhi);\r
+\r
+       //adding the conversion points from all gammas with e+e- daughters\r
+       if(particle->GetNDaughters() == 2){\r
+         TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
+         TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
+         if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode() == 11) ||\r
+            (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
+           \r
+           fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
+           fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
+           fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
+           fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
+\r
+         }\r
+       }\r
+      }// end else\r
+    }// end if (particle->GetPdgCode()== 22){\r
+    else if (TMath::Abs(particle->GetPdgCode())== 11){ // Means we have an electron or a positron\r
+      if(particle->GetMother(0)>-1){ // means we have a mother\r
+       if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==22 ){ // Means we have a gamma mother\r
+         if(particle->GetPdgCode() == 11){//electron \r
+           fHistograms->FillHistogram("MC_E_Energy", particle->Energy());\r
+           fHistograms->FillHistogram("MC_E_Pt", particle->Pt());\r
+           fHistograms->FillHistogram("MC_E_Eta", particle->Eta());\r
+           fHistograms->FillHistogram("MC_E_Phi", tmpPhi);\r
+         }\r
+         if(particle->GetPdgCode() == -11){//positron \r
+           fHistograms->FillHistogram("MC_P_Energy", particle->Energy());\r
+           fHistograms->FillHistogram("MC_P_Pt", particle->Pt());\r
+           fHistograms->FillHistogram("MC_P_Eta", particle->Eta());\r
+           fHistograms->FillHistogram("MC_P_Phi", tmpPhi);\r
+         }\r
+       }\r
+      }\r
+    } // end else if (TMath::Abs(particle->GetPdgCode())== 11)\r
+    else if(particle->GetNDaughters() == 2){\r
+\r
+      TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
+      TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
+      if(daughter0->GetPdgCode() == 22 && daughter1->GetPdgCode() == 22){//check for gamma gamma daughters\r
+       \r
+       if(particle->GetPdgCode()==111){//Pi0\r
+         \r
+         fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
+\r
+         if( iTracks >= fStack->GetNprimary()){\r
+           \r
+           fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
+\r
+           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
+           fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
+           fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
+           fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
+           fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
+           fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
+           fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
+         }\r
+         else{\r
+           fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());\r
+\r
+           fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
+           fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
+           fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
+           fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
+           fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
+           fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
+           fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
+         }\r
+       }\r
+       else if(particle->GetPdgCode()==221){//Eta\r
+         fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
+\r
+         fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
+         fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
+         fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
+         fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
+         fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
+         fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
+         fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
+       }\r
+       \r
+       //the match data should be filled no matter which mother the gamma-gamma comes from\r
+       fHistograms->FillHistogram("MC_Match_Gamma_R", particle->R());\r
+       fHistograms->FillHistogram("MC_Match_Gamma_ZR", particle->Vz(),particle->R());\r
+       fHistograms->FillHistogram("MC_Match_Gamma_XY", particle->Vx(),particle->Vy());\r
+       fHistograms->FillHistogram("MC_Match_Gamma_Mass", particle->GetCalcMass());\r
+       fHistograms->FillHistogram("MC_Match_Gamma_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
+       fHistograms->FillHistogram("MC_Match_Gamma_Energy", particle->Energy());\r
+       fHistograms->FillHistogram("MC_Match_Gamma_Pt", particle->Pt());\r
+       fHistograms->FillHistogram("MC_Match_Gamma_Eta", particle->Eta());\r
+       fHistograms->FillHistogram("MC_Match_Gamma_Phi",tmpPhi);\r
+      }\r
+    }// end else if(particle->GetNDaughters() == 2)\r
+  }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
+} // end ProcessMCData\r
+\r
+void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
+  // see header file for documentation\r
+\r
+  Int_t nSurvivingV0s=0;\r
+  while(fV0Reader->NextV0()){\r
+    nSurvivingV0s++;\r
+    //-------------------------- filling v0 information -------------------------------------\r
+    fHistograms->FillHistogram("ESD_EP_OpeningAngle", fV0Reader->GetOpeningAngle());    \r
+    fHistograms->FillHistogram("ESD_EP_R", fV0Reader->GetXYRadius());\r
+    fHistograms->FillHistogram("ESD_EP_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
+    fHistograms->FillHistogram("ESD_EP_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
+    \r
+    \r
+    fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
+    fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
+    fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
+    fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
+    \r
+    fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
+    fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
+    fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
+    fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
+    \r
+    fHistograms->FillHistogram("ESD_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
+    fHistograms->FillHistogram("ESD_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
+    fHistograms->FillHistogram("ESD_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
+    fHistograms->FillHistogram("ESD_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
+\r
+\r
+    // begin mapping\r
+    Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
+    Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
+    Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
+    \r
+    TString nameESDMappingPhiR="";\r
+    nameESDMappingPhiR.Form("ESD_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
+    fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
+\r
+    TString nameESDMappingPhi="";\r
+    nameESDMappingPhi.Form("ESD_EP_Mapping-Phi%02d",phiBin);\r
+    fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
+\r
+    TString nameESDMappingR="";\r
+    nameESDMappingR.Form("ESD_EP_Mapping-R%02d",rBin);\r
+    fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  \r
+\r
+    TString nameESDMappingPhiInR="";\r
+    nameESDMappingPhiInR.Form("ESD_EP_Mapping_Phi_R-%02d",rBin);\r
+    fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
+    // end mapping\r
+    \r
+    fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
+\r
+    //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
+    if(fDoMCTruth){\r
+      if(fV0Reader->HasSameMCMother() == kFALSE){\r
+       continue;\r
+      }\r
\r
+      TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
+      TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
+\r
+      if(negativeMC->GetPdgCode()!=11 || positiveMC->GetPdgCode()!=-11){\r
+       continue;\r
+      }\r
+      if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_OpeningAngle", fV0Reader->GetOpeningAngle());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
+\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_Mass", fV0Reader->GetMotherCandidateMass());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_Width", fV0Reader->GetMotherCandidateWidth());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_R", fV0Reader->GetXYRadius());\r
+       fHistograms->FillHistogram("ESD_Match_Gamma_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
+\r
+       //resolution\r
+       Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();\r
+       Double_t esdpt  = fV0Reader->GetMotherCandidatePt();\r
+       Double_t resdPt = 0;\r
+       if(mcpt != 0){\r
+         resdPt = ((esdpt - mcpt)/mcpt)*100;\r
+       }\r
+\r
+       fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
+       fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
+       fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
+       \r
+       Double_t resdZ = 0;\r
+       if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
+         resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
+       }\r
+       \r
+       fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
+       fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
+       fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
+       \r
+       Double_t resdR = 0;\r
+       if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
+         resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
+       }\r
+       fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
+       fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
+       fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
+       fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
+      }\r
+    }\r
+  }\r
+  fHistograms->FillHistogram("NumberOfSurvivingV0s", nSurvivingV0s);\r
+  fHistograms->FillHistogram("NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
+}\r
+\r
+void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
+  // see header file for documentation\r
+\r
+  for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
+    for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
+      AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
+      AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
+\r
+      \r
+      AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
+\r
+      Double_t massTwoGammaCandidate =0.;\r
+      Double_t widthTwoGammaCandidate = 0.;\r
+      Double_t chi2TwoGammaCandidate =10000.;  \r
+      twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
+      if(twoGammaCandidate->GetNDF()>0){\r
+       chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
+       if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){\r
+\r
+         TVector3 vectorTwoGammaCandidate(twoGammaCandidate->Px(),twoGammaCandidate->Py(),twoGammaCandidate->Pz());\r
+\r
+         Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);\r
+\r
+         //Calculating by hand the radius\r
+         Double_t tmpX= twoGammaCandidate->GetX();\r
+         Double_t tmpY= twoGammaCandidate->GetY();\r
+         \r
+         Double_t radiusTwoGammaCandidate = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
+\r
+         fHistograms->FillHistogram("ESD_TwoGammaCombination_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
+         fHistograms->FillHistogram("ESD_TwoGammaCombination_Energy", twoGammaCandidate->GetE());\r
+         fHistograms->FillHistogram("ESD_TwoGammaCombination_Pt", sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
+         fHistograms->FillHistogram("ESD_TwoGammaCombination_Eta", vectorTwoGammaCandidate.Eta());\r
+         fHistograms->FillHistogram("ESD_TwoGammaCombination_Phi", vectorTwoGammaCandidate.Phi());\r
+         fHistograms->FillHistogram("ESD_TwoGammaCombination_Mass", massTwoGammaCandidate);\r
+         fHistograms->FillHistogram("ESD_TwoGammaCombination_R", radiusTwoGammaCandidate);\r
+         fHistograms->FillHistogram("ESD_TwoGammaCombination_ZR", tmpY, radiusTwoGammaCandidate);\r
+         fHistograms->FillHistogram("ESD_TwoGammaCombination_XY", tmpX, tmpY);\r
+         fHistograms->FillHistogram("InvMass_vs_Pt_Spectra",massTwoGammaCandidate ,sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
+       }\r
+      }\r
+      delete twoGammaCandidate;\r
+    }\r
+  }\r
+}\r
+\r
+void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
+  // see header file for documentation\r
+\r
+  vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
+  vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
+  for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
+    AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
+    for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
+      AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
+\r
+      AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
+\r
+      Double_t massBG =0.;\r
+      Double_t widthBG = 0.;\r
+      Double_t chi2BG =10000.; \r
+      backgroundCandidate->GetMass(massBG,widthBG);\r
+      if(backgroundCandidate->GetNDF()>0){\r
+       chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
+       if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
+\r
+         TVector3 vectorBGCandidate(backgroundCandidate->Px(),backgroundCandidate->Py(),backgroundCandidate->Pz());\r
+\r
+         Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
+\r
+         //Calculating by hand the radius (find a better way)\r
+         Double_t tmpX= backgroundCandidate->GetX();\r
+         Double_t tmpY= backgroundCandidate->GetY();\r
+         \r
+         Double_t radiusBG = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
+\r
+         fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
+         fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
+         fHistograms->FillHistogram("ESD_Background_Pt", sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));\r
+         fHistograms->FillHistogram("ESD_Background_Eta", vectorBGCandidate.Eta());\r
+         fHistograms->FillHistogram("ESD_Background_Phi", vectorBGCandidate.Phi());\r
+         fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
+         fHistograms->FillHistogram("ESD_Background_R", radiusBG);\r
+         fHistograms->FillHistogram("ESD_Background_ZR", tmpY, radiusBG);\r
+         fHistograms->FillHistogram("ESD_Background_XY", tmpX, tmpY);\r
+         fHistograms->FillHistogram("Background_InvMass_vs_Pt_Spectra",massBG,sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));\r
+       }\r
+      }\r
+      delete backgroundCandidate;   \r
+    }\r
+  }\r
+}\r
+\r
+void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
+{\r
+  // Terminate analysis\r
+  //\r
+  AliDebug(1,"Do nothing in Terminate");\r
+}\r
+\r
+void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
+{\r
+  // Create the output container\r
+  if(fOutputContainer != NULL){\r
+    delete fOutputContainer;\r
+    fOutputContainer = NULL;\r
+  }\r
+  if(fOutputContainer == NULL){\r
+    fOutputContainer = new TList();\r
+  }\r
+  fHistograms->GetOutputContainer(fOutputContainer);\r
+  fOutputContainer->SetName(GetName());\r
+}\r
+\r
+Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* daughter0, TParticle* daughter1) const{\r
+  //helper function\r
+  TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
+  TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
+  return v3D0.Angle(v3D1);\r
+}\r
diff --git a/PWG4/GammaConv/AliAnalysisTaskGammaConversion.h b/PWG4/GammaConv/AliAnalysisTaskGammaConversion.h
new file mode 100644 (file)
index 0000000..9e48cdf
--- /dev/null
@@ -0,0 +1,95 @@
+#ifndef ALIANALYSISTASKGAMMACONVERSION_H\r
+#define ALIANALYSISTASKGAMMACONVERSION_H\r
\r
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice                               */\r
+\r
+////////////////////////////////////////////////\r
+//--------------------------------------------- \r
+// Class used to do analysis on conversion pairs\r
+//---------------------------------------------\r
+////////////////////////////////////////////////\r
\r
+#include "AliAnalysisTaskSE.h"\r
+#include <vector>\r
+#include "AliV0Reader.h"\r
+\r
+class AliGammaConversionHistograms;\r
+class AliESDv0;\r
+class AliKFParticle;\r
+class AliESDInputHandler;\r
+class AliESDEvent;\r
+class AliAODEvent;\r
+class TList;\r
+class AliStack;\r
+\r
+class AliAnalysisTaskGammaConversion : public AliAnalysisTaskSE\r
+{\r
+ public:\r
+    AliAnalysisTaskGammaConversion();\r
+    AliAnalysisTaskGammaConversion(const char* name);\r
+    virtual ~AliAnalysisTaskGammaConversion() ;// virtual destructor\r
\r
+    // Implementation of interface methods\r
+    virtual void UserCreateOutputObjects();\r
+    virtual void Init();\r
+    virtual void LocalInit() {Init();}\r
+    virtual void Exec(Option_t *option);\r
+    virtual void Terminate(Option_t *option);\r
+    virtual void ConnectInputData(Option_t *);\r
+       \r
+    void ProcessMCData();\r
+    void ProcessV0s();\r
+    void ProcessGammasForNeutralMesonAnalysis();\r
+    void SetHistograms(AliGammaConversionHistograms *histograms){fHistograms=histograms;}\r
+    void SetDoMCTruth(Bool_t flag){fDoMCTruth=flag;}\r
+    void SetElectronMass(Double_t electronMass){fElectronMass = electronMass;}\r
+    void SetGammaMass(Double_t gammaMass){fGammaMass = gammaMass;}\r
+    void SetGammaWidth(Double_t gammaWidth){fGammaWidth = gammaWidth;}\r
+    void SetPi0Mass(Double_t pi0Mass){fPi0Mass = pi0Mass;}\r
+    void SetPi0Width(Double_t pi0Width){fPi0Width = pi0Width;}\r
+    void SetEtaMass(Double_t etaMass){fEtaMass = etaMass;}\r
+    void SetEtaWidth(Double_t etaWidth){fEtaWidth = etaWidth;}\r
+    void SetV0Reader(AliV0Reader* reader){fV0Reader=reader;}\r
+    void SetCalculateBackground(Bool_t bg){fCalculateBackground=bg;}\r
+    void CalculateBackground();\r
+    Double_t GetMCOpeningAngle(TParticle* daughter0, TParticle* daughter1) const;\r
+\r
+ private:\r
+    AliAnalysisTaskGammaConversion(const AliAnalysisTaskGammaConversion&); // Not implemented\r
+    AliAnalysisTaskGammaConversion& operator=(const AliAnalysisTaskGammaConversion&); // Not implemented\r
+\r
+    AliV0Reader* fV0Reader;\r
+\r
+    AliStack * fStack;\r
+\r
+    TList * fOutputContainer ; // Histogram container\r
+\r
+    AliGammaConversionHistograms *fHistograms;\r
+\r
+    Bool_t fDoMCTruth;\r
+    \r
+    vector<TParticle*> fMCAllGammas;\r
+    vector<TParticle*> fMCPi0s;\r
+    vector<TParticle*> fMCEtas;\r
+    vector<TParticle*> fMCGammaChic;\r
+\r
+    vector<AliKFParticle> fKFReconstructedGammas;\r
+\r
+    //mass defines\r
+    Double_t fElectronMass;\r
+    Double_t fGammaMass;\r
+    Double_t fPi0Mass;\r
+    Double_t fEtaMass;\r
+\r
+    // width defines\r
+    Double_t fGammaWidth;\r
+    Double_t fPi0Width;\r
+    Double_t fEtaWidth;\r
+    Bool_t fCalculateBackground;\r
+\r
+\r
+    ClassDef(AliAnalysisTaskGammaConversion, 0); // Analysis task for gamma conversions\r
+};\r
\r
+#endif //ALIANALYSISTASKGAMMA_H\r
diff --git a/PWG4/GammaConv/AliGammaConversionHistograms.cxx b/PWG4/GammaConv/AliGammaConversionHistograms.cxx
new file mode 100644 (file)
index 0000000..08ada54
--- /dev/null
@@ -0,0 +1,304 @@
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt                        *\r
+ * Version 1.0                                                            *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+////////////////////////////////////////////////\r
+//--------------------------------------------- \r
+// Class used to do analysis on conversion pairs\r
+//---------------------------------------------\r
+////////////////////////////////////////////////\r
+\r
+#include "AliGammaConversionHistograms.h"\r
+#include "TMath.h"\r
+#include "TObjString.h"\r
+#include "TMap.h"\r
+#include "TList.h"\r
+#include "TH1F.h"\r
+#include "TH2F.h"\r
+\r
+\r
+using namespace std;\r
+\r
+ClassImp(AliGammaConversionHistograms)\r
+\r
+\r
+AliGammaConversionHistograms::AliGammaConversionHistograms() :\r
+  fHistogramMap(new TMap()),\r
+  fNPhiIndex(0),\r
+  fNRIndex(0),\r
+  fMinRadius(0.),\r
+  fMaxRadius(0.),\r
+  fDeltaR(0.),\r
+  fMinPhi(0.),\r
+  fMaxPhi(0.),\r
+  fDeltaPhi(0.),\r
+  fMappingContainer(NULL)\r
+{\r
+  // see header file for documenation\r
+}\r
+\r
+\r
+AliGammaConversionHistograms::AliGammaConversionHistograms(const AliGammaConversionHistograms & original) :\r
+  fHistogramMap(original.fHistogramMap),\r
+  fNPhiIndex(original.fNPhiIndex),\r
+  fNRIndex(original.fNRIndex),\r
+  fMinRadius(original.fMinRadius),\r
+  fMaxRadius(original.fMaxRadius),\r
+  fDeltaR(original.fDeltaR),\r
+  fMinPhi(original.fMinPhi),\r
+  fMaxPhi(original.fMaxPhi),\r
+  fDeltaPhi(original.fDeltaPhi),\r
+  fMappingContainer(original.fMappingContainer)\r
+{    \r
+  //see header file for documentation\r
+}\r
+\r
+\r
+AliGammaConversionHistograms & AliGammaConversionHistograms::operator = (const AliGammaConversionHistograms & /*original*/)\r
+{\r
+  // assignment operator\r
+  return *this;\r
+}\r
+\r
+\r
+AliGammaConversionHistograms::~AliGammaConversionHistograms() {\r
+  //destructor\r
+\r
+\r
+}\r
+\r
+void AliGammaConversionHistograms::AddHistogram(TString histogramName, TString histogramTitle, Int_t nXBins, Double_t firstX,Double_t lastX,TString xAxisTitle, TString yAxisTitle){\r
+  // see header file for documentation\r
+  TH1F *tmp = new TH1F(histogramName, histogramTitle,nXBins,firstX,lastX);\r
+  tmp->GetXaxis()->SetTitle(xAxisTitle);\r
+  tmp->GetYaxis()->SetTitle(yAxisTitle);\r
+  TObjString* tobjstring = new TObjString(histogramName.Data());\r
+  fHistogramMap->Add((TObject*)tobjstring,(TObject*)tmp);\r
+}\r
+\r
+void AliGammaConversionHistograms::AddHistogram(TString histogramName, TString histogramTitle, Int_t nXBins, Double_t firstX, Double_t lastX, Int_t nYBins, Double_t firstY, Double_t lastY, TString xAxisTitle, TString yAxisTitle){\r
+  // see header file for documentation\r
+  TH2F *tmp = new TH2F(histogramName, histogramTitle,nXBins,firstX,lastX,nYBins,firstY,lastY);\r
+  tmp->GetXaxis()->SetTitle(xAxisTitle);\r
+  tmp->GetYaxis()->SetTitle(yAxisTitle);\r
+  TObjString *tobjstring = new TObjString(histogramName.Data());\r
+  fHistogramMap->Add((TObject*)tobjstring,(TObject*)tmp);\r
+}\r
+\r
+void AliGammaConversionHistograms::FillHistogram(TString histogramName, Double_t xValue) const{\r
+  //see header file for documentation\r
+  TH1 *tmp = (TH1*)fHistogramMap->GetValue(histogramName.Data());\r
+  if(tmp){\r
+      tmp->Fill(xValue);\r
+  }\r
+  else{\r
+    cout<<"Histogram does not exist: "<<histogramName.Data()<<endl;\r
+  }\r
+}\r
+\r
+void AliGammaConversionHistograms::FillHistogram(TString histogramName, Double_t xValue, Double_t yValue) const{\r
+  //see header file for documentation\r
+  TH1 *tmp = (TH1*)fHistogramMap->GetValue(histogramName.Data());\r
+  if(tmp){\r
+    tmp->Fill(xValue, yValue);\r
+  }\r
+  else{\r
+    cout<<"Histogram does not exist: "<<histogramName.Data()<<endl;\r
+  }\r
+}\r
+\r
+void AliGammaConversionHistograms::GetOutputContainer(TList *fOutputContainer){\r
+  //checking if the container is alrerady created\r
+\r
+  if(fOutputContainer == NULL){\r
+    //print warning\r
+    return;\r
+  }\r
+\r
+  if(fHistogramMap != NULL){\r
+    TIter iter(fHistogramMap);\r
+    TObjString *histogramName;\r
+    while ((histogramName = (TObjString*) iter.Next())) {\r
+      TString histogramString = histogramName->GetString();\r
+      if(histogramString.Contains("Mapping")){// means it should be put in the mapping folder\r
+       if(fMappingContainer == NULL){\r
+         fMappingContainer = new TList();\r
+         fMappingContainer->SetName("Mapping histograms");\r
+       }\r
+       if(fMappingContainer != NULL){\r
+         fMappingContainer->Add((TH1*)fHistogramMap->GetValue(histogramString.Data()));\r
+       }\r
+      }\r
+      fOutputContainer->Add((TH1*)fHistogramMap->GetValue(histogramString.Data()));\r
+      histogramName = NULL;\r
+    } // end while\r
+    //   fOutputContainer->Add(fMappingContainer);\r
+  }\r
+}\r
+\r
+Int_t AliGammaConversionHistograms::GetRBin(Double_t radius) const{\r
+  // see header file for documentation\r
+  Int_t iResult=0;\r
+  if(fDeltaR>0){\r
+    iResult = (Int_t)((radius - fMinRadius)/fDeltaR);\r
+  }\r
+  return iResult;\r
+}\r
+\r
+Int_t AliGammaConversionHistograms::GetPhiBin(Double_t phi) const{\r
+  // see header file for documentation\r
+  Int_t iResult=0;\r
+  if(fDeltaPhi>0){\r
+    if(phi>TMath::Pi()){\r
+      phi-=2*TMath::Pi();\r
+    }\r
+    iResult = (Int_t)((phi - fMinPhi)/fDeltaPhi);\r
+  }\r
+  return iResult;\r
+}\r
+\r
+\r
+\r
+void AliGammaConversionHistograms::InitializeMappingValues(Int_t nPhiIndex, Int_t nRIndex, Int_t nBinsR, Double_t minRadius, Double_t maxRadius,Int_t nBinsPhi, Double_t minPhi, Double_t maxPhi){\r
+  // Initializing the valuse for the mapping\r
+\r
+  fNPhiIndex = nPhiIndex;\r
+  fNRIndex   = nRIndex;\r
+  fMinRadius      = minRadius;\r
+  fMaxRadius      = maxRadius;\r
+  if(nBinsR>0 && nRIndex!=0){\r
+    fDeltaR       = (fMaxRadius - fMinRadius)/nRIndex;\r
+  }\r
+  fMinPhi         = minPhi;\r
+  fMaxPhi         = maxPhi;\r
+  if(nBinsPhi>0 && nPhiIndex!=0){\r
+    fDeltaPhi     = (fMaxPhi-fMinPhi)/nPhiIndex;\r
+  }\r
+}\r
+\r
+\r
+//mapping\r
+void AliGammaConversionHistograms::AddMappingHistograms(Int_t nPhiIndex, Int_t nRIndex,Int_t nXBins, Double_t firstX, Double_t lastX, Int_t nYBins, Double_t firstY, Double_t lastY, TString xAxisTitle, TString yAxisTitle){\r
+  // see header file for documentation\r
+  \r
+  for(Int_t phi =0; phi<=fNPhiIndex;phi++){\r
+\r
+    for(Int_t r =0; r<fNRIndex;r++){\r
+\r
+      // setting axis to "" changes below\r
+      xAxisTitle="";\r
+      yAxisTitle="";\r
+      //Creating the axis titles\r
+      if(xAxisTitle.Length() == 0){\r
+       xAxisTitle.Form("Phi %02d",phi);\r
+      }\r
+      \r
+      if(yAxisTitle.Length() == 0){\r
+       yAxisTitle.Form("R %02d",phi);\r
+      }\r
+\r
+      //MC\r
+      TString nameMC="";\r
+      nameMC.Form("MC_EP_Mapping-Phi%02d-R%02d",phi,r);\r
+      TString titleMC="";\r
+      titleMC.Form("Electron-Positron MC Mapping-Phi%02d-R%02d",phi,r);\r
+\r
+      AddHistogram(nameMC, titleMC, nXBins, firstX, lastX, nYBins, firstY, lastY, xAxisTitle, yAxisTitle);\r
+\r
+      //ESD\r
+      TString nameESD="";\r
+      nameESD.Form("ESD_EP_Mapping-Phi%02d-R%02d",phi,r);\r
+      TString titleESD="";\r
+      titleESD.Form("Electron-Positron ESD Mapping-Phi%02d-R%02d",phi,r);\r
+\r
+      AddHistogram(nameESD, titleESD, nXBins, firstX, lastX, nYBins, firstY, lastY, xAxisTitle, yAxisTitle);\r
+    }\r
+  }\r
+\r
+\r
+  for(Int_t phi =0; phi<=nPhiIndex;phi++){ \r
+\r
+    // setting axis to "" changes below\r
+    xAxisTitle="";\r
+    yAxisTitle="";\r
+    //Creating the axis titles\r
+    if(xAxisTitle.Length() == 0){\r
+      xAxisTitle.Form("Phi %02d",phi);\r
+    }\r
+    if(yAxisTitle.Length() == 0){\r
+      yAxisTitle = "Counts";\r
+    }\r
+    \r
+    //MC\r
+    TString nameMC="";\r
+    nameMC.Form("MC_EP_Mapping-Phi%02d",phi);\r
+    TString titleMC="";\r
+    titleMC.Form("Electron-Positron MC Mapping-Phi%02d",phi);\r
+    \r
+    AddHistogram(nameMC, titleMC, nXBins, firstX, lastX, xAxisTitle, yAxisTitle);\r
+\r
+    //MC\r
+    TString nameESD="";\r
+    nameESD.Form("ESD_EP_Mapping-Phi%02d",phi);\r
+    TString titleESD="";\r
+    titleESD.Form("Electron-Positron ESD Mapping-Phi%02d",phi);\r
+    \r
+    AddHistogram(nameESD, titleESD, nXBins, firstX, lastX, xAxisTitle, yAxisTitle);\r
+  }\r
+\r
+\r
+  for(Int_t r =0; r<=nRIndex;r++){\r
+\r
+    // setting axis to "" changes below\r
+    xAxisTitle="";\r
+    yAxisTitle="";\r
+    //Creating the axis titles\r
+    if(xAxisTitle.Length() == 0){\r
+      xAxisTitle.Form("R %02d",r);\r
+    }\r
+    if(yAxisTitle.Length() == 0){\r
+      yAxisTitle = "Counts";\r
+    }\r
+    \r
+    //MC\r
+    TString nameMC="";\r
+    nameMC.Form("MC_EP_Mapping-R%02d",r);\r
+    TString titleMC="";\r
+    titleMC.Form("Electron-Positron MC Mapping-R%02d",r);\r
+    \r
+    AddHistogram(nameMC, titleMC, nXBins, firstX, lastX, xAxisTitle, yAxisTitle);\r
+\r
+    //ESD\r
+    TString nameESD="";\r
+    nameESD.Form("ESD_EP_Mapping-R%02d",r);\r
+    TString titleESD="";\r
+    titleESD.Form("Electron-Positron ESD Mapping-R%02d",r);\r
+    \r
+    AddHistogram(nameESD, titleESD, nXBins, firstX, lastX, xAxisTitle, yAxisTitle);\r
+\r
+    //Mapping Phi in R\r
+    TString nameMCPhiInR="";\r
+    nameMCPhiInR.Form("MC_EP_Mapping_Phi_R-%02d",r);\r
+    TString titleMCPhiInR="";\r
+    titleMCPhiInR.Form("MC Mapping of Phi in R%02d",r);\r
+    AddHistogram(nameMCPhiInR, titleMCPhiInR, nXBins, firstX, lastX, xAxisTitle, yAxisTitle);\r
+\r
+    //Mapping Phi in R\r
+    TString nameESDPhiInR="";\r
+    nameESDPhiInR.Form("ESD_EP_Mapping_Phi_R-%02d",r);\r
+    TString titleESDPhiInR="";\r
+    titleESDPhiInR.Form("ESD Mapping of Phi in R%02d",r);\r
+    AddHistogram(nameESDPhiInR, titleESDPhiInR, nXBins, firstX, lastX, xAxisTitle, yAxisTitle);    \r
+  }\r
+}\r
diff --git a/PWG4/GammaConv/AliGammaConversionHistograms.h b/PWG4/GammaConv/AliGammaConversionHistograms.h
new file mode 100644 (file)
index 0000000..8434d66
--- /dev/null
@@ -0,0 +1,83 @@
+#ifndef ALIGAMMACONVERSIONHISTOGRAMS_H\r
+#define ALIGAMMACONVERSIONHISTOGRAMS_H\r
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice     */\r
+\r
+////////////////////////////////////////////////\r
+//--------------------------------------------- \r
+// Class used to do analysis on conversion pairs\r
+//---------------------------------------------\r
+////////////////////////////////////////////////\r
+\r
+#include "TString.h"\r
+#include "Riostream.h"\r
+#include <vector>\r
+\r
+class TMap;\r
+class TList;\r
+class TH1F;\r
+class TH2F;\r
+\r
+class AliGammaConversionHistograms{\r
+\r
+ public: \r
+  \r
+  AliGammaConversionHistograms();                                                                  //constructor\r
+  AliGammaConversionHistograms(const AliGammaConversionHistograms & original);                     //copy constructor\r
+  AliGammaConversionHistograms & operator = (const AliGammaConversionHistograms & original);       //assignment operator\r
+  virtual ~AliGammaConversionHistograms();                                                         //virtual destructor\r
+  \r
+\r
+  //  TList * GetOutputContainer();\r
+  void GetOutputContainer(TList *fOutputContainer);\r
+  \r
+  Int_t GetRBin(Double_t radius) const;\r
+  Int_t GetPhiBin(Double_t phi) const;\r
+\r
+  void InitializeMappingValues(Int_t nPhiHistograms, Int_t nRHistograms, Int_t nBinsR, Double_t minRadius, Double_t maxRadius,Int_t nBinsPhi, Double_t minPhi, Double_t maxPhi);\r
+\r
+  void AddMappingHistograms(Int_t nPhiHistograms, Int_t nRHistograms,Int_t nXBins, Double_t firstX, Double_t lastX, Int_t nYBins, Double_t firstY, Double_t lastY, TString xAxisTitle="", TString yAxisTitle="");\r
+\r
+  /*\r
+   * Adds a TH1F histogram to the histogram map and create a key for it \r
+   */\r
+  void AddHistogram(TString histogramName, TString histogramTitle, Int_t nXBins, Double_t firstX,Double_t lastX,TString xAxisTitle="", TString yAxisTitle="");\r
+\r
+  /*\r
+   * Adds a TH2F histogram to the histogram map and create a key for it \r
+   */\r
+  void AddHistogram(TString histogramName, TString histogramTitle, Int_t nXBins, Double_t firstX, Double_t lastX, Int_t nYBins, Double_t firstY, Double_t lastY, TString xAxisTitle="", TString yAxisTitle="");\r
+\r
+  /*\r
+   * Fills a TH1F histogram with the given name with the given value \r
+   */\r
+  void FillHistogram(TString histogramName, Double_t xValue) const;\r
+\r
+  /*\r
+   * Fills a TH2F histogram with the given name with the given value \r
+   */\r
+  void FillHistogram(TString histogramName, Double_t xValue, Double_t yValue) const;\r
+\r
+ private:\r
+  TMap* fHistogramMap;\r
+\r
+  Int_t fNPhiIndex;\r
+  Int_t fNRIndex;\r
+  Double_t fMinRadius;\r
+  Double_t fMaxRadius;\r
+  Double_t fDeltaR;\r
+  Double_t fMinPhi;\r
+  Double_t fMaxPhi;\r
+  Double_t fDeltaPhi;\r
+\r
+  TList * fMappingContainer;\r
+\r
+  \r
+  ClassDef(AliGammaConversionHistograms,1)\r
+} ;\r
+\r
+\r
+#endif\r
+\r
+\r
+\r
diff --git a/PWG4/GammaConv/AliV0Reader.cxx b/PWG4/GammaConv/AliV0Reader.cxx
new file mode 100644 (file)
index 0000000..7e39f1f
--- /dev/null
@@ -0,0 +1,498 @@
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt                        *\r
+ * Version 1.0                                                            *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+////////////////////////////////////////////////\r
+//--------------------------------------------- \r
+// Class used to do analysis on conversion pairs\r
+//---------------------------------------------\r
+////////////////////////////////////////////////\r
+\r
+// --- ROOT system ---\r
+#include <TMath.h>\r
+\r
+//---- ANALYSIS system ----\r
+#include "AliV0Reader.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliESDInputHandler.h"\r
+#include "AliMCEvent.h"\r
+#include "AliKFVertex.h"\r
+\r
+#include "AliStack.h"\r
+#include "AliMCEventHandler.h"\r
+\r
+\r
+class iostream;\r
+class AliESDv0;\r
+class TFormula;\r
+\r
+using namespace std;\r
+\r
+ClassImp(AliV0Reader)\r
+\r
+\r
+\r
+  AliV0Reader::AliV0Reader() :\r
+    TObject(),\r
+    fMCStack(NULL),\r
+    fMCTruth(NULL),\r
+    fChain(NULL),\r
+    fESDHandler(NULL),\r
+    fESDEvent(NULL),\r
+    fHistograms(NULL),\r
+    fCurrentV0IndexNumber(0),\r
+    fCurrentV0(NULL),\r
+    fCurrentNegativeKFParticle(NULL),\r
+    fCurrentPositiveKFParticle(NULL),\r
+    fCurrentMotherKFCandidate(NULL),\r
+    fCurrentNegativeESDTrack(NULL),\r
+    fCurrentPositiveESDTrack(NULL),\r
+    fNegativeTrackLorentzVector(NULL),\r
+    fPositiveTrackLorentzVector(NULL),\r
+    fMotherCandidateLorentzVector(NULL),\r
+    fCurrentXValue(0),\r
+    fCurrentYValue(0),\r
+    fCurrentZValue(0),\r
+    fPositiveTrackPID(0),\r
+    fNegativeTrackPID(0),\r
+    fNegativeMCParticle(NULL),\r
+    fPositiveMCParticle(NULL),\r
+    fMotherMCParticle(NULL),\r
+    fMotherCandidateKFMass(0),\r
+    fMotherCandidateKFWidth(0),\r
+    fUseKFParticle(kTRUE),\r
+    fUseESDTrack(kFALSE),\r
+    fDoMC(kFALSE),\r
+    fMaxR(10000),// 100 meter(outside of ALICE)\r
+    fEtaCut(0.),\r
+    fPtCut(0.),\r
+    fChi2CutConversion(0.),\r
+    fChi2CutMeson(0.),\r
+    fPIDProbabilityCutNegativeParticle(0),\r
+    fPIDProbabilityCutPositiveParticle(0),\r
+    fXVertexCut(0.),\r
+    fYVertexCut(0.),\r
+    fZVertexCut(0.),\r
+    fNSigmaMass(0.),\r
+    fUseImprovedVertex(kFALSE),\r
+    fCurrentEventGoodV0s(),\r
+    fPreviousEventGoodV0s()\r
+{\r
+\r
+}\r
+\r
+\r
+AliV0Reader::AliV0Reader(const AliV0Reader & original) :\r
+  TObject(original),\r
+  fMCStack(original.fMCStack),\r
+  fMCTruth(original.fMCTruth),\r
+  fChain(original.fChain),\r
+  fESDHandler(original.fESDHandler),\r
+  fESDEvent(original.fESDEvent),\r
+  fHistograms(original.fHistograms),\r
+  fCurrentV0IndexNumber(original.fCurrentV0IndexNumber),\r
+  fCurrentV0(original.fCurrentV0),\r
+  fCurrentNegativeKFParticle(original.fCurrentNegativeKFParticle),\r
+  fCurrentPositiveKFParticle(original.fCurrentPositiveKFParticle),\r
+  fCurrentMotherKFCandidate(original.fCurrentMotherKFCandidate),\r
+  fCurrentNegativeESDTrack(original.fCurrentNegativeESDTrack),\r
+  fCurrentPositiveESDTrack(original.fCurrentPositiveESDTrack),\r
+  fNegativeTrackLorentzVector(original.fNegativeTrackLorentzVector),\r
+  fPositiveTrackLorentzVector(original.fPositiveTrackLorentzVector),\r
+  fMotherCandidateLorentzVector(original.fMotherCandidateLorentzVector),\r
+  fCurrentXValue(original.fCurrentXValue),\r
+  fCurrentYValue(original.fCurrentYValue),\r
+  fCurrentZValue(original.fCurrentZValue),\r
+  fPositiveTrackPID(original.fPositiveTrackPID),\r
+  fNegativeTrackPID(original.fNegativeTrackPID),\r
+  fNegativeMCParticle(original.fNegativeMCParticle),\r
+  fPositiveMCParticle(original.fPositiveMCParticle),\r
+  fMotherMCParticle(original.fMotherMCParticle),\r
+  fMotherCandidateKFMass(original.fMotherCandidateKFMass),\r
+  fMotherCandidateKFWidth(original.fMotherCandidateKFWidth),\r
+  fUseKFParticle(kTRUE),\r
+  fUseESDTrack(kFALSE),\r
+  fDoMC(kFALSE),\r
+  fMaxR(original.fMaxR),\r
+  fEtaCut(original.fEtaCut),\r
+  fPtCut(original.fPtCut),\r
+  fChi2CutConversion(original.fChi2CutConversion),\r
+  fChi2CutMeson(original.fChi2CutMeson),\r
+  fPIDProbabilityCutNegativeParticle(original.fPIDProbabilityCutNegativeParticle),\r
+  fPIDProbabilityCutPositiveParticle(original.fPIDProbabilityCutPositiveParticle),\r
+  fXVertexCut(original.fXVertexCut),\r
+  fYVertexCut(original.fYVertexCut),\r
+  fZVertexCut(original.fZVertexCut),\r
+  fNSigmaMass(original.fNSigmaMass),\r
+  fUseImprovedVertex(original.fUseImprovedVertex),\r
+  fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),\r
+  fPreviousEventGoodV0s(original.fPreviousEventGoodV0s)\r
+{\r
+\r
+}\r
+\r
+\r
+AliV0Reader & AliV0Reader::operator = (const AliV0Reader & /*source*/)\r
+{\r
+  // assignment operator\r
+  return *this;\r
+}\r
+\r
+void AliV0Reader::Initialize(){\r
+  //see header file for documentation\r
+\r
+  // Get the input handler from the manager\r
+  fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
+  if(fESDHandler == NULL){\r
+    //print warning here\r
+  }\r
+  \r
+  // Get pointer to esd event from input handler\r
+  fESDEvent = fESDHandler->GetEvent();\r
+  if(fESDEvent == NULL){\r
+    //print warning here\r
+  }\r
+\r
+  //Get pointer to MCTruth\r
+  fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());\r
+  if(fMCTruth == NULL){\r
+    //print warning here\r
+  }\r
+\r
+  //Get pointer to the mc stack\r
+  fMCStack = fMCTruth->MCEvent()->Stack();\r
+  if(fMCStack == NULL){\r
+    //print warning here\r
+  }\r
+\r
+  AliKFParticle::SetField(fESDEvent->GetMagneticField());\r
+\r
+}\r
+\r
+AliESDv0* AliV0Reader::GetV0(Int_t index){\r
+  //see header file for documentation\r
+\r
+  fCurrentV0 = fESDEvent->GetV0(index);\r
+  UpdateV0Information();\r
+  return fCurrentV0;\r
+}\r
+\r
+\r
+Bool_t AliV0Reader::NextV0(){\r
+  //see header file for documentation\r
+\r
+  Bool_t iResult=kFALSE;\r
+  while(fCurrentV0IndexNumber<fESDEvent->GetNumberOfV0s()){\r
+    fCurrentV0 = fESDEvent->GetV0(fCurrentV0IndexNumber);\r
+    \r
+    //checks if on the fly mode is set\r
+    if ( !fCurrentV0->GetOnFlyStatus() ){\r
+      fCurrentV0IndexNumber++;\r
+      fHistograms->FillHistogram("V0MassDebugCut1",GetMotherCandidateMass());\r
+      continue;\r
+    }\r
+\r
+    if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {//checks if we have a vertex\r
+      fCurrentV0IndexNumber++;\r
+      fHistograms->FillHistogram("V0MassDebugCut2",GetMotherCandidateMass());\r
+      continue;\r
+    }\r
+\r
+    if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){\r
+      fCurrentV0IndexNumber++;\r
+      fHistograms->FillHistogram("V0MassDebugCut3",GetMotherCandidateMass());\r
+      continue;\r
+    }\r
+\r
+\r
+    fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);\r
\r
+    if(GetXYRadius()>fMaxR){ // cuts on distance from collision point\r
+      fCurrentV0IndexNumber++;\r
+      fHistograms->FillHistogram("V0MassDebugCut4",GetMotherCandidateMass());\r
+      continue;\r
+    }\r
+\r
+    UpdateV0Information();\r
+        \r
+    if(fUseKFParticle){\r
+      if(fCurrentMotherKFCandidate->GetNDF()<=0){\r
+       fCurrentV0IndexNumber++;\r
+       fHistograms->FillHistogram("V0MassDebugCut5",GetMotherCandidateMass());\r
+       continue;\r
+      }\r
+      Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();\r
+      if(chi2V0 > fChi2CutConversion || chi2V0 <=0){\r
+       fCurrentV0IndexNumber++;\r
+       fHistograms->FillHistogram("V0MassDebugCut6",GetMotherCandidateMass());\r
+       continue;\r
+      }\r
+      \r
+      if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut){\r
+       fCurrentV0IndexNumber++;\r
+       fHistograms->FillHistogram("V0MassDebugCut7",GetMotherCandidateMass());\r
+       continue;\r
+      }\r
+      \r
+      if(fMotherCandidateLorentzVector->Pt()<fPtCut){\r
+       fCurrentV0IndexNumber++;\r
+       fHistograms->FillHistogram("V0MassDebugCut8",GetMotherCandidateMass());\r
+       continue;\r
+      }\r
+      \r
+    }\r
+    else if(fUseESDTrack){\r
+      //TODO\r
+    }\r
+\r
+    iResult=kTRUE;//means we have a v0 who survived all the cuts applied\r
+\r
+    fCurrentV0IndexNumber++;\r
+    \r
+    break;\r
+  }\r
+  return iResult; \r
+}\r
+\r
+void AliV0Reader::UpdateV0Information(){\r
+  //see header file for documentation\r
+  \r
+  if(fCurrentNegativeKFParticle != NULL){\r
+    delete fCurrentNegativeKFParticle;\r
+  }\r
+  fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fNegativeTrackPID);\r
+  \r
+  if(fCurrentPositiveKFParticle != NULL){\r
+    delete fCurrentPositiveKFParticle;\r
+  }\r
+  fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fPositiveTrackPID);\r
+    \r
+  if(fCurrentMotherKFCandidate != NULL){\r
+    delete fCurrentMotherKFCandidate;\r
+  }\r
+  fCurrentMotherKFCandidate = new AliKFParticle(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);\r
+\r
+  fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());\r
+\r
+  fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());\r
+\r
+  if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){\r
+    fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);\r
+  }\r
+\r
+  if(fUseImprovedVertex == kTRUE){\r
+    AliKFVertex primaryVertexImproved(*GetPrimaryVertex());\r
+    primaryVertexImproved+=*fCurrentMotherKFCandidate;\r
+    fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);\r
+  }\r
+\r
+  fCurrentMotherKFCandidate->GetMass(fMotherCandidateKFMass,fMotherCandidateKFWidth);\r
+\r
+\r
+  if(fNegativeTrackLorentzVector != NULL){\r
+    delete fNegativeTrackLorentzVector;\r
+  }\r
+  if(fUseKFParticle){\r
+    fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeKFParticle->Px(),fCurrentNegativeKFParticle->Py(),fCurrentNegativeKFParticle->Pz());\r
+  }\r
+  else if(fUseESDTrack){\r
+    fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeESDTrack->Px(),fCurrentNegativeESDTrack->Py(),fCurrentNegativeESDTrack->Pz());\r
+  }\r
+\r
+  if(fPositiveTrackLorentzVector != NULL){\r
+    delete fPositiveTrackLorentzVector;\r
+  }\r
+  if(fUseKFParticle){\r
+    fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveKFParticle->Px(),fCurrentPositiveKFParticle->Py(),fCurrentPositiveKFParticle->Pz());\r
+  }\r
+  else if(fUseESDTrack){\r
+    fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveESDTrack->Px(),fCurrentPositiveESDTrack->Py(),fCurrentPositiveESDTrack->Pz());\r
+  }\r
+\r
+  if(fMotherCandidateLorentzVector != NULL){\r
+    delete fMotherCandidateLorentzVector;\r
+  }\r
+  if(fUseKFParticle){\r
+    fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);\r
+  }\r
+  else if(fUseESDTrack){\r
+    fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);\r
+  }\r
+\r
+  if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){\r
+    fMotherCandidateLorentzVector->SetXYZM(fMotherCandidateLorentzVector->Px() ,fMotherCandidateLorentzVector->Py(),fMotherCandidateLorentzVector->Pz(),0.); \r
+  }\r
+    \r
+  if(fDoMC == kTRUE){\r
+    fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));\r
+    fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));\r
+  }\r
+  fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate);\r
+}\r
+\r
+Bool_t AliV0Reader::HasSameMCMother(){\r
+  //see header file for documentation\r
+\r
+  Bool_t iResult = kFALSE;\r
+  if(fDoMC == kTRUE){\r
+    if(fNegativeMCParticle != NULL && fPositiveMCParticle != NULL){\r
+      if(fNegativeMCParticle->GetMother(0) == fPositiveMCParticle->GetMother(0))\r
+       fMotherMCParticle = fMCStack->Particle(fPositiveMCParticle->GetMother(0));\r
+      iResult = kTRUE;\r
+    }\r
+  }\r
+  return iResult;\r
+}\r
+\r
+Bool_t AliV0Reader::CheckPIDProbability(Double_t negProbCut, Double_t posProbCut){\r
+  //see header file for documentation\r
+\r
+  Bool_t iResult=kFALSE;\r
+\r
+  Double_t *posProbArray = new Double_t[10];\r
+  Double_t *negProbArray = new Double_t[10];\r
+  AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());\r
+  AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());\r
+  \r
+  negTrack->GetTPCpid(negProbArray);\r
+  posTrack->GetTPCpid(posProbArray);\r
+\r
+  if(negProbArray!=NULL && posProbArray!=NULL){\r
+    if(negProbArray[GetSpeciesIndex(-1)]>=negProbCut && posProbArray[GetSpeciesIndex(1)]>=posProbCut){\r
+      iResult=kTRUE;\r
+    }\r
+  }\r
+  delete [] posProbArray;\r
+  delete [] negProbArray;\r
+  return iResult;\r
+}\r
+\r
+void AliV0Reader::UpdateEventByEventData(){\r
+  //see header file for documentation\r
+\r
+  if(fCurrentEventGoodV0s.size() >0 ){\r
+    fPreviousEventGoodV0s.clear();\r
+    fPreviousEventGoodV0s = fCurrentEventGoodV0s;\r
+  }\r
+  fCurrentEventGoodV0s.clear();\r
+  \r
+  fCurrentV0IndexNumber=0;\r
+}\r
+\r
+Double_t AliV0Reader::GetNegativeTrackPhi() const{\r
+  //see header file for documentation\r
+\r
+  Double_t offset=0;\r
+  if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){\r
+    offset = -2*TMath::Pi();\r
+  }\r
+  return fNegativeTrackLorentzVector->Phi()+offset;\r
+}\r
+\r
+Double_t AliV0Reader::GetPositiveTrackPhi() const{\r
+  //see header file for documentation\r
+\r
+  Double_t offset=0;\r
+  if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){\r
+    offset = -2*TMath::Pi();\r
+  }\r
+  return fPositiveTrackLorentzVector->Phi()+offset;\r
+}\r
+\r
+Double_t AliV0Reader::GetMotherCandidatePhi() const{\r
+  //see header file for documentation\r
+\r
+  Double_t offset=0;\r
+  if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){\r
+    offset = -2*TMath::Pi();\r
+  }\r
+  return fMotherCandidateLorentzVector->Phi()+offset;\r
+}\r
+\r
+Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){\r
+  //see header file for documentation\r
+\r
+  Int_t iResult = 10; // Unknown particle\r
+\r
+  if(chargeOfTrack==-1){ //negative track\r
+    switch(abs(fNegativeTrackPID)){\r
+    case 11:       //electron\r
+      iResult = 0;\r
+      break;\r
+    case 13:       //muon\r
+      iResult = 1;\r
+      break;\r
+    case 211:      //pion\r
+      iResult = 2;\r
+      break;\r
+    case 321:      //kaon\r
+      iResult = 3;\r
+      break;\r
+    case 2212:     //proton\r
+      iResult = 4;\r
+      break;\r
+    case 22:       //photon\r
+      iResult = 5;\r
+      break;\r
+    case 111:      //pi0\r
+      iResult = 6;\r
+      break;\r
+    case 2112:     //neutron\r
+      iResult = 7;\r
+      break;\r
+    case 311:      //K0\r
+      iResult = 8;\r
+      break;\r
+      \r
+      //Put in here for kSPECIES::kEleCon  ????\r
+    }\r
+  }\r
+  else if(chargeOfTrack==1){ //positive track\r
+    switch(abs(fPositiveTrackPID)){\r
+    case 11:       //electron\r
+      iResult = 0;\r
+      break;\r
+    case 13:       //muon\r
+      iResult = 1;\r
+      break;\r
+    case 211:      //pion\r
+      iResult = 2;\r
+      break;\r
+    case 321:      //kaon\r
+      iResult = 3;\r
+      break;\r
+    case 2212:     //proton\r
+      iResult = 4;\r
+      break;\r
+    case 22:       //photon\r
+      iResult = 5;\r
+      break;\r
+    case 111:      //pi0\r
+      iResult = 6;\r
+      break;\r
+    case 2112:     //neutron\r
+      iResult = 7;\r
+      break;\r
+    case 311:      //K0\r
+      iResult = 8;\r
+      break;\r
+\r
+      //Put in here for kSPECIES::kEleCon  ????\r
+    }\r
+  }\r
+  else{\r
+    //Wrong parameter.. Print warning\r
+  }\r
+  return iResult;\r
+}\r
diff --git a/PWG4/GammaConv/AliV0Reader.h b/PWG4/GammaConv/AliV0Reader.h
new file mode 100644 (file)
index 0000000..41a3ab9
--- /dev/null
@@ -0,0 +1,308 @@
+#ifndef ALIV0READER_H\r
+#define ALIV0READER_H\r
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice     */\r
+\r
+////////////////////////////////////////////////\r
+//--------------------------------------------- \r
+// Class used to do analysis on conversion pairs\r
+//---------------------------------------------\r
+////////////////////////////////////////////////\r
+\r
+// --- ROOT system ---\r
+#include "TObject.h" \r
+#include "AliESDv0.h"\r
+#include "AliESDEvent.h"\r
+#include "AliKFParticle.h"\r
+#include "TParticle.h"\r
+#include "AliGammaConversionHistograms.h"\r
+#include <vector>\r
+\r
+class TClonesArray; \r
+class TFormula;\r
+class Riostream;\r
+class TChain;\r
+//--- AliRoot system ---\r
+\r
+class AliStack;\r
+class AliESDEvent; \r
+class AliMCEventHandler;\r
+class AliESDInputHandler;\r
+class AliESDVertex;\r
+class AliLog;\r
+class TChain;\r
+class TChain;\r
+\r
+\r
+\r
+class AliV0Reader : public TObject {\r
+\r
+ public: \r
+\r
+  AliV0Reader();                                        //constructor\r
+  AliV0Reader(const AliV0Reader & g);                   //copy constructor\r
+  AliV0Reader & operator = (const AliV0Reader & g);     //assignment operator\r
+  virtual ~AliV0Reader() {;}                            //virtual destructor\r
+  /*\r
+   *Initialize the reader\r
+   */\r
+  void Initialize();\r
+\r
+  /*\r
+   *Returns the number of v0s in the event, no cuts applied.\r
+   */\r
+  Int_t GetNumberOfV0s() const{return fESDEvent->GetNumberOfV0s();}\r
+\r
+  /*\r
+   * Check if there are any more good v0s left in the v0 stack\r
+   * if so, fCurrent v0 is set to this v0 and can be retrieved\r
+   * by GetCurrentV0 function.\r
+   * returns kFALSE if there is no more good v0s in the v0 stack\r
+   */\r
+  Bool_t NextV0();\r
+  \r
+  /*\r
+   * Returns the v0 at the given index, no checks are done on the v0. \r
+   */\r
+  AliESDv0* GetV0(Int_t index);\r
+\r
+  /*\r
+   * Returns the current v0\r
+   */\r
+  AliESDv0* GetCurrentV0() const{return fCurrentV0;}\r
+\r
+  /*\r
+   * Returns the negative ESD track which belongs to fCurrentV0\r
+   */\r
+  AliESDtrack* GetNegativeESDTrack(){return fESDEvent->GetTrack(fCurrentV0->GetNindex());}\r
+\r
+  /*\r
+   * Returns the positive ESD track which belongs to fCurrentV0\r
+   */\r
+  AliESDtrack* GetPositiveESDTrack(){return fESDEvent->GetTrack(fCurrentV0->GetPindex());}\r
+\r
+  /*\r
+   * Returns the negative KF particle which belongs to fCurrentV0\r
+   */\r
+  AliKFParticle* GetNegativeKFParticle() const{return fCurrentNegativeKFParticle;}\r
+\r
+  /*\r
+   * Returns the positive KF particle which belongs to fCurrentV0\r
+   */\r
+  AliKFParticle* GetPositiveKFParticle() const{return fCurrentPositiveKFParticle;}\r
+  /*\r
+   * Returns the KFParticle object of the 2 tracks.\r
+   */\r
+  AliKFParticle* GetMotherCandidateKFCombination() const{return fCurrentMotherKFCandidate;}\r
+  /*\r
+   * Checks the probablity that the PID of the particle is what we want it to be.\r
+   */\r
+  Bool_t CheckPIDProbability(Double_t negProbCut, Double_t posProbCut);\r
+\r
+  /*\r
+   *Get the negative MC TParticle from the stack \r
+   */\r
+  TParticle * GetNegativeMCParticle() const{return fNegativeMCParticle;}\r
+\r
+  /*\r
+   *Get the positive MC TParticle from the stack \r
+   */\r
+  TParticle * GetPositiveMCParticle() const{return fPositiveMCParticle;}\r
+\r
+  /*\r
+   *Get the mother MC TParticle from the stack \r
+   */\r
+  TParticle * GetMotherMCParticle() const{return fMotherMCParticle;}\r
+\r
+  Bool_t HasSameMCMother();\r
+\r
+  /*\r
+   *Get the MC stack \r
+   */\r
+  AliStack* GetMCStack() const{return fMCStack;}\r
+\r
+  /*\r
+   *Get the magnetic field from the ESD event \r
+   */\r
+  Double_t GetMagneticField() const{return fESDEvent->GetMagneticField();}\r
+\r
+  /*\r
+   *Get the primary vertex from the esd event\r
+   */\r
+  const AliESDVertex *GetPrimaryVertex() const {return fESDEvent->GetPrimaryVertex();}\r
+\r
+  /*\r
+   * Set the PID of the negative track\r
+   */\r
+  void SetNegativeTrackPID(Int_t negTrackPID){fNegativeTrackPID=negTrackPID;}\r
+\r
+  /*\r
+   * Set the PID of the positive track\r
+   */\r
+  void SetPositiveTrackPID(Int_t posTrackPID){fPositiveTrackPID=posTrackPID;}\r
+\r
+  /*\r
+   * Set the flag to use the kfparticle class. Will also disable the use of esd tracks\r
+   */\r
+  void UseKFParticle(){fUseKFParticle = kTRUE; fUseESDTrack = kFALSE;}\r
+\r
+  /*\r
+   *  Set the flag to use the esd track class. Will also disable the use of kf particles\r
+   */\r
+  void UseESDTrack(){fUseESDTrack = kTRUE; fUseKFParticle = kFALSE;}\r
+\r
+  /*\r
+   *  Set the flag to use improved vertex or not\r
+   */\r
+  void SetUseImprovedVertex(Bool_t useImprovedVertex){fUseImprovedVertex=useImprovedVertex;}\r
+\r
+  /*\r
+   * Return the number in the species array belonging to the negative or positive track pid.\r
+   */\r
+  Int_t GetSpeciesIndex(Int_t chargeOfTrack);\r
+\r
+  /*\r
+   * Return the x coordinate of the v0\r
+   */\r
+  Double_t GetX() const{return fCurrentXValue;}\r
+\r
+  /*\r
+   * Return the y coordinate of the v0\r
+   */\r
+  Double_t GetY() const{return fCurrentYValue;}\r
+\r
+  /*\r
+   * Return the Z coordinate of the v0\r
+   */\r
+  Double_t GetZ() const{return fCurrentZValue;}\r
+\r
+  /*\r
+   * Return the radius of the v0\r
+   */\r
+  Double_t GetXYRadius() const{return sqrt((Double_t)(fCurrentXValue*fCurrentXValue + fCurrentYValue*fCurrentYValue));}\r
+\r
+  /*\r
+   * Get the opening angle between the two tracks\r
+   */\r
+  Double_t GetOpeningAngle(){return fNegativeTrackLorentzVector->Angle(fPositiveTrackLorentzVector->Vect());}\r
+\r
+  Double_t GetNegativeTrackEnergy() const{return fCurrentNegativeKFParticle->E();}\r
+  Double_t GetPositiveTrackEnergy() const{return fCurrentPositiveKFParticle->E();}\r
+  Double_t GetMotherCandidateEnergy() const{return fCurrentMotherKFCandidate->E();}\r
+\r
+  Double_t GetNegativeTrackPt() const{return fNegativeTrackLorentzVector->Pt();}\r
+  Double_t GetPositiveTrackPt() const{return fPositiveTrackLorentzVector->Pt();}\r
+  Double_t GetMotherCandidatePt() const{return fMotherCandidateLorentzVector->Pt();}\r
+\r
+  Double_t GetNegativeTrackEta() const{return fNegativeTrackLorentzVector->Eta();}\r
+  Double_t GetPositiveTrackEta() const{return fPositiveTrackLorentzVector->Eta();}\r
+  Double_t GetMotherCandidateEta() const{return fMotherCandidateLorentzVector->Eta();}\r
+\r
+  Double_t GetMotherCandidateNDF() const{return fCurrentMotherKFCandidate->GetNDF();}\r
+  Double_t GetMotherCandidateChi2() const{return fCurrentMotherKFCandidate->GetChi2();}\r
+  Double_t GetMotherCandidateMass() const{return fMotherCandidateKFMass;}\r
+  Double_t GetMotherCandidateWidth() const{return fMotherCandidateKFWidth;}\r
+\r
+  Double_t GetNegativeTrackPhi() const;\r
+  Double_t GetPositiveTrackPhi() const;\r
+  Double_t GetMotherCandidatePhi() const;\r
+\r
+  void UpdateEventByEventData();\r
+  \r
+  Double_t GetMaxRCut() const{return fMaxR;}\r
+  Double_t GetEtaCut() const{return fEtaCut;}\r
+  Double_t GetPtCut() const{return fPtCut;}\r
+  Double_t GetChi2CutConversion() const{return fChi2CutConversion;}\r
+  Double_t GetChi2CutMeson() const{return fChi2CutMeson;}\r
+\r
+  void SetMaxRCut(Double_t maxR){fMaxR=maxR;}\r
+  void SetEtaCut(Double_t etaCut){fEtaCut=etaCut;}\r
+  void SetPtCut(Double_t ptCut){fPtCut=ptCut;}\r
+  void SetChi2CutConversion(Double_t chi2){fChi2CutConversion=chi2;}\r
+  void SetChi2CutMeson(Double_t chi2){fChi2CutMeson=chi2;}\r
+  \r
+  void SetXVertexCut(Double_t xVtx){fCurrentXValue=xVtx;}\r
+  void SetYVertexCut(Double_t yVtx){fCurrentYValue=yVtx;}\r
+  void SetZVertexCut(Double_t zVtx){fCurrentZValue=zVtx;}\r
+  void SetPIDProbability(Double_t pidProb){fPIDProbabilityCutPositiveParticle=pidProb; fPIDProbabilityCutNegativeParticle=pidProb;}\r
+  void SetPIDProbabilityNegativeParticle(Double_t pidProb){fPIDProbabilityCutNegativeParticle=pidProb;}\r
+  void SetPIDProbabilityPositiveParticle(Double_t pidProb){fPIDProbabilityCutPositiveParticle=pidProb;}\r
+  void SetSigmaMass(Double_t sigmaMass){fNSigmaMass=sigmaMass;}\r
+\r
+  void SetDoMCTruth(Bool_t doMC){fDoMC = doMC;}\r
+\r
+  void UpdateV0Information();\r
+\r
+  void SetHistograms(AliGammaConversionHistograms *histograms){fHistograms=histograms;}\r
+\r
+  vector<AliKFParticle> GetCurrentEventGoodV0s() const{return fCurrentEventGoodV0s;}\r
+  vector<AliKFParticle> GetPreviousEventGoodV0s() const{return fPreviousEventGoodV0s;}\r
+\r
+ private:\r
+  AliStack * fMCStack;           // pointer to MonteCarlo particle stack \r
+  AliMCEventHandler* fMCTruth;   // pointer to the MC event handler\r
+  TChain * fChain;               // pointer to the TChain\r
+  \r
+  AliESDInputHandler* fESDHandler;      //! pointer to esd object\r
+  AliESDEvent *fESDEvent;               //! pointer to esd object\r
+\r
+  AliGammaConversionHistograms *fHistograms;\r
+  \r
+  Int_t fCurrentV0IndexNumber;\r
+  AliESDv0 * fCurrentV0;                //! pointer to the current v0\r
+  AliKFParticle * fCurrentNegativeKFParticle;  //! pointer to the negative KF particle\r
+  AliKFParticle * fCurrentPositiveKFParticle;  //! pointer to the positive KF particle\r
+  AliKFParticle * fCurrentMotherKFCandidate;   //! pointer to the positive KF particle\r
+\r
+  AliESDtrack * fCurrentNegativeESDTrack;      //! pointer to the negative ESD track\r
+  AliESDtrack * fCurrentPositiveESDTrack;      //! pointer to the positive ESD track\r
\r
+  TLorentzVector * fNegativeTrackLorentzVector; //! pointer to the negative Track Lorentz Vector\r
+  TLorentzVector * fPositiveTrackLorentzVector; //! pointer to the positive Track Lorentz Vector\r
+  TLorentzVector * fMotherCandidateLorentzVector;   //! pointer to the mother candidate Track Lorentz Vector\r
+\r
+  Double_t fCurrentXValue;\r
+  Double_t fCurrentYValue;\r
+  Double_t fCurrentZValue;\r
+\r
+  Int_t fPositiveTrackPID;\r
+  Int_t fNegativeTrackPID;\r
+\r
+  TParticle *fNegativeMCParticle;      //!\r
+  TParticle *fPositiveMCParticle;      //!\r
+  TParticle *fMotherMCParticle;        //!\r
+\r
+  Double_t fMotherCandidateKFMass;\r
+  Double_t fMotherCandidateKFWidth;\r
+\r
+  Bool_t fUseKFParticle;\r
+  Bool_t fUseESDTrack;\r
+  Bool_t fDoMC;\r
+\r
+  //cuts\r
+  Double_t fMaxR;\r
+  Double_t fEtaCut;\r
+  Double_t fPtCut;\r
+  Double_t fChi2CutConversion;\r
+  Double_t fChi2CutMeson;\r
+  Double_t fPIDProbabilityCutNegativeParticle;\r
+  Double_t fPIDProbabilityCutPositiveParticle;\r
+  Double_t fXVertexCut;\r
+  Double_t fYVertexCut;\r
+  Double_t fZVertexCut;\r
+  \r
+  Double_t fNSigmaMass;\r
+  \r
+  Bool_t fUseImprovedVertex;\r
+  \r
+  vector<AliKFParticle> fCurrentEventGoodV0s;\r
+  vector<AliKFParticle> fPreviousEventGoodV0s;\r
+\r
+  ClassDef(AliV0Reader,0)\r
+};\r
+\r
+\r
+#endif\r
+\r
+\r
+\r
diff --git a/PWG4/PWG4GammaConvLinkDef.h b/PWG4/PWG4GammaConvLinkDef.h
new file mode 100644 (file)
index 0000000..22d160a
--- /dev/null
@@ -0,0 +1,11 @@
+#ifdef __CINT__
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliAnalysisTaskGammaConversion+;
+#pragma link C++ class AliV0Reader+;
+#pragma link C++ class AliGammaConversionHistograms+;
+
+#endif
diff --git a/PWG4/libPWG4GammaConv.pkg b/PWG4/libPWG4GammaConv.pkg
new file mode 100644 (file)
index 0000000..0141f57
--- /dev/null
@@ -0,0 +1,15 @@
+#-*- Mode: Makefile -*-
+
+SRCS = GammaConv/AliV0Reader.cxx GammaConv/AliAnalysisTaskGammaConversion.cxx GammaConv/AliGammaConversionHistograms.cxx
+      
+HDRS:= $(SRCS:.cxx=.h) 
+
+DHDR:= PWG4GammaConvLinkDef.h
+
+EXPORT:=$(SRCS:.cxx=.h)
+
+ifeq (win32gcc,$(ALICE_TARGET))
+PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) -lSTEERBase \
+                         -lESD -lAOD -lANALYSIS -lANALYSISalice \
+                         -L$(shell root-config --libdir) -lEG
+endif