task to do the isolation and correlation analysis only at generator level with the...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Jul 2012 17:18:39 +0000 (17:18 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Jul 2012 17:18:39 +0000 (17:18 +0000)
PWGGA/CMakelibPWGGACaloTrackCorrelations.pkg
PWGGA/CaloTrackCorrelations/AliAnaGeneratorKine.cxx [new file with mode: 0755]
PWGGA/CaloTrackCorrelations/AliAnaGeneratorKine.h [new file with mode: 0755]
PWGGA/CaloTrackCorrelations/macros/AddTaskGenKine.C [new file with mode: 0644]
PWGGA/CaloTrackCorrelations/macros/anaGenKine.C [new file with mode: 0644]
PWGGA/PWGGACaloTrackCorrelationsLinkDef.h

index b8c42c1..26e286e 100644 (file)
@@ -42,6 +42,7 @@ set ( SRCS
     CaloTrackCorrelations/AliAnaPhotonConvInCalo.cxx 
     CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx
     CaloTrackCorrelations/AliAnaRandomTrigger.cxx
+    CaloTrackCorrelations/AliAnaGeneratorKine.cxx
     )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
diff --git a/PWGGA/CaloTrackCorrelations/AliAnaGeneratorKine.cxx b/PWGGA/CaloTrackCorrelations/AliAnaGeneratorKine.cxx
new file mode 100755 (executable)
index 0000000..28dd79e
--- /dev/null
@@ -0,0 +1,859 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//_________________________________________________________________________
+// Do photon/pi0 analysis for isolation and correlation
+// at the generator level. Only for kine stack (ESDs)
+//
+//
+// -- Author: Gustavo Conesa (LPSC-CNRS-Grenoble) 
+//////////////////////////////////////////////////////////////////////////////
+
+// --- ROOT system ---
+#include "TH2F.h"
+#include "TParticle.h"
+#include "TDatabasePDG.h"
+
+//---- ANALYSIS system ----
+#include "AliAnaGeneratorKine.h" 
+#include "AliStack.h"  
+#include "AliGenPythiaEventHeader.h"
+
+ClassImp(AliAnaGeneratorKine)
+
+
+//__________________________________________
+AliAnaGeneratorKine::AliAnaGeneratorKine() : 
+AliAnaCaloTrackCorrBaseClass(), 
+fStack(0),
+fParton2(0),         fParton3(0), 
+fParton6(0),         fParton7(0),   
+fJet6(),             fJet7(),
+fPtHard(0),
+fhPtHard(0),         fhPtParton(0),    fhPtJet(0),
+fhPtPartonPtHard(0), fhPtJetPtHard(0), fhPtJetPtParton(0),
+fhPtPhoton(0),       fhPtPi0(0)
+{
+  //Default Ctor
+  
+  //Initialize parameters
+  InitParameters();
+  
+  for(Int_t i = 0; i < 4; i++)
+  {
+    fhPtPhotonLeading[i]              = fhPtPi0Leading[i]              = 0; 
+    fhPtPhotonLeadingIsolated[i]      = fhPtPi0LeadingIsolated[i]      = 0; 
+    fhZHardPhotonLeading[i]           = fhZHardPi0Leading[i]           = 0;            
+    fhZHardPhotonLeadingIsolated[i]   = fhZHardPi0LeadingIsolated[i]   = 0; 
+    fhZPartonPhotonLeading[i]         = fhZPartonPi0Leading[i]         = 0;            
+    fhZPartonPhotonLeadingIsolated[i] = fhZPartonPi0LeadingIsolated[i] = 0; 
+    fhZJetPhotonLeading[i]            = fhZJetPi0Leading[i]            = 0;            
+    fhZJetPhotonLeadingIsolated[i]    = fhZJetPi0LeadingIsolated[i]    = 0; 
+    fhXEPhotonLeading[i]              = fhXEPi0Leading[i]              = 0;            
+    fhXEPhotonLeadingIsolated[i]      = fhXEPi0LeadingIsolated[i]      = 0; 
+    fhXEUEPhotonLeading[i]            = fhXEUEPi0Leading[i]            = 0;            
+    fhXEUEPhotonLeadingIsolated[i]    = fhXEUEPi0LeadingIsolated[i]    = 0; 
+  }
+  
+}
+
+//_______________________________________________________________________________
+void  AliAnaGeneratorKine::CorrelateWithPartonOrJet(const TLorentzVector trigger, 
+                                                    const Int_t   indexTrig,                     
+                                                    const Int_t   pdgTrig, 
+                                                    const Bool_t  leading[4],
+                                                    const Bool_t  isolated[4], 
+                                                    Int_t & iparton )  
+{
+  //Correlate trigger with partons or jets, get z
+  
+  //Get the index of the mother
+  iparton =  (fStack->Particle(indexTrig))->GetFirstMother();
+  TParticle * mother = fStack->Particle(iparton);
+  while (iparton > 7) 
+  {
+    iparton   = mother->GetFirstMother();
+    mother = fStack->Particle(iparton);
+  }
+  
+  //printf("Mother is parton %d with pdg %d\n",imom,mother->GetPdgCode());
+  
+  if(iparton < 6)
+  {
+    //printf("This particle is not from hard process - pdg %d, parton index %d\n",pdgTrig, iparton);
+    return; 
+  }
+  
+  Float_t ptTrig   = trigger.Pt(); 
+  Float_t partonPt = fParton6->Pt();
+  Float_t jetPt    = fJet6.Pt();
+  if(iparton==7)
+  {
+    partonPt = fParton6->Pt();
+    jetPt    = fJet6.Pt();
+  }
+  
+  
+  fhPtPartonPtHard->Fill(fPtHard, partonPt/fPtHard);
+  fhPtJetPtHard   ->Fill(fPtHard, jetPt/fPtHard);
+  fhPtJetPtParton ->Fill(fPtHard, jetPt/partonPt);
+
+  Float_t zHard = ptTrig / fPtHard;
+  Float_t zPart = ptTrig / partonPt;
+  Float_t zJet  = ptTrig / jetPt;
+
+  //if(zHard > 1 ) printf("*** Particle energy larger than pT hard z=%f\n",zHard); 
+  
+  //printf("Z : hard %2.2f, parton %2.2f, jet %2.2f\n",zHard,zPart,zJet);
+  
+  for( Int_t i = 0; i < 4; i++ )
+  {
+    if(pdgTrig==111)
+    {
+      if(leading[i])
+      { 
+        fhZHardPi0Leading[i]  ->Fill(ptTrig,zHard);
+        fhZPartonPi0Leading[i]->Fill(ptTrig,zPart);
+        fhZJetPi0Leading[i]   ->Fill(ptTrig,zJet );
+        
+        if(isolated[i])
+        {
+          fhZHardPi0LeadingIsolated[i]  ->Fill(ptTrig,zHard);
+          fhZPartonPi0LeadingIsolated[i]->Fill(ptTrig,zPart);
+          fhZJetPi0LeadingIsolated[i]   ->Fill(ptTrig,zJet);
+        }
+      }
+    }// pi0
+    else if(pdgTrig==22)
+    {
+      if(leading[i])
+      { 
+        fhZHardPhotonLeading[i]  ->Fill(ptTrig,zHard);
+        fhZPartonPhotonLeading[i]->Fill(ptTrig,zPart);
+        fhZJetPhotonLeading[i]   ->Fill(ptTrig,zJet );
+        
+        if(isolated[i])
+        {
+          fhZHardPhotonLeadingIsolated[i]  ->Fill(ptTrig,zHard);
+          fhZPartonPhotonLeadingIsolated[i]->Fill(ptTrig,zPart);
+          fhZJetPhotonLeadingIsolated[i]   ->Fill(ptTrig,zJet);
+        }
+      }      
+    } // photon
+  } // conditions loop  
+}  
+
+
+//____________________________________________________
+TList *  AliAnaGeneratorKine::GetCreateOutputObjects()
+{  
+  // Create histograms to be saved in output file 
+  
+  TList * outputContainer = new TList() ; 
+  outputContainer->SetName("GenKineHistos") ; 
+  
+  Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();  
+  Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  
+  Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin(); 
+  
+  fhPtHard  = new TH1F("hPtHard"," pt hard for selected triggers",nptbins,ptmin,ptmax); 
+  fhPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
+  outputContainer->Add(fhPtHard);
+  
+  fhPtParton  = new TH1F("hPtParton"," pt parton for selected triggers",nptbins,ptmin,ptmax); 
+  fhPtParton->SetXTitle("p_{T}^{parton} (GeV/c)");
+  outputContainer->Add(fhPtParton);
+  
+  fhPtJet  = new TH1F("hPtJet"," pt jet for selected triggers",nptbins,ptmin,ptmax); 
+  fhPtJet->SetXTitle("p_{T}^{jet} (GeV/c)");
+  outputContainer->Add(fhPtJet);
+  
+  fhPtPartonPtHard  = new TH2F("hPtPartonPtHard","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2); 
+  fhPtPartonPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
+  fhPtPartonPtHard->SetYTitle("p_{T}^{parton}/p_{T}^{hard}");
+  outputContainer->Add(fhPtPartonPtHard);
+  
+  fhPtJetPtHard  = new TH2F("hPtJetPtHard","jet pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2); 
+  fhPtJetPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
+  fhPtJetPtHard->SetYTitle("p_{T}^{jet}/p_{T}^{hard}");
+  outputContainer->Add(fhPtJetPtHard);
+  
+  fhPtJetPtParton  = new TH2F("hPtJetPtParton","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2); 
+  fhPtJetPtParton->SetXTitle("p_{T}^{hard} (GeV/c)");
+  fhPtJetPtParton->SetYTitle("p_{T}^{jet}/p_{T}^{parton}");
+  outputContainer->Add(fhPtJetPtParton);
+  
+  
+  fhPtPhoton  = new TH1F("hPtPhoton","Input Photon",nptbins,ptmin,ptmax); 
+  fhPtPhoton->SetXTitle("p_{T} (GeV/c)");
+  outputContainer->Add(fhPtPhoton);
+
+  fhPtPi0  = new TH1F("hPtPi0","Input Pi0",nptbins,ptmin,ptmax); 
+  fhPtPi0->SetXTitle("p_{T} (GeV/c)");
+  outputContainer->Add(fhPtPi0);
+  
+  TString name[]  = {"","_EMC","_Photon","_EMC_Photon"};
+  TString title[] = {"",", neutral in EMCal",", neutral only photon like",", neutral in EMCal and only photon like"};
+
+  for(Int_t i = 0; i < 4; i++)
+  {
+    
+    // Pt
+    
+    fhPtPhotonLeading[i]  = new TH1F(Form("hPtPhotonLeading%s",name[i].Data()),
+                                  Form("Photon : Leading of all particles%s",title[i].Data()),
+                                  nptbins,ptmin,ptmax); 
+    fhPtPhotonLeading[i]->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhPtPhotonLeading[i]);
+    
+    fhPtPi0Leading[i]  = new TH1F(Form("hPtPi0Leading%s",name[i].Data()),
+                               Form("Pi0 : Leading of all particles%s",title[i].Data()),
+                               nptbins,ptmin,ptmax); 
+    fhPtPi0Leading[i]->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhPtPi0Leading[i]);  
+    
+    fhPtPhotonLeadingIsolated[i]  = new TH1F(Form("hPtPhotonLeadingIsolated%s",name[i].Data()),
+                                     Form("Photon : Leading of all particles%s, isolated",title[i].Data()),
+                                     nptbins,ptmin,ptmax); 
+    fhPtPhotonLeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhPtPhotonLeadingIsolated[i]);
+    
+    fhPtPi0LeadingIsolated[i]  = new TH1F(Form("hPtPi0LeadingIsolated%s",name[i].Data()),
+                                  Form("Pi0 : Leading of all particles%s, isolated",title[i].Data()),
+                                  nptbins,ptmin,ptmax); 
+    fhPtPi0LeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhPtPi0LeadingIsolated[i]);  
+    
+    
+    
+    // zHard
+    
+    fhZHardPhotonLeading[i]  = new TH2F(Form("hZHardPhotonLeading%s",name[i].Data()),
+                                     Form("Z-Hard of Photon : Leading of all particles%s",title[i].Data()),
+                                     nptbins,ptmin,ptmax,200,0,2); 
+    fhZHardPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZHardPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZHardPhotonLeading[i]);
+    
+    fhZHardPi0Leading[i]  = new TH2F(Form("hZHardPi0Leading%s",name[i].Data()),
+                                  Form("Z-Hard of Pi0 : Leading of all particles%s",title[i].Data()),
+                                  nptbins,ptmin,ptmax,200,0,2); 
+    fhZHardPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZHardPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZHardPi0Leading[i]);  
+    
+    fhZHardPhotonLeadingIsolated[i]  = new TH2F(Form("hZHardPhotonLeadingIsolated%s",name[i].Data()),
+                                             Form("Z-Hard of Photon : Leading of all particles%s, isolated",title[i].Data()),
+                                             nptbins,ptmin,ptmax,200,0,2); 
+    fhZHardPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZHardPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZHardPhotonLeadingIsolated[i]);
+    
+    fhZHardPi0LeadingIsolated[i]  = new TH2F(Form("hZHardPi0LeadingIsolated%s",name[i].Data()),
+                                          Form("Z-Hard of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
+                                          nptbins,ptmin,ptmax,200,0,2); 
+    fhZHardPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZHardPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZHardPi0LeadingIsolated[i]);  
+    
+    // zHard
+    
+    fhZPartonPhotonLeading[i]  = new TH2F(Form("hZPartonPhotonLeading%s",name[i].Data()),
+                                        Form("Z-Parton of Photon : Leading of all particles%s",title[i].Data()),
+                                        nptbins,ptmin,ptmax,200,0,2); 
+    fhZPartonPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZPartonPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZPartonPhotonLeading[i]);
+    
+    fhZPartonPi0Leading[i]  = new TH2F(Form("hZPartonPi0Leading%s",name[i].Data()),
+                                     Form("Z-Parton of Pi0 : Leading of all particles%s",title[i].Data()),
+                                     nptbins,ptmin,ptmax,200,0,2); 
+    fhZPartonPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZPartonPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZPartonPi0Leading[i]);  
+    
+    fhZPartonPhotonLeadingIsolated[i]  = new TH2F(Form("hZPartonPhotonLeadingIsolated%s",name[i].Data()),
+                                                Form("Z-Parton of Photon : Leading of all particles%s, isolated",title[i].Data()),
+                                                nptbins,ptmin,ptmax,200,0,2); 
+    fhZPartonPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZPartonPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZPartonPhotonLeadingIsolated[i]);
+    
+    fhZPartonPi0LeadingIsolated[i]  = new TH2F(Form("hZPartonPi0LeadingIsolated%s",name[i].Data()),
+                                             Form("Z-Parton of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
+                                             nptbins,ptmin,ptmax,200,0,2); 
+    fhZPartonPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZPartonPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZPartonPi0LeadingIsolated[i]);  
+    
+    
+    // zJet
+    
+    fhZJetPhotonLeading[i]  = new TH2F(Form("hZJetPhotonLeading%s",name[i].Data()),
+                                        Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
+                                        nptbins,ptmin,ptmax,200,0,2); 
+    fhZJetPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZJetPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZJetPhotonLeading[i]);
+    
+    fhZJetPi0Leading[i]  = new TH2F(Form("hZJetPi0Leading%s",name[i].Data()),
+                                     Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
+                                     nptbins,ptmin,ptmax,200,0,2); 
+    fhZJetPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZJetPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZJetPi0Leading[i]);  
+    
+    fhZJetPhotonLeadingIsolated[i]  = new TH2F(Form("hZJetPhotonLeadingIsolated%s",name[i].Data()),
+                                                Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
+                                                nptbins,ptmin,ptmax,200,0,2); 
+    fhZJetPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZJetPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZJetPhotonLeadingIsolated[i]);
+    
+    fhZJetPi0LeadingIsolated[i]  = new TH2F(Form("hZJetPi0LeadingIsolated%s",name[i].Data()),
+                                             Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
+                                             nptbins,ptmin,ptmax,200,0,2); 
+    fhZJetPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhZJetPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhZJetPi0LeadingIsolated[i]);    
+    
+    
+    // XE
+    
+    fhXEPhotonLeading[i]  = new TH2F(Form("hXEPhotonLeading%s",name[i].Data()),
+                                       Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
+                                       nptbins,ptmin,ptmax,200,0,2); 
+    fhXEPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhXEPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhXEPhotonLeading[i]);
+    
+    fhXEPi0Leading[i]  = new TH2F(Form("hXEPi0Leading%s",name[i].Data()),
+                                    Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
+                                    nptbins,ptmin,ptmax,200,0,2); 
+    fhXEPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhXEPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhXEPi0Leading[i]);  
+    
+    fhXEPhotonLeadingIsolated[i]  = new TH2F(Form("hXEPhotonLeadingIsolated%s",name[i].Data()),
+                                               Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
+                                               nptbins,ptmin,ptmax,200,0,2); 
+    fhXEPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhXEPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhXEPhotonLeadingIsolated[i]);
+    
+    fhXEPi0LeadingIsolated[i]  = new TH2F(Form("hXEPi0LeadingIsolated%s",name[i].Data()),
+                                            Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
+                                            nptbins,ptmin,ptmax,200,0,2); 
+    fhXEPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhXEPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhXEPi0LeadingIsolated[i]);          
+    
+    
+    // XE from UE
+    
+    fhXEUEPhotonLeading[i]  = new TH2F(Form("hXEUEPhotonLeading%s",name[i].Data()),
+                                       Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
+                                       nptbins,ptmin,ptmax,200,0,2); 
+    fhXEUEPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhXEUEPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhXEUEPhotonLeading[i]);
+    
+    fhXEUEPi0Leading[i]  = new TH2F(Form("hXEUEPi0Leading%s",name[i].Data()),
+                                    Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
+                                    nptbins,ptmin,ptmax,200,0,2); 
+    fhXEUEPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhXEUEPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhXEUEPi0Leading[i]);  
+    
+    fhXEUEPhotonLeadingIsolated[i]  = new TH2F(Form("hXEUEPhotonLeadingIsolated%s",name[i].Data()),
+                                               Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
+                                               nptbins,ptmin,ptmax,200,0,2); 
+    fhXEUEPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhXEUEPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhXEUEPhotonLeadingIsolated[i]);
+    
+    fhXEUEPi0LeadingIsolated[i]  = new TH2F(Form("hXEUEPi0LeadingIsolated%s",name[i].Data()),
+                                            Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
+                                            nptbins,ptmin,ptmax,200,0,2); 
+    fhXEUEPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
+    fhXEUEPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
+    outputContainer->Add(fhXEUEPi0LeadingIsolated[i]);          
+    
+  }
+  
+  return outputContainer;
+  
+}
+
+//____________________________________________
+void  AliAnaGeneratorKine::GetPartonsAndJets() 
+{
+  // Fill data members with partons,jets and generated pt hard 
+  
+  fStack =  GetMCStack() ;
+  
+  if(!fStack) 
+  {
+    printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - No Stack available, STOP\n");
+    abort();
+  }
+  
+  fParton2 = fStack->Particle(2) ;
+  fParton3 = fStack->Particle(3) ;
+  fParton6 = fStack->Particle(6) ;
+  fParton7 = fStack->Particle(7) ;
+  
+  Float_t p6phi = fParton6->Phi();
+  if(p6phi < 0) p6phi +=TMath::TwoPi();
+  Float_t p7phi = fParton7->Phi();
+  if(p7phi < 0) p7phi +=TMath::TwoPi();  
+  
+  //printf("parton6: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",fParton6->Pt(),fParton6->Eta(),p6phi, fParton6->GetPdgCode());
+  //printf("parton7: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",fParton7->Pt(),fParton7->Eta(),p7phi, fParton7->GetPdgCode());
+  
+  // Get the jets, only for pythia
+  if(!strcmp(GetReader()->GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
+  {
+    AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetReader()->GetGenEventHeader();
+    
+    fPtHard = pygeh->GetPtHard();
+    
+    //printf("pt Hard %2.2f\n",fPtHard);
+    
+    const Int_t nTriggerJets =  pygeh->NTriggerJets();
+        
+    Float_t tmpjet[]={0,0,0,0};
+    
+    // select the closest jet to parton
+    Float_t jet7R = 100;
+    Float_t jet6R = 100;
+    
+    for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
+    {
+      pygeh->TriggerJet(ijet, tmpjet);
+      
+      TLorentzVector jet(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
+      Float_t jphi = jet.Phi();
+      if(jphi < 0) jphi +=TMath::TwoPi();
+      
+      Double_t radius6 = GetIsolationCut()->Radius(fParton6->Eta(), p6phi, jet.Eta() , jphi) ;
+      Double_t radius7 = GetIsolationCut()->Radius(fParton7->Eta(), p7phi, jet.Eta() , jphi) ;
+      
+      //printf("jet %d: pt %2.2f, eta %2.2f, phi %2.2f, r6 %2.2f, r7 %2.2f\n",ijet,jet.Pt(),jet.Eta(),jphi,radius6, radius7);
+      
+      if (radius6 < jet6R)
+      {
+        jet6R = radius6;
+        fJet6 = jet;
+        
+      }
+      if (radius7 < jet7R) 
+      {
+        jet7R = radius7;
+        fJet7 = jet;
+      }
+            
+    } // jet loop
+    
+    //printf("jet6: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet6.Pt(),fJet6.Eta(),fJet6.Phi());
+    //printf("jet7: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet7.Pt(),fJet7.Eta(),fJet6.Phi());
+    
+  } // pythia header
+  
+  fhPtHard   ->Fill(fPtHard);
+  fhPtJet    ->Fill(fJet6.Pt());
+  fhPtJet    ->Fill(fJet7.Pt());
+  fhPtParton ->Fill(fParton6->Pt());
+  fhPtParton ->Fill(fParton7->Pt());
+
+}
+
+//___________________________________________________________
+void AliAnaGeneratorKine::GetXE(const TLorentzVector trigger,  
+                                const Int_t   indexTrig,                     
+                                const Int_t   pdgTrig, 
+                                const Bool_t  leading[4], 
+                                const Bool_t  isolated[4], 
+                                const Int_t   iparton)     
+{
+
+  // Calculate the real XE and the UE XE
+
+  Float_t ptThresTrack = 0.2;
+
+  Float_t ptTrig  = trigger.Pt();
+  Float_t etaTrig = trigger.Eta();
+  Float_t phiTrig = trigger.Phi();
+  if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
+  
+  //Loop on primaries, start from position 8, no partons
+  for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
+  {
+    
+    TParticle * particle = fStack->Particle(ipr) ;
+    
+    if(ipr==indexTrig) continue;
+    
+    Int_t   pdg    = particle->GetPdgCode();
+    Int_t   status = particle->GetStatusCode();
+        
+    // Compare trigger with final state particles
+    if( status != 1) continue ;
+    
+    Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+    
+    if(charge==0) continue; // construct xe only with charged        
+    
+    Float_t pt     = particle->Pt();
+    Float_t eta    = particle->Eta();
+    Float_t phi    = particle->Phi();
+    if(phi < 0 ) phi += TMath::TwoPi();
+    
+    if( pt < ptThresTrack)    continue ;
+    
+    if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
+    
+    //Isolation
+    Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
+    
+    Float_t xe = -pt/ptTrig*TMath::Cos(phi-phiTrig);
+    
+    //Get the index of the mother
+    Int_t ipartonAway =  particle->GetFirstMother();
+    TParticle * mother = fStack->Particle(ipartonAway);
+    while (ipartonAway > 7) 
+    {
+      ipartonAway   = mother->GetFirstMother();
+      mother = fStack->Particle(ipartonAway);
+    }
+    
+    if((ipartonAway==6 || ipartonAway==7) && iparton!=ipartonAway) 
+    {
+      //printf("xE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
+      if(radius > 1 ) continue; // avoid particles too far from trigger
+      
+      for( Int_t i = 0; i < 4; i++ )
+      {
+        if(pdgTrig==111)
+        {
+          if(leading[i])
+          { 
+            fhXEPi0Leading[i]  ->Fill(ptTrig,xe);
+            
+            if(isolated[i])
+            {
+              fhXEPi0LeadingIsolated[i]  ->Fill(ptTrig,xe);
+            }
+          }
+        }// pi0
+        else if(pdgTrig==22)
+        {
+          if(leading[i])
+          { 
+            fhXEPhotonLeading[i]  ->Fill(ptTrig,xe);
+            
+            if(isolated[i])
+            {
+              fhXEPhotonLeadingIsolated[i]  ->Fill(ptTrig,xe);
+            }
+          }      
+        } // photon
+      } // conditions loop  
+    } // Away side
+
+    if(ipartonAway!=6 && ipartonAway!=7) 
+    {
+      
+      //printf("xE UE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
+
+      for( Int_t i = 0; i < 4; i++ )
+      {
+        if(pdgTrig==111)
+        {
+          if(leading[i])
+          { 
+            fhXEUEPi0Leading[i]  ->Fill(ptTrig,xe);
+            
+            if(isolated[i])
+            {
+              fhXEUEPi0LeadingIsolated[i]  ->Fill(ptTrig,xe);
+            }
+          }
+        }// pi0
+        else if(pdgTrig==22)
+        {
+          if(leading[i])
+          { 
+            fhXEUEPhotonLeading[i]  ->Fill(ptTrig,xe);
+            
+            if(isolated[i])
+            {
+              fhXEUEPhotonLeadingIsolated[i]  ->Fill(ptTrig,xe);
+            }
+          }      
+        } // photon
+      } // conditions loop  
+    } // Away side    
+    
+  } // primary loop
+
+
+}
+
+//________________________________________
+void AliAnaGeneratorKine::InitParameters()
+{
+  
+  //Initialize the parameters of the analysis.
+  AddToHistogramsName("AnaGenKine_");
+  
+}
+
+//___________________________________________________________________________
+void  AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger, 
+                                                const Int_t indexTrig, 
+                                                const Int_t pdgTrig, 
+                                                Bool_t leading[4],
+                                                Bool_t isolated[4]) 
+{
+  // Check if the trigger is the leading particle
+  // In case of neutral particles check all neutral or neutral in EMCAL acceptance
+  
+  Float_t ptMaxCharged       = 0; // all charged
+  Float_t ptMaxNeutral       = 0; // all neutral
+  Float_t ptMaxNeutEMCAL     = 0; // for neutral, select them in EMCAL acceptance
+  Float_t ptMaxNeutPhot      = 0; // for neutral, take only photons
+  Float_t ptMaxNeutEMCALPhot = 0; // for neutral, take only photons in EMCAL acceptance 
+  
+  leading[0] = 0;
+  leading[1] = 0;
+  leading[2] = 0;
+  leading[3] = 0;
+  
+  isolated[0] = 0;
+  isolated[1] = 0;
+  isolated[2] = 0;
+  isolated[3] = 0;
+  
+  Float_t ptTrig  = trigger.Pt();
+  Float_t etaTrig = trigger.Eta();
+  Float_t phiTrig = trigger.Phi();
+  if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
+
+  // Minimum track or cluster energy
+  Float_t ptThresTrack = 0.2;
+  Float_t ptThresCalo  = 0.3;
+
+  //Isolation cuts
+  Float_t ptThresIC    = 0.5;
+  Float_t rThresIC     = 0.4;
+  Int_t   nICTrack     = 0;
+  Int_t   nICNeutral   = 0;
+  Int_t   nICNeutEMCAL = 0;
+  Int_t   nICNeutPhot  = 0;
+  Int_t   nICNeutEMCALPhot = 0;
+  
+  //Loop on primaries, start from position 8, no partons
+  for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
+  {
+    if(ipr == indexTrig) continue;
+    TParticle * particle = fStack->Particle(ipr) ;
+    
+    Int_t   imother= particle->GetFirstMother();
+    //printf("Leading ipr %d - mother %d\n",ipr, imother);
+    
+    if(imother==indexTrig)  continue ;
+    
+    Int_t   pdg    = particle->GetPdgCode();
+    Int_t   status = particle->GetStatusCode();
+     
+    // Compare trigger with final state particles
+    if( status != 1) continue ;
+    
+    Float_t pt     = particle->Pt();
+    Float_t eta    = particle->Eta();
+    Float_t phi    = particle->Phi();
+    if(phi < 0 ) phi += TMath::TwoPi();
+    
+    if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
+
+    Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+    
+    //Isolation
+    Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
+    
+    if(charge==0)
+    {
+      if(pt < ptThresCalo)  continue ;
+      
+      if( ptMaxNeutral < pt )                   ptMaxNeutral = pt;
+      if( pt > ptThresIC && radius < rThresIC ) nICNeutral++ ;
+
+      Bool_t phPDG = kFALSE;
+      if(pdg==22 || pdg==111) phPDG = kTRUE;
+    
+      //if(pt > ptTrig) printf(" --- pdg %d, phPDG %d pT %2.2f, pTtrig %2.2f, eta %2.2f, phi %2.2f\n",pdg,phPDG,pt,ptTrig,particle->Eta(), particle->Phi()*TMath::RadToDeg());
+      if(phPDG)
+      {
+        if( ptMaxNeutPhot < pt)                   ptMaxNeutPhot = pt;
+        if( pt > ptThresIC && radius < rThresIC ) nICNeutPhot++ ;
+      }
+      
+      //EMCAL acceptance
+      Bool_t inEMCAL = kTRUE;
+      if(TMath::Abs(particle->Eta()) > 0.7
+         || particle->Phi() > TMath::DegToRad()*180
+         || particle->Phi() < TMath::DegToRad()*80 ) inEMCAL = kFALSE ;
+      
+      if(inEMCAL)
+      {
+        if( ptMaxNeutEMCAL < pt )                 ptMaxNeutEMCAL = pt;
+        if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCAL++ ;
+
+        if(phPDG)
+        {
+          if( ptMaxNeutEMCALPhot < pt )             ptMaxNeutEMCALPhot = pt;
+          if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCALPhot++ ;
+        }
+      }
+    }
+    else  
+    {
+      if( pt < ptThresTrack)  continue ;
+
+      if( ptMaxCharged < pt )   ptMaxCharged   = pt;
+      
+      if( pt > ptThresIC && radius < rThresIC )  
+      {
+        //printf("UE track? pTtrig %2.2f, pt %2.2f, etaTrig %2.2f,  eta %2.2f, phiTrig %2.2f,  phi %2.2f, radius %2.2f\n",
+        //       ptTrig, pt,etaTrig, eta, phiTrig*TMath::RadToDeg(), phi*TMath::RadToDeg(),radius);
+        nICTrack++ ;
+      }
+    }
+
+  } // particle loop
+  
+  // Leding decision
+  if(ptTrig > ptMaxCharged)
+  {
+    //printf("pt charged %2.2f, pt neutral %2.2f, pt neutral emcal %2.2f, pt photon %2.2f, pt photon emcal %2.2f\n", 
+    //       ptMaxCharged, ptMaxNeutral, ptMaxNeutEMCAL,ptMaxNeutPhot, ptMaxNeutEMCALPhot);
+    if(ptTrig > ptMaxNeutral      ) leading[0] = kTRUE ;
+    if(ptTrig > ptMaxNeutEMCAL    ) leading[1] = kTRUE ;
+    if(ptTrig > ptMaxNeutPhot     ) leading[2] = kTRUE ;
+    if(ptTrig > ptMaxNeutEMCALPhot) leading[3] = kTRUE ;
+  }
+  
+  //printf("N in cone over threshold : tracks  %d, neutral %d, neutral emcal %d, photon %d, photon emcal %d\n", 
+  //       nICTrack, nICNeutral ,nICNeutEMCAL,nICNeutPhot, nICNeutEMCALPhot);
+  
+  // Isolation decision
+  if(nICTrack == 0)
+  {
+    if(nICNeutral       == 0 ) isolated[0] = kTRUE ;
+    if(nICNeutEMCAL     == 0 ) isolated[1] = kTRUE ;
+    if(nICNeutPhot      == 0 ) isolated[2] = kTRUE ;
+    if(nICNeutEMCALPhot == 0 ) isolated[3] = kTRUE ;
+  }
+  
+  // Fill histograms if conditions apply for all 4 cases
+  for( Int_t i = 0; i < 4; i++ )
+  {
+    if(pdgTrig==111)
+    {
+      if(leading[i])
+      { 
+        fhPtPi0Leading[i]->Fill(ptTrig);
+        if(isolated[i]) fhPtPi0LeadingIsolated[i]->Fill(ptTrig);
+      }
+    }// pi0
+    else if(pdgTrig==22)
+    {
+      if(leading[i])
+      { 
+        fhPtPhotonLeading[i]->Fill(ptTrig);
+        if(isolated[i]) fhPtPhotonLeadingIsolated[i]->Fill(ptTrig);
+      }      
+    } // photon
+  } // conditions loop
+    
+}
+  
+//_____________________________________________________
+void  AliAnaGeneratorKine::MakeAnalysisFillHistograms() 
+{
+  //Particle-Parton Correlation Analysis, fill histograms
+  
+  TLorentzVector trigger;
+  
+  GetPartonsAndJets();
+  
+  for(Int_t ipr = 0; ipr < fStack->GetNprimary(); ipr ++ )
+  {
+    TParticle * particle = fStack->Particle(ipr) ;
+    
+    Int_t   pdgTrig    = particle->GetPdgCode();
+    Int_t   statusTrig = particle->GetStatusCode();
+    Int_t   imother    = particle->GetFirstMother();
+    Float_t ptTrig     = particle->Pt(); 
+
+    // Select final state photons (prompt, fragmentation) or pi0s
+    
+    //Check the origin of the photon, accept if prompt or initial/final state radiation
+    if(pdgTrig == 22 && statusTrig == 1)
+    {
+      if(imother > 8) continue;
+    }
+    // If not photon, trigger on pi0
+    else if(pdgTrig != 111) continue;
+    
+    // Acceptance and kinematical cuts
+    if( ptTrig < 8 )    continue ;
+    
+    //EMCAL acceptance, a bit less
+    if(TMath::Abs(particle->Eta()) > 0.6) continue ;
+    if(particle->Phi() > TMath::DegToRad()*176) continue ;
+    if(particle->Phi() < TMath::DegToRad()*74 ) continue ;
+
+//    printf("Particle %d : pdg %d status %d, mother index %d, pT %2.2f, eta %2.2f, phi %2.2f \n",
+//           ipr, pdgTrig, statusTrig, imother, ptTrig, particle->Eta(), particle->Phi()*TMath::RadToDeg());
+    
+//    if(pdgTrig==111)
+//    {
+//      printf("\t pi0 daughters %d, %d\n", particle->GetDaughter(0), particle->GetDaughter(1));
+//    }
+    
+
+    if     (pdgTrig==22 ) fhPtPhoton->Fill(ptTrig);
+    else if(pdgTrig==111) fhPtPi0   ->Fill(ptTrig);
+    
+    // Check if it is leading
+    Bool_t leading[4] ;
+    Bool_t isolated[4] ;
+
+    particle->Momentum(trigger);
+    
+    IsLeadingAndIsolated(trigger, ipr, pdgTrig, leading, isolated);
+    
+    Int_t iparton = -1;
+    CorrelateWithPartonOrJet(trigger, ipr, pdgTrig, leading, isolated, iparton); 
+    
+    GetXE(trigger,ipr,pdgTrig,leading,isolated,iparton) ;    
+    
+  }
+  
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - End fill histograms \n");
+  
+} 
diff --git a/PWGGA/CaloTrackCorrelations/AliAnaGeneratorKine.h b/PWGGA/CaloTrackCorrelations/AliAnaGeneratorKine.h
new file mode 100755 (executable)
index 0000000..91ae308
--- /dev/null
@@ -0,0 +1,127 @@
+#ifndef ALIANAGENERATORKINE_H
+#define ALIANAGENERATORKINE_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+
+//___________________________________________________________________________
+// Do photon/pi0 analysis for isolation and correlation
+// at the generator level. Only for kine stack (ESDs)
+//
+//
+//-- Author: Gustavo Conesa (LPSC-CNRS-Grenoble)
+
+// --- ROOT ---
+class TH2F ;
+class TParticle ;
+class AliStack ;
+class TLorentzVector ;
+
+// --- ANALYSIS ---
+#include "AliAnaCaloTrackCorrBaseClass.h"
+
+class AliAnaGeneratorKine : public AliAnaCaloTrackCorrBaseClass {
+       
+public:
+  
+  AliAnaGeneratorKine() ; // default ctor
+  virtual ~AliAnaGeneratorKine() { ; } //virtual dtor              
+  
+  void CorrelateWithPartonOrJet(const TLorentzVector trigger,  
+                                const Int_t   indexTrig,                     
+                                const Int_t   pdgTrig, 
+                                const Bool_t  leading[4], 
+                                const Bool_t  isolated[4], 
+                                Int_t & iparton) ; 
+  
+    TList * GetCreateOutputObjects() ;
+  
+  void    GetPartonsAndJets() ;
+    
+  void GetXE(const TLorentzVector trigger,  
+             const Int_t   indexTrig,                     
+             const Int_t   pdgTrig, 
+             const Bool_t  leading[4], 
+             const Bool_t  isolated[4], 
+             const Int_t   iparton) ;    
+  
+  void    InitParameters() ;
+  
+  void    IsLeadingAndIsolated(const TLorentzVector trigger, 
+                               const Int_t   indexTrig,                     
+                               const Int_t   pdgTrig, 
+                               Bool_t leading[4],     
+                               Bool_t isolated[4]) ;
+  
+  void    MakeAnalysisFillAOD()  { ; }
+  
+  void    MakeAnalysisFillHistograms() ; 
+    
+private:
+  
+  AliStack  * fStack;                       //! access stack
+  
+  TParticle * fParton2;                     //! Initial state Parton
+  TParticle * fParton3;                     //! Initial state Parton
+  
+  TParticle * fParton6;                     //! Final state Parton
+  TParticle * fParton7;                     //! Final state Parton
+  
+  TLorentzVector fJet6;                     //! Pythia jet close to parton in position 6
+  TLorentzVector fJet7;                     //! Pythia jet close to parton in position 7
+
+  Float_t     fPtHard;                      //! Generated pT hard
+  
+  TH1F      * fhPtHard;                     //! pt of parton 
+  TH1F      * fhPtParton;                   //! pt of parton  
+  TH1F      * fhPtJet;                      //! pt of jet 
+  
+  TH2F      * fhPtPartonPtHard;             //! pt of parton divided to pt hard, trigger is photon 
+  TH2F      * fhPtJetPtHard;                //! pt of jet divided to pt hard, trigger is photon 
+  TH2F      * fhPtJetPtParton;              //! pt of parton divided to pt parton, trigger is photon 
+
+  TH1F      * fhPtPhoton;                   //! Input photon
+  TH1F      * fhPtPi0;                      //! Input pi0
+  
+  TH1F      * fhPtPhotonLeading[4];         //! Leading photon
+  TH1F      * fhPtPi0Leading[4];            //! Leading pi0
+  
+  TH1F      * fhPtPhotonLeadingIsolated[4]; //! Leading photon, isolated
+  TH1F      * fhPtPi0LeadingIsolated[4];    //! Leading pi0, isolated
+
+  TH2F      * fhZHardPhotonLeading[4];           //! Leading photon, zHard
+  TH2F      * fhZHardPi0Leading[4];              //! Leading pi0, zHard
+  TH2F      * fhZHardPhotonLeadingIsolated[4];   //! Leading photon, isolated, zHard
+  TH2F      * fhZHardPi0LeadingIsolated[4];      //! Leading pi0, isolated, zHard
+  
+  TH2F      * fhZPartonPhotonLeading[4];         //! Leading photon, zHard
+  TH2F      * fhZPartonPi0Leading[4];            //! Leading pi0, zHard
+  TH2F      * fhZPartonPhotonLeadingIsolated[4]; //! Leading photon, isolated, zHard
+  TH2F      * fhZPartonPi0LeadingIsolated[4];    //! Leading pi0, isolated, zHard
+
+  TH2F      * fhZJetPhotonLeading[4];            //! Leading photon, zHard
+  TH2F      * fhZJetPi0Leading[4];               //! Leading pi0, zHard
+  TH2F      * fhZJetPhotonLeadingIsolated[4];    //! Leading photon, isolated, zHard
+  TH2F      * fhZJetPi0LeadingIsolated[4];       //! Leading pi0, isolated, zHard
+  
+  TH2F      * fhXEPhotonLeading[4];              //! Leading photon, xE away side
+  TH2F      * fhXEPi0Leading[4];                 //! Leading pi0, xE away side
+  TH2F      * fhXEPhotonLeadingIsolated[4];      //! Leading photon, xE away side
+  TH2F      * fhXEPi0LeadingIsolated[4];         //! Leading pi0, isolated, xE away side
+  
+  TH2F      * fhXEUEPhotonLeading[4];              //! Leading photon, xE away side
+  TH2F      * fhXEUEPi0Leading[4];                 //! Leading pi0, xE away side
+  TH2F      * fhXEUEPhotonLeadingIsolated[4];      //! Leading photon, xE away side
+  TH2F      * fhXEUEPi0LeadingIsolated[4];         //! Leading pi0, isolated, xE away side
+  
+  AliAnaGeneratorKine              (const AliAnaGeneratorKine & gk) ; // cpy ctor
+  AliAnaGeneratorKine & operator = (const AliAnaGeneratorKine & gk) ; // cpy assignment
+  
+  ClassDef(AliAnaGeneratorKine,1)
+  
+} ;
+
+
+#endif //ALIANAGENERATORKINE_H
+
+
+
diff --git a/PWGGA/CaloTrackCorrelations/macros/AddTaskGenKine.C b/PWGGA/CaloTrackCorrelations/macros/AddTaskGenKine.C
new file mode 100644 (file)
index 0000000..a3e4c06
--- /dev/null
@@ -0,0 +1,249 @@
+
+TString  kCalorimeter = "EMCAL";
+Int_t    kYears       = 2011;
+Int_t    kDebug       = -1;
+Bool_t   kPrint       = kFALSE;
+
+AliAnalysisTaskCaloTrackCorrelation *AddTaskGenKine(TString outputfile, const Double_t scaleFactor   = -1)
+{
+  // Creates a CaloTrackCorr task, configures it and adds it to the analysis manager.
+    
+  // Get the pointer to the existing analysis manager via the static access method.
+  
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) 
+  {
+    ::Error("AddTask", "No analysis manager to connect to.");
+    return NULL;
+  }  
+    
+  // #### Configure analysis ####
+    
+  AliAnaCaloTrackCorrMaker * maker = new AliAnaCaloTrackCorrMaker();
+  printf("SCALE FACTOR %e\n",scaleFactor);
+  maker->SetScaleFactor(scaleFactor); // for MC, negative (not scaled) by default
+  
+  // General frame setting and configuration
+  maker->SetReader   (ConfigureReader()   ); 
+  maker->SetCaloUtils(ConfigureCaloUtils()); 
+
+  // Analysis tasks setting and configuration
+  Int_t n = 0;//Analysis number, order is important  
+  maker->AddAnalysis(ConfigureGenKine(), n++); // Photon cluster selection
+   
+  maker->SetAnaDebug(kDebug)  ;
+  maker->SwitchOnHistogramsMaker()  ;
+  maker->SwitchOffAODsMaker() ;
+  
+  if(kPrint) maker->Print("");
+  
+  // Create task
+  
+  AliAnalysisTaskCaloTrackCorrelation * task = new AliAnalysisTaskCaloTrackCorrelation ("GeneratorKineAnalysis");
+  task->SetConfigFileName(""); //Don't configure the analysis via configuration file.
+  task->SetDebugLevel(kDebug);
+  task->SetBranches("ESD:AliESDRun.,AliESDHeader"); 
+  task->SetAnalysisMaker(maker);
+  mgr->AddTask(task);
+  
+  //Create containers
+
+  if(outputfile.Length()==0) outputfile = AliAnalysisManager::GetCommonFileName(); 
+  
+  AliAnalysisDataContainer *cout_pc   = mgr->CreateContainer("GenKine", TList::Class(), 
+                                                             AliAnalysisManager::kOutputContainer, 
+                                                             Form("%s",outputfile.Data()));
+       
+  AliAnalysisDataContainer *cout_cuts = mgr->CreateContainer("Param_GemKine", TList::Class(), 
+                                                             AliAnalysisManager::kParamContainer, 
+                                                             "AnalysisParameters.root");
+  
+  // Create ONLY the output containers for the data produced by the task.
+  // Get and connect other common input/output containers via the manager as below
+  //==============================================================================
+  mgr->ConnectInput  (task, 0, mgr->GetCommonInputContainer());
+  // AOD output slot will be used in a different way in future
+  mgr->ConnectOutput (task, 1, cout_pc);
+  mgr->ConnectOutput (task, 2, cout_cuts);
+      
+  return task;
+}
+
+//____________________________________
+AliCaloTrackReader * ConfigureReader()
+{
+  
+  AliCaloTrackReader * reader = new AliCaloTrackMCReader();
+  
+  reader->SetDebug(kDebug);//10 for lots of messages
+    
+//  reader->SetPtHardAndJetPtComparison(kTRUE);
+//  reader->SetPtHardAndJetPtFactor(4);
+//  
+//  reader->SetPtHardAndClusterPtComparison(kTRUE);
+//  reader->SetPtHardAndClusterPtFactor(1.);
+
+  
+  reader->SwitchOffWriteDeltaAOD()  ;
+  
+  reader->SwitchOnStack();          
+  reader->SwitchOffAODMCParticles(); 
+  
+  //------------------------
+  // Detector input filling
+  //------------------------
+  
+  //Min cluster/track E
+  reader->SetEMCALEMin(0.3); 
+  reader->SetEMCALEMax(1000); 
+  reader->SetPHOSEMin(0.3);
+  reader->SetPHOSEMax(1000);
+  reader->SetCTSPtMin(0.2);
+  reader->SetCTSPtMax(1000);
+
+  reader->SwitchOnFiducialCut();
+  //reader->GetFiducialCut()->SetSimpleCTSFiducialCut(0.8, 0, 360) ;
+
+  reader->SwitchOffCTS();
+  reader->SwitchOffEMCALCells();  
+  reader->SwitchOffEMCAL();
+  reader->SwitchOffPHOSCells();  
+  reader->SwitchOffPHOS();
+  
+  reader->SetZvertexCut(10.);                // Open cut
+  reader->SwitchOnPrimaryVertexSelection(); // and besides primary vertex
+
+  reader->SwitchOffEventSelection();         // remove pileup by default
+  reader->SwitchOffV0ANDSelection() ;        // and besides v0 AND
+      
+  reader->SetImportGeometryFromFile(kFALSE);
+  
+  if(kPrint)reader->Print("");
+  
+  return reader;
+  
+}
+
+
+//_______________________________________
+AliCalorimeterUtils* ConfigureCaloUtils()
+{
+  
+  AliCalorimeterUtils *cu = new AliCalorimeterUtils;
+  cu->SetDebug(kDebug);
+  
+  
+  cu->SwitchOffClusterPlot();
+  
+  cu->SwitchOffRecalculateClusterTrackMatching(); // Done in clusterization
+  
+  cu->SwitchOffBadChannelsRemoval() ;
+  
+  cu->SwitchOffLoadOwnEMCALGeometryMatrices();
+  
+  AliEMCALRecoUtils * recou = cu->GetEMCALRecoUtils();
+  
+  cu->SwitchOffRecalibration(); 
+  cu->SwitchOffTimeRecalibration();     
+  cu->SwitchOffRunDepCorrection(); 
+  cu->SwitchOffCorrectClusterLinearity();
+  cu->SwitchOffEMCALOADB() ; 
+  
+  // PHOS 
+  cu->SwitchOffLoadOwnPHOSGeometryMatrices();
+  
+  if(kPrint) cu->Print("");
+  
+  return cu;
+  
+}
+
+
+//_____________________________________
+AliAnaGeneratorKine* ConfigureGenKine()
+{
+  
+  AliAnaGeneratorKine *ana = new AliAnaGeneratorKine();
+  ana->SetDebug(kDebug); //10 for lots of messages
+  
+  // cluster selection cuts
+  
+  ana->SwitchOffFiducialCut();
+
+    
+  ana->AddToHistogramsName("AnaGenKine_");
+  SetHistoRangeAndNBins(ana->GetHistogramRanges()); // see method below
+    
+  if(kPrint)ana->Print("");
+  
+  return ana;
+  
+}
+
+//________________________________________________________
+void SetHistoRangeAndNBins (AliHistogramRanges* histoRanges)
+{
+  // Set common bins for all analysis and MC histograms filling
+  
+  histoRanges->SetHistoPtRangeAndNBins(0, 200, 400) ; // Energy and pt histograms
+  
+  if(kCalorimeter=="EMCAL")
+  {
+    if(kYears==2010)
+    {
+      histoRanges->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 122*TMath::DegToRad(), 78) ;
+      histoRanges->SetHistoXRangeAndNBins(-230,90,120); // QA
+      histoRanges->SetHistoYRangeAndNBins(370,450,40);  // QA
+    }
+    else if(kYears==2011)
+    {           
+      histoRanges->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 182*TMath::DegToRad(), 108) ;
+      histoRanges->SetHistoXRangeAndNBins(-600,90,200); // QA
+      histoRanges->SetHistoYRangeAndNBins(100,450,100); // QA
+    }
+    else
+    {
+      histoRanges->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 190*TMath::DegToRad(), 122) ;
+      histoRanges->SetHistoXRangeAndNBins(-100,90,200); // QA
+      histoRanges->SetHistoYRangeAndNBins(50,450,100);  // QA
+    }
+    
+    histoRanges->SetHistoEtaRangeAndNBins(-0.72, 0.72, 144) ;
+  }
+  else
+  {
+    histoRanges->SetHistoPhiRangeAndNBins(260*TMath::DegToRad(), 320*TMath::DegToRad(), 60) ;
+    histoRanges->SetHistoEtaRangeAndNBins(-0.13, 0.13, 130) ;
+  }
+  
+  histoRanges->SetHistoShowerShapeRangeAndNBins(-0.1, 4.9, 500);
+  
+  // Invariant mass histoRangeslysis
+  histoRanges->SetHistoMassRangeAndNBins(0., 1., 200) ;
+  histoRanges->SetHistoAsymmetryRangeAndNBins(0., 1. , 100) ;
+  
+  // check if time calibration is on
+  histoRanges->SetHistoTimeRangeAndNBins(-1000.,1000,1000);
+  histoRanges->SetHistoDiffTimeRangeAndNBins(-200, 200, 800);
+  
+  // track-cluster residuals
+  histoRanges->SetHistoTrackResidualEtaRangeAndNBins(-0.15,0.15,300);
+  histoRanges->SetHistoTrackResidualPhiRangeAndNBins(-0.15,0.15,300);
+  histoRanges->SetHistodRRangeAndNBins(0.,0.15,150);//QA
+  
+  // QA, electron, charged
+  histoRanges->SetHistoPOverERangeAndNBins(0,10.,100);
+  histoRanges->SetHistodEdxRangeAndNBins(0.,200.,200);
+  
+  // QA
+  histoRanges->SetHistoFinePtRangeAndNBins(0, 10, 200) ; // bining for fhAmpId
+  histoRanges->SetHistoRatioRangeAndNBins(0.,2.,100);
+  histoRanges->SetHistoVertexDistRangeAndNBins(0.,500.,500);
+  histoRanges->SetHistoNClusterCellRangeAndNBins(0,500,500);
+  histoRanges->SetHistoZRangeAndNBins(-400,400,200);
+  histoRanges->SetHistoRRangeAndNBins(400,450,25);
+  histoRanges->SetHistoV0SignalRangeAndNBins(0,5000,500);
+  histoRanges->SetHistoV0MultiplicityRangeAndNBins(0,5000,500);
+  histoRanges->SetHistoTrackMultiplicityRangeAndNBins(0,5000,500);
+  
+}
diff --git a/PWGGA/CaloTrackCorrelations/macros/anaGenKine.C b/PWGGA/CaloTrackCorrelations/macros/anaGenKine.C
new file mode 100644 (file)
index 0000000..0a31bdf
--- /dev/null
@@ -0,0 +1,871 @@
+/* $Id:  $ */
+//--------------------------------------------------
+// Example macro to do analysis with the 
+// analysis classes in CaloTrackCorrelations
+// Can be executed with Root and AliRoot
+//
+// Pay attention to the options and definitions
+// set in the lines below
+//
+//  Author : Gustavo Conesa Balbastre (INFN-LNF)
+//
+//-------------------------------------------------
+enum anaModes {mLocal=0, mPROOF=1, mPlugin=2, mGRID=3};
+//mLocal    = 0: Analyze locally files in your computer
+//mPROOF    = 1: Analyze files on GRID with Plugin
+//mPlugin   = 2: Analyze files on GRID with Plugin
+//mGRID     = 3: Analyze files on GRID, jobs launched from aliensh
+
+//---------------------------------------------------------------------------
+// Settings to read locally several files, only for "mLocal" mode
+// The different values are default, they can be set with environmental 
+// variables: INDIR, PATTERN, NFILES, respectivelly
+
+char * kInDir   = "/user/data/files/"; 
+char * kPattern = ""; // Data are in files kInDir/kPattern+i 
+Int_t  kFile    = 6; 
+
+//---------------------------------------------------------------------------
+// Dataset for proof analysis, mode=mPROOF
+// char * kDataset = "/alice/vernet/PbPb_LHC10h_ESD";
+
+char *  kDatasetPROOF     = "/alice/vernet/LHC11b_149646";
+Int_t   kDatasetNMaxFiles = 20;
+TString ccin2p3UserName   = "arbor" ;
+TString alienUserName     = "narbor" ;
+
+//---------------------------------------------------------------------------
+// Collection file for grid analysis
+
+char * kXML = "collection.xml";
+
+//---------------------------------------------------------------------------
+//Scale histograms from file. Change to kTRUE when xsection file exists
+//Put name of file containing xsection 
+//Put number of events per ESD file
+//This is an specific case for normalization of Pythia files.
+const char * kXSFileName = "pyxsec.root";
+
+//---------------------------------------------------------------------------
+
+//Set some default values, but used values are set in the code!
+
+Bool_t  kMC        = kFALSE; //With real data kMC = kFALSE
+TString kInputData = "ESD"; //ESD, AOD, MC, deltaAOD
+Int_t   kYear      = 2011;
+TString kCollision = "pp";
+Bool_t  outAOD     = kFALSE; //Some tasks doesnt need it.
+TString kTreeName;
+TString kPass      = "";
+char    kTrigger[1024];
+Int_t   kRun       = 0;
+
+//________________________
+void anaGenKine(Int_t mode=mGRID)
+{
+  // Main
+  
+  //--------------------------------------------------------------------
+  // Load analysis libraries
+  // Look at the method below, 
+  // change whatever you need for your analysis case
+  // ------------------------------------------------------------------
+  
+  LoadLibraries(mode) ;
+  //gSystem->ListLibraries();
+  
+  //-------------------------------------------------------------------------------------------------
+  //Create chain from ESD and from cross sections files, look below for options.
+  //-------------------------------------------------------------------------------------------------
+  
+  // Set kInputData and kTreeName looking to the kINDIR
+  
+  //CheckInputData(mode);
+  
+  // Check global analysis settings  
+  
+  //CheckEnvironmentVariables();
+  
+  //printf("*********************************************\n");
+  //printf("*** Input data < %s >, pass %s, tree < %s >, MC?  < %d > ***\n",kInputData.Data(),kPass.Data(),kTreeName.Data(),kMC);
+  
+  //printf("*********************************************\n");
+  kTreeName  = "TE";
+  kInputData = "MC";
+  TChain * chain   = new TChain(kTreeName) ;
+  TChain * chainxs = new TChain("Xsection") ;
+  
+  chain  ->AddFile("galice.root");
+  chainxs->AddFile("pyxsec.root");
+
+  //CreateChain(mode, chain, chainxs); 
+  
+  Double_t scale = -1;
+  printf("===== kMC %d, chainxs %p\n",kMC,chainxs);
+  if(chainxs && chainxs->GetEntries() > 0)
+  {
+    Int_t nfiles = chainxs->GetEntries();
+    
+    //Get the cross section
+    Double_t xsection = 0; 
+    Float_t ntrials   = 0;
+    
+    GetAverageXsection(chainxs, xsection, ntrials);
+    
+    Int_t    nEventsPerFile = chain->GetEntries() / nfiles;
+    
+    Double_t trials = ntrials / nEventsPerFile ;
+    
+    scale = xsection/trials;
+    
+    printf("Get Cross section : nfiles  %d, nevents %d, nevents per file %d \n",nfiles, chain->GetEntries(),nEventsPerFile);     
+    printf("                    ntrials %d, trials %2.2f, xs %2.2e, scale factor %2.2e\n", ntrials,trials,xsection,scale);
+    
+  } 
+  
+  printf("*********************************************\n");
+  printf("number of entries # %lld, skipped %d\n", chain->GetEntries()) ;      
+  printf("*********************************************\n");
+  
+  if(!chain)
+  { 
+    printf("STOP, no chain available\n"); 
+    return;
+  }
+  
+  AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
+  
+  //------------------------------------------
+  //  Alien handler part
+  //------------------------------------------
+  AliAnalysisGrid *alienHandler=0x0;
+  if(mode==mPlugin)
+  {
+    // Create and configure the alien handler plugin
+    gROOT->LoadMacro("CreateAlienHandler.C");
+    alienHandler = CreateAlienHandler();
+    if (!alienHandler) return;
+  }  
+  
+  //--------------------------------------
+  // Make the analysis manager
+  //-------------------------------------
+  AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
+  //AliAnalysisManager::SetUseProgressBar(kTRUE);
+  //mgr->SetSkipTerminate(kTRUE);
+  //mgr->SetNSysInfo(1);
+  
+//  if(mode==mPlugin)
+//  {
+//    // Connect plugin to the analysis manager
+//    mgr->SetGridHandler(alienHandler);
+//  }
+  
+  // MC handler
+  if((kMC || kInputData == "MC") && !kInputData.Contains("AOD"))
+  {
+    AliMCEventHandler* mcHandler = new AliMCEventHandler();
+    mcHandler->SetReadTR(kFALSE);//Do not search TrackRef file
+    mgr->SetMCtruthEventHandler(mcHandler);
+    if( kInputData == "MC") 
+    {
+      cout<<"MC INPUT EVENT HANDLER"<<endl;
+      mgr->SetInputEventHandler(NULL);
+    }
+  }
+  
+  
+//  // AOD output handler
+//  if(kInputData!="deltaAOD" && outAOD)
+//  {
+//    cout<<"Init output handler"<<endl;
+//    AliAODHandler* aodoutHandler   = new AliAODHandler();
+//    aodoutHandler->SetOutputFileName("aod.root");
+//    ////aodoutHandler->SetCreateNonStandardAOD();
+//    mgr->SetOutputEventHandler(aodoutHandler);
+//  }
+//  
+//  //input
+//  
+//  if(kInputData == "ESD")
+//  {
+//    // ESD handler
+//    AliESDInputHandler *esdHandler = new AliESDInputHandler();
+//    esdHandler->SetReadFriends(kFALSE);
+//    mgr->SetInputEventHandler(esdHandler);
+//    cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
+//  }
+//  else if(kInputData.Contains("AOD"))
+//  {
+//    // AOD handler
+//    AliAODInputHandler *aodHandler = new AliAODInputHandler();
+//    mgr->SetInputEventHandler(aodHandler);
+//    if(kInputData == "deltaAOD") aodHandler->AddFriend("deltaAODCaloTrackCorr.root");
+//    cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
+//  }
+  
+  //mgr->RegisterExternalFile("deltaAODCaloTrackCorr.root");
+  mgr->SetDebugLevel(1); // For debugging, do not uncomment if you want no messages.
+  
+  TString outputFile = AliAnalysisManager::GetCommonFileName(); 
+   
+  gROOT->LoadMacro("AddTaskGenKine.C");  
+  
+  AliAnalysisTaskCaloTrackCorrelation *ana   = AddTaskGenKine(outputFile, scale);
+    
+  
+  //-----------------------
+  // Run the analysis
+  //-----------------------    
+  mgr->InitAnalysis();
+  mgr->PrintStatus();
+  
+  if      (mode == mPlugin) mgr->StartAnalysis("grid");
+  else if (mode == mPROOF ) mgr->StartAnalysis("proof",chain);
+  else                      mgr->StartAnalysis("local",chain);
+  
+  cout <<" Analysis ended sucessfully "<< endl ;
+  
+}
+
+//_____________________________
+void  LoadLibraries(Int_t mode)
+{
+  
+  if (mode == mPROOF) {
+    //TProof::Mgr("ccalpmaster")->SetROOTVersion("ALICE_v5-27-06b");
+    gROOT->LoadMacro("/afs/in2p3.fr/group/alice/laf/EnableAliRootForLAF.C");
+    TProof* proof = EnableAliRootForLAF("ccaplmaster",nPROOFWorkers.Data(),ccin2p3UserName.Data(),alienUserName.Data(),"",kFALSE,kTRUE,kTRUE,"OADB:ANALYSIS:ANALYSISalice:AOD:ESD:CORRFW:STEERBase:EMCALUtils:PHOSUtils:PWGCaloTrackCorrBase:PWGGACaloTrackCorrelations:PWGGAEMCALTasks");
+    
+    //  TProof* proof = TProof::Open("ccaplmaster",Form("workers=%s",nPROOFWorkers.Data()));
+    
+    //     //proof->ClearPackages();
+    //     proof->UploadPackage("STEERBase");
+    //     proof->UploadPackage("ESD");
+    //     proof->UploadPackage("AOD");
+    //     proof->UploadPackage("ANALYSIS");
+    //     proof->UploadPackage("OADB");
+    //     proof->UploadPackage("ANALYSISalice");
+    //     proof->UploadPackage("CORRFW");
+    //     //proof->UploadPackage("JETAN");
+    //     proof->UploadPackage("PHOSUtils");
+    //     proof->UploadPackage("EMCALUtils");
+    //     proof->UploadPackage("PWGCaloTrackCorrBase");
+    //     proof->UploadPackage("PWGGACaloTrackCorrelations");
+    //     proof->UploadPackage("PWGGAEMCALTasks");
+    
+    //     proof->EnablePackage("STEERBase");
+    //     proof->EnablePackage("ESD");
+    //     proof->EnablePackage("AOD");
+    //     proof->EnablePackage("ANALYSIS");
+    //     proof->EnablePackage("OADB");
+    //     proof->EnablePackage("ANALYSISalice");
+    //     proof->EnablePackage("CORRFW");
+    //     //proof->EnablePackage("JETAN");
+    //     proof->EnablePackage("PHOSUtils");
+    //     proof->EnablePackage("EMCALUtils");
+    //     proof->EnablePackage("PWGCaloTrackCorrBase");
+    //     proof->EnablePackage("PWGGACaloTrackCorrelations");
+    //     proof->EnablePackage("PWGGAEMCALTasks");
+    return;
+  }  
+  
+  //--------------------------------------
+  // Load the needed libraries most of them already loaded by aliroot
+  //--------------------------------------
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libXMLIO.so");
+  gSystem->Load("libMatrix.so");
+  gSystem->Load("libPhysics.so");
+  gSystem->Load("libMinuit.so"); // Root + libraries to if reclusterization is done
+  
+  gSystem->Load("libSTEERBase.so");
+  gSystem->Load("libGui.so"); // Root + libraries to if reclusterization is done
+  gSystem->Load("libCDB.so"); // Root + libraries to if reclusterization is done
+  gSystem->Load("libESD.so"); // Root + libraries to if reclusterization is done
+  gSystem->Load("libAOD.so");
+  gSystem->Load("libRAWDatabase.so"); // Root + libraries to if reclusterization is done
+  gSystem->Load("libProof.so"); 
+  gSystem->Load("libOADB");
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libSTEER.so"); // Root + libraries to if reclusterization is done
+  
+  gSystem->Load("libRAWDatarec.so"); // Root + libraries to if reclusterization is done
+  gSystem->Load("libRAWDatasim.so"); // Root + libraries to if reclusterization is done
+  gSystem->Load("libVZERObase.so");  // Root + libraries to if reclusterization is done
+  gSystem->Load("libVZEROrec.so");   // Root + libraries to if reclusterization is done
+  
+  gSystem->Load("libEMCALUtils");
+  //SetupPar("EMCALUtils");
+  gSystem->Load("libEMCALraw");  // Root + libraries to if reclusterization is done
+  gSystem->Load("libEMCALbase"); // Root + libraries to if reclusterization is done
+  gSystem->Load("libEMCALsim");  // Root + libraries to if reclusterization is done
+  gSystem->Load("libEMCALrec");  // Root + libraries to if reclusterization is done
+  //SetupPar("EMCALraw");
+  //SetupPar("EMCALbase");
+  //SetupPar("EMCALsim");
+  //SetupPar("EMCALrec");
+  
+  gSystem->Load("libANALYSISalice.so");
+  //gSystem->Load("libTENDER.so"); 
+  //gSystem->Load("libTENDERSupplies.so");
+  
+  gSystem->Load("libPHOSUtils");
+  gSystem->Load("libEMCALUtils");
+  gSystem->Load("libPWGCaloTrackCorrBase");
+  gSystem->Load("libPWGGACaloTrackCorrelations");
+  gSystem->Load("libPWGGAEMCALTasks");
+  //SetupPar("PWGCaloTrackCorrBase");
+  //SetupPar("PWGGACaloTrackCorrelations");
+  //SetupPar("PWGGAEMCALTasks");
+  
+  //gSystem->Load("libJETAN");
+  //gSystem->Load("FASTJETAN");
+  //gSystem->Load("PWGJE");
+
+  //gSystem->Load("libCORRFW.so");
+  //gSystem->Load("libPWGGAGammaConv.so"); 
+  //SetupPar("PWGGAGammaConv"); 
+  
+  // needed for plugin?
+  gSystem->AddIncludePath("-I$ALICE_ROOT");
+  gSystem->AddIncludePath("-I./");     
+  
+}
+
+//_________________________________
+void SetupPar(char* pararchivename)
+{
+  //Load par files, create analysis libraries
+  //For testing, if par file already decompressed and modified
+  //classes then do not decompress.
+  
+  TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
+  TString parpar(Form("%s.par", pararchivename)) ; 
+  
+  if ( gSystem->AccessPathName(pararchivename) ) {  
+    TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
+    gROOT->ProcessLine(processline.Data());
+  }
+  
+  TString ocwd = gSystem->WorkingDirectory();
+  gSystem->ChangeDirectory(pararchivename);
+  
+  // check for BUILD.sh and execute
+  if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
+    printf("*******************************\n");
+    printf("*** Building PAR archive    ***\n");
+    cout<<pararchivename<<endl;
+    printf("*******************************\n");
+    
+    if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
+      Error("runProcess","Cannot Build the PAR Archive! - Abort!");
+      return -1;
+    }
+  }
+  // check for SETUP.C and execute
+  if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
+    printf("*******************************\n");
+    printf("*** Setup PAR archive       ***\n");
+    cout<<pararchivename<<endl;
+    printf("*******************************\n");
+    gROOT->Macro("PROOF-INF/SETUP.C");
+  }
+  
+  gSystem->ChangeDirectory(ocwd.Data());
+  printf("Current dir: %s\n", ocwd.Data());
+}
+
+//______________________________________
+void CheckInputData(const anaModes mode)
+{
+  //Sets input data and tree
+  
+  TString ocwd = gSystem->WorkingDirectory();
+  
+  //---------------------------------------
+  //Local files analysis
+  //---------------------------------------
+  if(mode == mLocal){    
+    //If you want to add several ESD files sitting in a common directory INDIR
+    //Specify as environmental variables the directory (INDIR), the number of files 
+    //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
+    
+    if(gSystem->Getenv("INDIR"))  
+      kInDir = gSystem->Getenv("INDIR") ; 
+    else cout<<"INDIR not set, use default: "<<kInDir<<endl;   
+    
+    TString sindir(kInDir);
+    if     (sindir.Contains("pass1")) kPass = "pass1";
+    else if(sindir.Contains("pass2")) kPass = "pass2";
+    else if(sindir.Contains("pass3")) kPass = "pass3";
+    
+    if(gSystem->Getenv("PATTERN"))   
+      kPattern = gSystem->Getenv("PATTERN") ; 
+    else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
+    
+    cout<<"INDIR   : "<<kInDir<<endl;
+    cout<<"NFILES  : "<<kFile<<endl;
+    
+    char fileE[120] ;   
+    char fileA[120] ;   
+    char fileG[120] ;
+    char fileEm[120] ;   
+    for (Int_t event = 0 ; event < kFile ; event++) {
+      sprintf(fileE,  "%s/%s%d/AliESDs.root",    kInDir,kPattern,event) ; 
+      sprintf(fileA,  "%s/%s%d/AliAOD.root",     kInDir,kPattern,event) ; 
+      sprintf(fileG,  "%s/%s%d/galice.root",     kInDir,kPattern,event) ; 
+      sprintf(fileEm, "%s/%s%d/embededAOD.root", kInDir,kPattern,event) ; 
+      
+      TFile * fESD = TFile::Open(fileE) ; 
+      TFile * fAOD = TFile::Open(fileA) ; 
+      
+      //Check if file exists and add it, if not skip it
+      if (fESD) 
+      {
+        kTreeName  = "esdTree";
+        kInputData = "ESD";
+        TFile * fG = TFile::Open(fileG);
+        if(fG) { kMC = kTRUE; fG->Close();}
+        else     kMC = kFALSE;
+        
+        // Get run number
+        TTree* esdTree = (TTree*)fESD->Get("esdTree");
+        AliESDEvent* esd = new AliESDEvent();
+        esd->ReadFromTree(esdTree);
+        esdTree->GetEvent(0);
+        kRun = esd->GetRunNumber();
+        
+        return;
+      }
+      else if(fAOD)
+      {
+        kTreeName  = "aodTree";
+        kInputData = "AOD";
+        if(((TTree*) fAOD->Get("aodTree"))->GetBranch("mcparticles")) kMC=kTRUE;
+        else kMC = kFALSE;
+        
+        // Get run number
+        TTree* aodTree = (TTree*)fAOD->Get("aodTree");
+        AliAODEvent* aod = new AliAODEvent();
+        aod->ReadFromTree(aodTree);
+        aodTree->GetEvent(0);
+        kRun = aod->GetRunNumber();
+        return;
+      }
+      else if(TFile::Open(fileEm))
+      {
+        kTreeName  = "aodTree";
+        kInputData = "AOD";
+        kMC        = kTRUE;
+        
+        return;
+      }
+      else if(TFile::Open(fileG))
+      {
+        kTreeName  = "TE";
+        kInputData = "MC";
+        kMC        = kTRUE;
+        return;
+      }
+    }
+    
+    if(fESD) fESD->Close();
+    if(fAOD) fAOD->Close();
+    
+  }// local files analysis
+  
+  //------------------------------
+  //GRID xml files
+  //-----------------------------
+  else if(mode == mGRID){
+    //Get colection file. It is specified by the environmental
+    //variable XML
+    
+    if(gSystem->Getenv("XML") )
+      kXML = gSystem->Getenv("XML");
+    else
+      sprintf(kXML, "collection.xml") ; 
+    
+    if (!TFile::Open(kXML)) {
+      printf("No collection file with name -- %s -- was found\n",kXML);
+      return ;
+    }
+    else cout<<"XML file "<<kXML<<endl;
+    
+    //Load necessary libraries and connect to the GRID
+    gSystem->Load("libNetx.so") ; 
+    gSystem->Load("libRAliEn.so"); 
+    TGrid::Connect("alien://") ;
+    
+    //Feed Grid with collection file
+    TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
+    if (! collection) {
+      AliError(Form("%s not found", kXML)) ; 
+      return kFALSE ; 
+    }
+    TGridResult* result = collection->GetGridResult("",0 ,0);
+    
+    for (Int_t index = 0; index < result->GetEntries(); index++) {
+      TString alienURL = result->GetKey(index, "turl") ; 
+      cout << "================== " << alienURL << endl ; 
+      
+      if     (alienURL.Contains("pass1")) kPass = "pass1";
+      else if(alienURL.Contains("pass2")) kPass = "pass2";
+      else if(alienURL.Contains("pass3")) kPass = "pass3";
+      
+      kRun = AliAnalysisManager::GetRunFromAlienPath(alienURL.Data());
+      printf("Run number from alien path = %d\n",kRun);
+      
+      TFile * fAOD = 0 ; 
+      //Check if file exists and add it, if not skip it
+      if (alienURL.Contains("AliESDs.root"))  
+      {
+        kTreeName  = "esdTree";
+        kInputData = "ESD";
+        alienURL.ReplaceAll("AliESDs.root","galice.root");
+        if(TFile::Open(alienURL)) kMC=kTRUE;
+        else kMC = kFALSE;
+        return;
+      }
+      else if(alienURL.Contains("AliAOD.root"))
+      {
+        kTreeName  = "aodTree";
+        kInputData = "AOD";
+        fAOD = TFile::Open(alienURL);
+        if(((TTree*) fAOD->Get("aodTree"))->GetBranch("mcparticles")) kMC=kTRUE;
+        else kMC = kFALSE;
+        return;
+      }
+      else if(alienURL.Contains("embededAOD.root"))
+      {
+        kTreeName  = "aodTree";
+        kInputData = "AOD";
+        kMC=kTRUE;
+        return;
+      }
+      else if(alienURL.Contains("galice.root"))
+      {
+        kTreeName  = "TE";
+        kInputData = "MC";
+        kMC=kTRUE;
+        return;
+      } 
+    }
+  }// xml analysis
+  //------------------------------
+  //PROOF files
+  //-----------------------------
+  else if(mode == mPROOF){
+    
+    TFileCollection* coll  = gProof->GetDataSet(kDatasetPROOF)->GetStagedSubset();
+    
+    TIter iter(coll->GetList());
+    
+    TFileInfo* fileInfo = 0;
+    while ((fileInfo = dynamic_cast<TFileInfo*> (iter())))
+    {
+      if (fileInfo->GetFirstUrl()) {
+        TString ProofURL = fileInfo->GetFirstUrl()->GetUrl();
+        cout << "================== " << ProofURL << endl ; 
+        
+        if     (ProofURL.Contains("pass1")) kPass = "pass1";
+        else if(ProofURL.Contains("pass2")) kPass = "pass2";
+        else if(ProofURL.Contains("pass3")) kPass = "pass3";
+        
+        kRun = AliAnalysisManager::GetRunFromAlienPath(ProofURL.Data());
+        printf("Run number from alien path = %d\n",kRun);
+        
+        TFile * fAOD = 0 ; 
+        //Check if file exists and add it, if not skip it
+        if (ProofURL.Contains("AliESDs.root"))  
+        {
+          kTreeName  = "esdTree";
+          kInputData = "ESD";
+          alienURL.ReplaceAll("AliESDs.root","galice.root");
+          if(TFile::Open(ProofURL)) kMC=kTRUE;
+          else kMC = kFALSE;
+          
+          return;
+        }
+        else if(ProofURL.Contains("AliAOD.root"))
+        {
+          kTreeName  = "aodTree";
+          kInputData = "AOD";
+          fAOD = TFile::Open(ProofURL);
+          if(((TTree*) fAOD->Get("aodTree"))->GetBranch("mcparticles")) kMC=kTRUE;
+          else kMC = kFALSE;
+          return;
+        }
+        else if(ProofURL.Contains("embededAOD.root"))
+        {
+          kTreeName  = "aodTree";
+          kInputData = "AOD";
+          kMC=kTRUE;
+          return;
+        }
+        else if(ProofURL.Contains("galice.root"))
+        {
+          kTreeName  = "TE";
+          kInputData = "MC";
+          kMC=kTRUE;
+          return;
+        } 
+      }
+    }
+  }// proof analysis
+  
+  gSystem->ChangeDirectory(ocwd.Data());
+  
+}
+
+//_____________________________________________________________________
+void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs)
+{
+  //Fills chain with data
+  TString ocwd = gSystem->WorkingDirectory();
+  
+  //---------------------------------------
+  // Local files analysis
+  //---------------------------------------
+  if(mode == mLocal){    
+    //If you want to add several ESD files sitting in a common directory INDIR
+    //Specify as environmental variables the directory (INDIR), the number of files 
+    //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
+    
+    if(gSystem->Getenv("INDIR"))  
+      kInDir = gSystem->Getenv("INDIR") ; 
+    else cout<<"INDIR not set, use default: "<<kInDir<<endl;   
+    
+    if(gSystem->Getenv("PATTERN"))   
+      kPattern = gSystem->Getenv("PATTERN") ; 
+    else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
+    
+    if(gSystem->Getenv("NFILES"))
+      kFile = atoi(gSystem->Getenv("NFILES")) ;
+    else cout<<"NFILES not set, use default: "<<kFile<<endl;
+    
+    //Check if env variables are set and are correct
+    if ( kInDir  && kFile) {
+      printf("Get %d files from directory %s\n",kFile,kInDir);
+      if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
+        printf("%s does not exist\n", kInDir) ;
+        return ;
+      }
+      
+      //if(gSystem->Getenv("XSFILE"))  
+      //kXSFileName = gSystem->Getenv("XSFILE") ; 
+      //else cout<<" XS file name not set, use default: "<<kXSFileName<<endl;  
+      char * kGener = gSystem->Getenv("GENER");
+      if(kGener) {
+        cout<<"GENER "<<kGener<<endl;
+        if     (!strcmp(kGener,"PYTHIA")) kXSFileName = "pyxsec.root";
+        else if(!strcmp(kGener,"HERWIG")) kXSFileName = "hexsec.root";
+        else cout<<" UNKNOWN GENER, use default: "<<kXSFileName<<endl;
+      }
+      else cout<<" GENER not set, use default xs file name: "<<kXSFileName<<endl;
+      
+      cout<<"INDIR   : "<<kInDir     <<endl;
+      cout<<"NFILES  : "<<kFile      <<endl;
+      cout<<"PATTERN : "<<kPattern   <<endl;
+      cout<<"XSFILE  : "<<kXSFileName<<endl;
+      
+      TString datafile="";
+      if     (kInputData == "ESD")        datafile = "AliESDs.root" ;
+      else if(kInputData.Contains("AOD")) datafile = "AliAOD.root"  ;
+      else if(kInputData == "MC")         datafile = "galice.root"  ;
+      
+      //Loop on ESD/AOD/MC files, add them to chain
+      Int_t event =0;
+      Int_t skipped=0 ; 
+      char file[120] ;
+      char filexs[120] ;
+      
+      for (event = 0 ; event < kFile ; event++) {
+        sprintf(file,   "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ; 
+        sprintf(filexs, "%s/%s%d/%s", kInDir,kPattern,event,kXSFileName) ; 
+        TFile * fData = 0 ; 
+        //Check if file exists and add it, if not skip it
+        if ( fData = TFile::Open(file)) {
+          if ( fData->Get(kTreeName) ) { 
+            printf("++++ Adding %s\n", file) ;
+            chain->AddFile(file);
+            chainxs->Add(filexs) ; 
+          }
+        }
+        else { 
+          printf("---- Skipping %s\n", file) ;
+          skipped++ ;
+        }
+      }
+    }
+    else {
+      TString input = "AliESDs.root" ;
+      cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
+      chain->AddFile(input);
+    }
+    
+  }// local files analysis
+  
+  //------------------------------
+  // GRID xml files
+  //------------------------------
+  else if(mode == mGRID){
+    //Get colection file. It is specified by the environmental
+    //variable XML
+    
+    //Feed Grid with collection file
+    TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
+    if (! collection) {
+      AliError(Form("%s not found", kXML)) ; 
+      return kFALSE ; 
+    }
+    
+    TGridResult* result = collection->GetGridResult("",0 ,0);
+    
+    // Makes the ESD chain 
+    printf("*** Getting the Chain       ***\n");
+    for (Int_t index = 0; index < result->GetEntries(); index++) {
+      TString alienURL = result->GetKey(index, "turl") ; 
+      cout << "================== " << alienURL << endl ; 
+      chain->Add(alienURL) ; 
+      alienURL.ReplaceAll("AliESDs.root",kXSFileName);
+      chainxs->Add(alienURL) ; 
+    }
+  }// xml analysis
+  
+  //------------------------------
+  // PROOF
+  //------------------------------
+  else if (mode == mPROOF) {
+    
+    TFileCollection* ds= gProof->GetDataSet(kDatasetPROOF)->GetStagedSubset();
+    
+    gROOT->LoadMacro("/afs/in2p3.fr/group/alice/laf/dataset_management/CreateChainFromDataSet.C");
+    chain = CreateChainFromDataSet(ds, kTreeName , kDatasetNMaxFiles);
+    printf("chain has %d entries\n",chain->GetEntries());
+  }
+  
+  gSystem->ChangeDirectory(ocwd.Data());
+  
+}
+
+//______________________________
+void CheckEnvironmentVariables()
+{
+  
+  sprintf(kTrigger,"");
+  
+  Bool_t bRecalibrate = kFALSE;
+  Bool_t bBadChannel = kFALSE;
+  
+  for (int i=0; i< gApplication->Argc();i++){
+    
+#ifdef VERBOSEARGS
+    
+    printf("Arg  %d:  %s\n",i,gApplication->Argv(i));
+    
+#endif
+    
+    TString sRun = "";
+    
+    if (!(strcmp(gApplication->Argv(i),"--trigger")))
+      sprintf(trigger,gApplication->Argv(i+1));
+    
+    if (!(strcmp(gApplication->Argv(i),"--recalibrate")))
+      bRecalibrate = atoi(gApplication->Argv(i+1));
+    
+    if (!(strcmp(gApplication->Argv(i),"--badchannel")))
+      bBadChannel = atoi(gApplication->Argv(i+1));
+    
+    if (!(strcmp(gApplication->Argv(i),"--run"))){
+      sRun = gApplication->Argv(i+1);
+      if(sRun.Contains("LHC10")) {
+        kYear = 2010;
+      }
+      else {
+        if(kRun <=0){
+          kRun = atoi(gApplication->Argv(i+1));
+        }
+        else printf("** Run number already set  to %d, do not set to %d\n",kRun,atoi(gApplication->Argv(i+1)));
+      }//numeric run
+    }//--run available
+    
+  }// args loop
+  
+  if(!sRun.Contains("LHC10")){
+    if     ( kRun < 140000) 
+    {
+      kYear = 2010;
+      if( kRun >= 136851 ) kCollision = "PbPb";
+    }
+    else if( kRun < 170600)
+    {
+      kYear = 2011;
+      if( kRun >= 166500 ) kCollision = "PbPb";
+    }
+    else 
+    {
+      kYear = 2012;
+
+    }
+  }
+  
+  if(kMC) sprintf(kTrigger,"");
+  
+  printf("*********************************************\n");
+  //printf("*** Settings trigger %s, recalibrate %d, remove bad channels %d, year %d, collision %s, run %d ***\n",
+  //       kTrigger,bRecalibrate,bBadChannel, kYear,kCollision.Data(), kRun);
+  printf("*** Settings year %d, collision %s, run %d ***\n",kYear,kCollision.Data(), kRun);
+  printf("*********************************************\n");
+  
+}
+
+
+//_________________________________________________________________
+void GetAverageXsection(TTree * tree, Double_t & xs, Float_t & ntr)
+{
+  // Read the PYTHIA statistics from the file pyxsec.root created by
+  // the function WriteXsection():
+  // integrated cross section (xsection) and
+  // the  number of Pyevent() calls (ntrials)
+  // and calculate the weight per one event xsection/ntrials
+  // The spectrum calculated by a user should be
+  // multiplied by this weight, something like this:
+  // TH1F *userSpectrum ... // book and fill the spectrum
+  // userSpectrum->Scale(weight)
+  //
+  // Yuri Kharlov 19 June 2007
+  // Gustavo Conesa 15 April 2008
+  Double_t xsection = 0;
+  UInt_t    ntrials = 0;
+  xs = 0;
+  ntr = 0;
+  
+  Int_t nfiles =  tree->GetEntries()  ;
+  if (tree && nfiles > 0) {
+    tree->SetBranchAddress("xsection",&xsection);
+    tree->SetBranchAddress("ntrials" ,&ntrials );
+    for(Int_t i = 0; i < nfiles; i++){
+      tree->GetEntry(i);
+      xs  += xsection ;
+      ntr += ntrials ;
+      cout << "xsection " <<xsection<<" ntrials "<<ntrials<<endl; 
+    }
+    
+    xs =   xs /  nfiles;
+    ntr =  ntr / nfiles;
+    cout << "-----------------------------------------------------------------"<<endl;
+    cout << "Average of "<< nfiles<<" files: xsection " <<xs<<" ntrials "<<ntr<<endl; 
+    cout << "-----------------------------------------------------------------"<<endl;
+  } 
+  else cout << " >>>> Empty tree !!!! <<<<< "<<endl;
+  
+}
+
+
+
index bab25d2..766637b 100644 (file)
@@ -19,5 +19,6 @@
 #pragma link C++ class AliAnaPhotonConvInCalo+;
 #pragma link C++ class AliAnaInsideClusterInvariantMass+;
 #pragma link C++ class AliAnaRandomTrigger+;
+#pragma link C++ class AliAnaGeneratorKine+;
 
 #endif