{
// Load the ESD library
- gSystem->Load("libAnalysisCheck");
+ gSystem->Load("libGamma");
// Set the Include paths
- gSystem->SetIncludePath("-I$ROOTSYS/include -IAnalysisCheck");
- gROOT->ProcessLine(".include AnalysisCheck");
+ gSystem->SetIncludePath("-I$ROOTSYS/include -IPWG4");
+ gROOT->ProcessLine(".include PWG4");
// Set our location, so that other packages can find us
- gSystem->Setenv("AnalysisCheck_INCLUDE", "AnalysisCheck");
+ gSystem->Setenv("Gamma_INCLUDE", "PWG4");
}
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+// Class for the analysis of gamma correlations (gamma-jet,
+// gamma-hadron and isolation cut.
+// This class contains 3 main methods: one to fill lists of particles (ESDs) comming
+// from the CTS (ITS+TPC) and the calorimeters; the other one tags a candidate
+// cluster as isolated; the last one search in the
+// corresponing calorimeter for the highest energy cluster, identified it as
+// prompt photon;
+//
+// Class created from old AliPHOSGammaJet
+// (see AliRoot versions previous Release 4-09)
+//
+//*-- Author: Gustavo Conesa (LNF-INFN)
+//////////////////////////////////////////////////////////////////////////////
+
+
+// --- ROOT system ---
+
+#include <TFile.h>
+#include <TParticle.h>
+#include <TH2.h>
+
+#include "AliAnaGammaDirect.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliESDCaloCluster.h"
+#include "Riostream.h"
+#include "AliLog.h"
+
+ClassImp(AliAnaGammaDirect)
+
+//____________________________________________________________________________
+ AliAnaGammaDirect::AliAnaGammaDirect(const char *name) :
+ AliAnalysisTask(name,""), fChain(0), fESD(0),
+ fOutputContainer(new TObjArray(100)),
+ fPrintInfo(0), fMinGammaPt(0.),
+ fCalorimeter(""), fEMCALPID(0),fPHOSPID(0),
+ fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0),
+ fNCones(0),fNPtThres(0),fMakeICMethod(0)
+{
+ //Ctor
+ TList * list = gDirectory->GetListOfKeys() ;
+ TIter next(list) ;
+ TH2F * h = 0 ;
+ Int_t index ;
+ for (index = 0 ; index < list->GetSize()-1 ; index++) {
+ //-1 to avoid GammaJet Task
+ h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
+ fOutputContainer->Add(h) ;
+ }
+
+ for(Int_t i = 0; i < 10 ; i++){
+ fConeSizes[i]=0;
+ fPtThresholds[i]=0;
+ }
+
+ // Input slot #0 works with an Ntuple
+ DefineInput(0, TChain::Class());
+ // Output slot #0 writes into a TH1 container
+ DefineOutput(0, TObjArray::Class()) ;
+
+}
+
+
+//____________________________________________________________________________
+AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
+ AliAnalysisTask(g), fChain(g.fChain), fESD(g.fESD),
+ fOutputContainer(g. fOutputContainer), fPrintInfo(g.fPrintInfo),
+ fMinGammaPt(g.fMinGammaPt), fCalorimeter(g.fCalorimeter),
+ fEMCALPID(g.fEMCALPID),fPHOSPID(g.fPHOSPID),
+ fConeSize(g.fConeSize),
+ fPtThreshold(g.fPtThreshold),
+ fPtSumThreshold(g.fPtSumThreshold),
+ fNCones(g.fNCones),fNPtThres(g.fNPtThres),
+ fMakeICMethod(g.fMakeICMethod)
+{
+ // cpy ctor
+ SetName (g.GetName()) ;
+ SetTitle(g.GetTitle()) ;
+
+ for(Int_t i = 0; i < 10 ; i++){
+ fConeSizes[i]= g.fConeSizes[i];
+ fPtThresholds[i]= g.fPtThresholds[i];
+ }
+}
+
+//____________________________________________________________________________
+AliAnaGammaDirect::~AliAnaGammaDirect()
+{
+ // Remove all pointers
+ fOutputContainer->Clear() ;
+ delete fOutputContainer ;
+ delete fhNGamma ;
+ delete fhPhiGamma ;
+ delete fhEtaGamma ;
+ delete fhPtCandidate ;
+ delete [] fhPtThresIsolated ;
+ delete [] fhPtSumIsolated ;
+
+}
+
+//____________________________________________________________________________
+void AliAnaGammaDirect::CreateParticleList(TClonesArray * pl,
+ TClonesArray * plCTS,
+ TClonesArray * plEMCAL,
+ TClonesArray * plPHOS){
+
+ //Create a list of particles from the ESD. These particles have been measured
+ //by the Central Tracking system (TPC+ITS), PHOS and EMCAL
+
+ Int_t index = pl->GetEntries() ;
+ Int_t npar = 0 ;
+ Float_t *pid = new Float_t[AliPID::kSPECIESN];
+ AliDebug(3,"Fill particle lists");
+ if(fPrintInfo)
+ AliInfo(Form("fCalorimeter %s",fCalorimeter.Data()));
+
+ Double_t v[3] ; //vertex ;
+ fESD->GetVertex()->GetXYZ(v) ;
+
+ //########### PHOS ##############
+ if(fCalorimeter == "PHOS"){
+
+ Int_t begphos = fESD->GetFirstPHOSCluster();
+ Int_t endphos = fESD->GetFirstPHOSCluster() +
+ fESD->GetNumberOfPHOSClusters() ;
+ Int_t indexNePHOS = plPHOS->GetEntries() ;
+ AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
+
+ if(fCalorimeter == "PHOS"){
+ for (npar = begphos; npar < endphos; npar++) {//////////////PHOS track loop
+ AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
+
+ //Create a TParticle to fill the particle list
+
+ Float_t en = clus->GetClusterEnergy() ;
+ Float_t *p = new Float_t();
+ clus->GetGlobalPosition(p) ;
+ TVector3 pos(p[0],p[1],p[2]) ;
+ Double_t phi = pos.Phi();
+ Double_t theta= pos.Theta();
+ Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
+ Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
+ Double_t pz = en*TMath::Cos(theta);
+
+ TParticle * particle = new TParticle() ;
+ particle->SetMomentum(px,py,pz,en) ;
+ AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
+
+ //Select only photons
+
+ pid=clus->GetPid();
+ //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
+ if( !fPHOSPID)
+ new((*plPHOS)[indexNePHOS++]) TParticle(*particle) ;
+ else if( pid[AliPID::kPhoton] > 0.75)
+ new((*plPHOS)[indexNePHOS++]) TParticle(*particle) ;
+ }
+ }
+ }
+
+ //########### CTS (TPC+ITS) #####################
+ Int_t begtpc = 0 ;
+ Int_t endtpc = fESD->GetNumberOfTracks() ;
+ Int_t indexCh = plCTS->GetEntries() ;
+ AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
+
+ for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop
+ AliESDtrack * track = fESD->GetTrack(npar) ; // retrieve track from esd
+ //We want tracks fitted in the detectors:
+ ULong_t status=AliESDtrack::kTPCrefit;
+ status|=AliESDtrack::kITSrefit;
+
+ //We want tracks whose PID bit is set:
+ // ULong_t status =AliESDtrack::kITSpid;
+ // status|=AliESDtrack::kTPCpid;
+
+ if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
+ // Do something with the tracks which were successfully
+ // re-fitted
+ Double_t en = 0; //track ->GetTPCsignal() ;
+ Double_t mom[3];
+ track->GetPxPyPz(mom) ;
+ Double_t px = mom[0];
+ Double_t py = mom[1];
+ Double_t pz = mom[2]; //Check with TPC people if this is correct.
+ Int_t pdg = 11; //Give any charged PDG code, in this case electron.
+ //I just want to tag the particle as charged
+ TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
+ px, py, pz, en, v[0], v[1], v[2], 0);
+
+ //TParticle * particle = new TParticle() ;
+ //particle->SetMomentum(px,py,pz,en) ;
+
+ new((*plCTS)[indexCh++]) TParticle(*particle) ;
+ new((*pl)[index++]) TParticle(*particle) ;
+ }
+ }
+
+ //################ EMCAL ##############
+
+ Int_t begem = fESD->GetFirstEMCALCluster();
+ Int_t endem = fESD->GetFirstEMCALCluster() +
+ fESD->GetNumberOfEMCALClusters() ;
+ Int_t indexNe = plEMCAL->GetEntries() ;
+
+ AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
+
+ for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop
+ AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
+ Int_t clustertype= clus->GetClusterType();
+ if(clustertype == AliESDCaloCluster::kClusterv1){
+ Float_t en = clus->GetClusterEnergy() ;
+ Float_t *p = new Float_t();
+ clus->GetGlobalPosition(p) ;
+ TVector3 pos(p[0],p[1],p[2]) ;
+ Double_t phi = pos.Phi();
+ Double_t theta= pos.Theta();
+ Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
+ Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
+ Double_t pz = en*TMath::Cos(theta);
+
+ pid=clus->GetPid();
+ if(fCalorimeter == "EMCAL")
+ {
+ TParticle * particle = new TParticle() ;
+ particle->SetMomentum(px,py,pz,en) ;
+ AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
+ if(!fEMCALPID) //Only identified particles
+ new((*plEMCAL)[indexNe++]) TParticle(*particle) ;
+ else if(pid[AliPID::kPhoton] > 0.75)
+ new((*plEMCAL)[indexNe++]) TParticle(*particle) ;
+ }
+ else
+ {
+ Int_t pdg = 0;
+ if(fEMCALPID)
+ {
+ if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
+ pdg = 22;
+ else if( pid[AliPID::kPi0] > 0.75)
+ pdg = 111;
+ }
+ else
+ pdg = 22; //No PID, assume all photons
+
+ TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
+ px, py, pz, en, v[0], v[1], v[2], 0);
+ AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
+
+ new((*plEMCAL)[indexNe++]) TParticle(*particle) ;
+ new((*pl)[index++]) TParticle(*particle) ;
+ }
+ }
+ }
+
+ AliDebug(3,"Particle lists filled");
+
+}
+
+
+
+//____________________________________________________________________________
+void AliAnaGammaDirect::Exec(Option_t *)
+{
+
+ // Processing of one event
+
+ //Get ESDs
+ Long64_t entry = fChain->GetReadEntry() ;
+
+ if (!fESD) {
+ AliError("fESD is not connected to the input!") ;
+ return ;
+ }
+
+ if (fPrintInfo)
+ AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
+
+ //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
+
+ TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
+ TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
+ TClonesArray * plNe = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter (EMCAL)
+ TClonesArray * plPHOS = new TClonesArray("TParticle",1000); // All particles measured in PHOS as Gamma calorimeter
+ TClonesArray * plEMCAL = new TClonesArray("TParticle",1000); // All particles measured in EMCAL as Gamma calorimeter
+
+ TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
+ TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
+
+ //Fill lists with photons, neutral particles and charged particles
+ //look for the highest energy photon in the event inside fCalorimeter
+
+ AliDebug(2, "Fill particle lists, get prompt gamma");
+
+ //Fill particle lists
+ CreateParticleList(particleList, plCTS,plEMCAL,plPHOS);
+
+ if(fCalorimeter == "PHOS")
+ plNe = plPHOS;
+ if(fCalorimeter == "EMCAL")
+ plNe = plEMCAL;
+
+
+ //_______________Analysis 1__________________________
+ //Look for the prompt photon in the selected calorimeter
+ if(fMakeICMethod < 3){
+
+ Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
+ //Search highest energy prompt gamma in calorimeter
+ GetPromptGamma(plNe, plCTS, pGamma, iIsInPHOSorEMCAL) ;
+
+ AliDebug(1, Form("Is Gamma in %s? %d",fCalorimeter.Data(),iIsInPHOSorEMCAL));
+
+ //If there is any photon candidate in fCalorimeter
+ if(iIsInPHOSorEMCAL){
+ if (fPrintInfo)
+ AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+
+ }//Gamma in Calo
+
+ }//Analysis 1
+
+
+ //_______________Analysis 2__________________________
+ //Look for the prompt photon in the selected calorimeter
+ //Isolation Cut Analysis for both methods and different pt cuts and cones
+ if(fMakeICMethod == 3){
+
+ for(Int_t ipr = 0; ipr < plNe->GetEntries() ; ipr ++ ){
+ TParticle * pCandidate = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
+
+ if(pCandidate->Pt() > fMinGammaPt){
+
+ Bool_t icPtThres = kFALSE;
+ Bool_t icPtSum = kFALSE;
+
+ Float_t ptC = pCandidate->Pt() ;
+
+ fhPtCandidate->Fill(ptC);
+
+ for(Int_t icone = 0; icone<fNCones; icone++){
+ fConeSize = fConeSizes[icone] ;
+ Float_t coneptsum = 0 ;
+ for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
+ fPtThreshold = fPtThresholds[ipt] ;
+ MakeIsolationCut(plCTS,plNe, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
+ AliDebug(4,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
+ pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
+
+ fhPtThresIsolated[icone][ipt]->Fill(ptC);
+ }//pt thresh loop
+ fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
+ }//cone size loop
+ }//min pt candidate
+ }//candidate loop
+ }//Analysis 2
+
+ AliDebug(2, "End of analysis, delete pointers");
+
+ particleList->Delete() ;
+ plCTS->Delete() ;
+ plNe->Delete() ;
+ plPHOS->Delete() ;
+ pLeading->Delete();
+ pGamma->Delete();
+
+ delete plNe ;
+ delete plCTS ;
+ delete particleList ;
+ // delete pLeading;
+ // delete pGamma;
+
+ PostData(0, fOutputContainer);
+}
+
+
+//____________________________________________________________________________
+void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &Is) const
+{
+ //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
+
+ Double_t pt = 0;
+ Int_t index = -1;
+ for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
+ TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
+
+ if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
+ index = ipr ;
+ pt = particle->Pt();
+ pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
+ AliDebug(4,Form("Cluster in calo: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+ Is = kTRUE;
+ }
+ }
+
+ //Do Isolation?
+ if(fMakeICMethod && Is)
+ {
+ Float_t coneptsum = 0 ;
+ Bool_t icPtThres = kFALSE;
+ Bool_t icPtSum = kFALSE;
+ MakeIsolationCut(plCTS,pl, pGamma, index,
+ icPtThres, icPtSum,coneptsum);
+ if(fMakeICMethod == 1) //Pt thres method
+ Is = icPtThres ;
+ if(fMakeICMethod == 2) //Pt cone sum method
+ Is = icPtSum ;
+ }
+
+ if(Is){
+ AliDebug(3,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
+ AliDebug(3,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+ //Fill prompt gamma histograms
+ fhNGamma ->Fill(pGamma->Pt());
+ fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
+ fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
+ }
+ else
+ AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
+}
+
+ //____________________________________________________________________________
+void AliAnaGammaDirect::Init(const Option_t * )
+{
+ // Initialisation of branch container
+
+ AliDebug(2,Form("*** Initialization of %s", GetName())) ;
+
+ // Get input data
+ fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
+ if (!fChain) {
+ AliError(Form("Input 0 for %s not found\n", GetName()));
+ return ;
+ }
+
+ if (!fESD) {
+ // One should first check if the branch address was taken by some other task
+ char ** address = (char **)GetBranchAddress(0, "ESD") ;
+ if (address)
+ fESD = (AliESD *)(*address) ;
+ if (!fESD)
+ fChain->SetBranchAddress("ESD", &fESD) ;
+ }
+ // The output objects will be written to
+ TDirectory * cdir = gDirectory ;
+ // Open a file for output #0
+ char outputName[1024] ;
+ sprintf(outputName, "%s.root", GetName() ) ;
+ OpenFile(0, outputName , "RECREATE") ;
+ if (cdir)
+ cdir->cd() ;
+
+ // //Initialize the parameters of the analysis.
+ fCalorimeter="PHOS";
+ fPrintInfo = kTRUE;
+ fMinGammaPt = 10. ;
+
+ //Fill particle lists when PID is ok
+ fEMCALPID = kFALSE;
+ fPHOSPID = kFALSE;
+
+ fConeSize = 0.2 ;
+ fPtThreshold = 2.0;
+ fPtSumThreshold = 1.;
+
+ fNCones = 4 ;
+ fNPtThres = 4 ;
+ fConeSizes[0] = 0.1; fConeSizes[0] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
+ fPtThresholds[0]=1.; fPtThresholds[0]=2.; fPtThresholds[0]=3.; fPtThresholds[0]=4.;
+
+ fMakeICMethod = 1; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method, 3 make isolation study
+
+ //Initialization of histograms
+ MakeHistos() ;
+}
+
+//__________________________________________________________________
+void AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
+ TClonesArray * plNe,
+ TParticle * pCandidate,
+ Int_t index,
+ Bool_t &icmpt, Bool_t &icms,
+ Float_t &coneptsum) const
+{
+ //Search in cone around a candidate particle if it is isolated
+ Float_t phiC = pCandidate->Phi() ;
+ Float_t etaC = pCandidate->Eta() ;
+ Float_t pt = -100. ;
+ Float_t eta = -100. ;
+ Float_t phi = -100. ;
+ Float_t rad = -100 ;
+ Int_t n = 0 ;
+ TParticle * particle = new TParticle();
+
+ coneptsum = 0;
+ icmpt = kFALSE;
+ icms = kFALSE;
+
+ //Check charged particles in cone.
+ for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
+ particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
+ pt = particle->Pt();
+ eta = particle->Eta();
+ phi = particle->Phi() ;
+
+ //Check if there is any particle inside cone with pt larger than fPtThreshold
+ rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
+ TMath::Power((phi-phiC),2));
+ if(rad<fConeSize){
+ AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
+ coneptsum+=pt;
+ if(pt > fPtThreshold ) n++;
+ }
+ }// charged particle loop
+
+ //Check neutral particles in cone.
+ for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
+ if(ipr != index){//Do not count the candidate
+ particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
+ pt = particle->Pt();
+ eta = particle->Eta();
+ phi = particle->Phi() ;
+
+ //Check if there is any particle inside cone with pt larger than fPtThreshold
+ rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
+ TMath::Power((phi-phiC),2));
+ if(rad<fConeSize){
+ AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
+ coneptsum+=pt;
+ if(pt > fPtThreshold ) n++;
+ }
+ }
+ }// neutral particle loop
+
+ if(n == 0)
+ icmpt = kTRUE ;
+ if(coneptsum < fPtSumThreshold)
+ icms = kTRUE ;
+
+}
+
+//___________________________________________________________________
+void AliAnaGammaDirect::MakeHistos()
+{
+ // Create histograms to be saved in output file and
+ // stores them in fOutputContainer
+
+ fOutputContainer = new TObjArray(10000) ;
+
+ //Histograms of highest gamma identified in Event
+ fhNGamma = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120);
+ fhNGamma->SetYTitle("N");
+ fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
+ fOutputContainer->Add(fhNGamma) ;
+
+ fhPhiGamma = new TH2F
+ ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
+ fhPhiGamma->SetYTitle("#phi");
+ fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhPhiGamma) ;
+
+ fhEtaGamma = new TH2F
+ ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
+ fhEtaGamma->SetYTitle("#eta");
+ fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhEtaGamma) ;
+
+ if( fMakeICMethod== 3 ){
+
+ //Isolation cut histograms
+ fhPtCandidate = new TH1F
+ ("PtCandidate","p_{T} of candidate particles for isolation",240,0,120);
+ fhPtCandidate->SetXTitle("p_{T} (GeV/c)");
+ fOutputContainer->Add(fhPtCandidate) ;
+
+ char name[128];
+ char title[128];
+ for(Int_t icone = 0; icone<fNCones; icone++){
+ sprintf(name,"PtSumIsolated_Cone_%d",icone);
+ sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+ fhPtSumIsolated[icone] = new TH2F(name, title,240,0,120,120,0,10);
+ fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+ fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
+ fOutputContainer->Add(fhPtSumIsolated[icone]) ;
+
+ for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
+ sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
+ sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+ fhPtThresIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
+ fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+ fOutputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
+ }//icone loop
+ }//ipt loop
+ }
+}
+
+void AliAnaGammaDirect::Print(const Option_t * opt) const
+{
+
+ //Print some relevant parameters set for the analysis
+ if(! opt)
+ return;
+
+ Info("Print", "%s %s", GetName(), GetTitle() ) ;
+ printf("Cone Size = %f\n", fConeSize) ;
+ printf("pT threshold = %f\n", fPtThreshold) ;
+ printf("pT sum threshold = %f\n", fPtSumThreshold) ;
+ printf("Min Gamma pT = %f\n", fMinGammaPt) ;
+ printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
+}
+
+void AliAnaGammaDirect::Terminate(Option_t *)
+{
+ // The Terminate() function is the last function to be called during
+ // a query. It always runs on the client, it can be used to present
+ // the results graphically or save the results to file.
+
+
+}
--- /dev/null
+#ifndef ALIANAGAMMADIRECT_H
+#define ALIANAGAMMADIRECT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+
+// Class for the analysis of gamma (gamma-jet,
+// gamma-hadron(Arleo, TODO))
+// This class only contains 3 methods: one to fill lists of particles (ESDs) comming
+// from the CTS (ITS+TPC) and the calorimeters: The other search in the
+// corresponing calorimeter for the highest energy cluster, identify it as
+// prompt photon(Shower Shape (done) and Isolation Cut (TODO in new class?)).
+//
+// Class created from old AliPHOSGammaJet
+// (see AliRoot versions previous Release 4-09)
+
+//*-- Author: Gustavo Conesa (INFN-LNF)
+
+// --- ROOT system ---
+#include <TROOT.h>
+#include <TChain.h>
+#include "TTask.h"
+#include "TArrayD.h"
+#include "TChain.h"
+#include <TH2F.h>
+#include <TTree.h>
+#include <TParticle.h>
+#include "AliAnalysisTask.h"
+
+class AliESD ;
+
+class AliAnaGammaDirect : public AliAnalysisTask {
+
+public:
+
+ AliAnaGammaDirect(const char *name) ; // default ctor
+ AliAnaGammaDirect(const AliAnaGammaDirect & g) ; // cpy ctor
+ virtual ~AliAnaGammaDirect() ; //virtual dtor
+ virtual void Exec(Option_t * opt = "") ;
+ virtual void Init(Option_t * opt = "");
+ virtual void Terminate(Option_t * opt = "");
+
+ TTree * GetChain() const {return fChain ; }
+ AliESD * GetESD() const {return fESD ; }
+ TObjArray * GetOutputContainer() const {return fOutputContainer ; }
+ Double_t GetMinGammaPt() const {return fMinGammaPt ; }
+ TString GetCalorimeter() const {return fCalorimeter ; }
+ Bool_t GetPrintInfo() const {return fPrintInfo ; }
+ Float_t GetConeSize() const {return fConeSize ; }
+ Float_t GetPtThreshold() const {return fPtThreshold ; }
+ Float_t GetPtSumThres() const {return fPtSumThreshold ; }
+ Int_t GetICMethod() const {return fMakeICMethod ; }
+
+ Bool_t IsEMCALPIDOn() const {return fEMCALPID ; }
+ Bool_t IsPHOSPIDOn() const {return fPHOSPID ; }
+
+ void Print(const Option_t * opt)const;
+
+ void SetMinGammaPt(Double_t ptcut){fMinGammaPt =ptcut;}
+ void SetCalorimeter(TString calo){ fCalorimeter= calo ; }
+ void SetPrintInfo(Bool_t print){ fPrintInfo = print ; }
+ void SetConeSize(Float_t r) {fConeSize = r ; }
+ void SetPtThreshold(Float_t pt) {fPtThreshold = pt; };
+ void SetPtSumThreshold(Float_t pt) {fPtSumThreshold = pt; };
+ void SetICMethod(Int_t i ) {fMakeICMethod = i ; }
+
+ void SetConeSizes(Int_t i, Float_t r) {fConeSizes[i] = r ; }
+ void SetPtThresholds(Int_t i, Float_t pt) {fPtThresholds[i] = pt; };
+
+ void SetEMCALPIDOn(Bool_t pid){ fEMCALPID= pid ; }
+ void SetPHOSPIDOn(Bool_t pid){ fPHOSPID= pid ; }
+
+ void CreateParticleList(TClonesArray * particleList,
+ TClonesArray * plCh, TClonesArray * plNe,
+ TClonesArray * plNePHOS);
+
+
+ void GetPromptGamma(TClonesArray * plNe, TClonesArray * plCTS, TParticle * pGamma, Bool_t &Is) const;
+
+ void MakeIsolationCut(TClonesArray * plCTS, TClonesArray * plNe,
+ TParticle *pCandidate, Int_t index,
+ Bool_t &imcpt, Bool_t &icms, Float_t &ptsum) const ;
+
+ void MakeHistos() ;
+
+
+ private:
+
+ TTree *fChain ; //!pointer to the analyzed TTree or TChain
+ AliESD *fESD ; //! Declaration of leave types
+ TObjArray *fOutputContainer ; //! output data container
+ Bool_t fPrintInfo ; //Print most interesting information on screen
+ Double_t fMinGammaPt ; // Min pt in Calorimeter
+ TString fCalorimeter ; //PHOS or EMCAL detects Gamma
+ Bool_t fEMCALPID ;//Fill EMCAL particle lists with particles with corresponding pid
+ Bool_t fPHOSPID; //Fill PHOS particle lists with particles with corresponding pid
+ Float_t fConeSize ; //Size of the isolation cone
+ Float_t fPtThreshold ; //Mimium pt of the particles in the cone to set isolation
+ Float_t fPtSumThreshold ; //Mimium pt sum of the particles in the cone to set isolation
+ Int_t fNCones ; //Number of cone sizes to test
+ Int_t fNPtThres ; //Number of ptThres to test
+ Float_t fConeSizes[10] ; // Arrat with cones to test
+ Float_t fPtThresholds[10] ; // Array with pt thresholds to test
+ Int_t fMakeICMethod ; //Isolation cut method to be used
+ // 0: No isolation
+ // 1: Pt threshold method
+ // 2: Cone pt sum method
+ // 3: Study both methods for several cones and pt.
+ //Histograms
+ TH1F * fhNGamma ;
+ TH2F * fhPhiGamma ;
+ TH2F * fhEtaGamma ;
+ TH1F * fhPtCandidate ;
+ TH1F* fhPtThresIsolated[10][10] ;
+ TH2F* fhPtSumIsolated[10] ;
+
+ ClassDef(AliAnaGammaDirect,0)
+} ;
+
+
+#endif //ALIANAGAMMADIRECT_H
+
+
+
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+// Class for the analysis of gamma-jet correlations
+// Basically it seaches for a prompt photon in the Calorimeters acceptance,
+// if so we construct a jet around the highest pt particle in the opposite
+// side in azimuth, inside the Central Tracking System (ITS+TPC) and
+// EMCAL acceptances (only when PHOS detects the gamma). First the leading
+// particle and then the jet have to fullfill several conditions
+// (energy, direction ..) to be accepted. Then the fragmentation function
+// of this jet is constructed
+// Class created from old AliPHOSGammaPion
+// (see AliRoot versions previous Release 4-09)
+//
+//*-- Author: Gustavo Conesa (LNF-INFN)
+//////////////////////////////////////////////////////////////////////////////
+
+
+// --- ROOT system ---
+
+#include <TFile.h>
+#include <TParticle.h>
+#include <TH2.h>
+
+#include "AliAnaGammaHadron.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliESDCaloCluster.h"
+#include "Riostream.h"
+#include "AliLog.h"
+
+ClassImp(AliAnaGammaHadron)
+
+//____________________________________________________________________________
+AliAnaGammaHadron::AliAnaGammaHadron(const char *name) :
+ AliAnaGammaDirect(name),
+ fPhiMaxCut(0.), fPhiMinCut(0.),
+ fInvMassMaxCut(0.), fInvMassMinCut(0.),
+ fMinPtPion(0),
+ fAngleMaxParam()
+{
+
+ // ctor
+ fAngleMaxParam.Set(4) ;
+ fAngleMaxParam.Reset(0.);
+
+ TList * list = gDirectory->GetListOfKeys() ;
+ TIter next(list) ;
+ TH2F * h = 0 ;
+ Int_t index ;
+ for (index = 0 ; index < list->GetSize()-1 ; index++) {
+ //-1 to avoid GammaPion Task
+ h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
+ fOutputContainer->Add(h) ;
+ }
+
+ // Input slot #0 works with an Ntuple
+ DefineInput(0, TChain::Class());
+ // Output slot #0 writes into a TH1 container
+ DefineOutput(0, TObjArray::Class()) ;
+
+}
+
+
+//____________________________________________________________________________
+AliAnaGammaHadron::AliAnaGammaHadron(const AliAnaGammaHadron & gj) :
+ AliAnaGammaDirect(gj),
+ fPhiMaxCut(gj.fPhiMaxCut), fPhiMinCut(gj.fPhiMinCut),
+ fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut),
+ fMinPtPion(gj.fMinPtPion),
+ fOutputContainer(0), fAngleMaxParam(gj.fAngleMaxParam)
+
+{
+ // cpy ctor
+ SetName (gj.GetName()) ;
+ SetTitle(gj.GetTitle()) ;
+
+}
+
+//____________________________________________________________________________
+AliAnaGammaHadron::~AliAnaGammaHadron()
+{
+ // Remove all pointers
+ fOutputContainer->Clear() ;
+ delete fOutputContainer ;
+
+ delete fhPhiCharged ;
+ delete fhPhiNeutral ;
+ delete fhEtaCharged ;
+ delete fhEtaNeutral ;
+ delete fhDeltaPhiGammaCharged ;
+ delete fhDeltaPhiGammaNeutral ;
+ delete fhDeltaEtaGammaCharged ;
+ delete fhDeltaEtaGammaNeutral ;
+
+ delete fhCorrelationGammaNeutral ;
+ delete fhCorrelationGammaCharged ;
+
+ delete fhAnglePairNoCut ;
+ delete fhAnglePairAzimuthCut ;
+ delete fhAnglePairOpeningAngleCut ;
+ delete fhAnglePairAllCut ;
+ delete fhInvMassPairNoCut ;
+ delete fhInvMassPairAzimuthCut ;
+ delete fhInvMassPairOpeningAngleCut ;
+ delete fhInvMassPairAllCut ;
+
+}
+
+
+
+//____________________________________________________________________________
+void AliAnaGammaHadron::Exec(Option_t *)
+{
+
+ // Processing of one event
+
+ //Get ESDs
+ Long64_t entry = GetChain()->GetReadEntry() ;
+
+ if (!GetESD()) {
+ AliError("fESD is not connected to the input!") ;
+ return ;
+ }
+
+ if (GetPrintInfo())
+ AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ;
+
+ //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
+
+ TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
+ TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
+ TClonesArray * plNe = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter
+ TClonesArray * plCalo = new TClonesArray("TParticle",1000); // All particles measured in Prompt Gamma calorimeter
+
+
+ TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
+
+
+ Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
+
+ //Fill lists with photons, neutral particles and charged particles
+ //look for the highest energy photon in the event inside fCalorimeter
+ //Fill particle lists
+ AliDebug(2, "Fill particle lists, get prompt gamma");
+
+ //Fill particle lists
+ if(GetCalorimeter() == "PHOS")
+ CreateParticleList(particleList, plCTS,plNe,plCalo);
+ else if(GetCalorimeter() == "EMCAL")
+ CreateParticleList(particleList, plCTS,plCalo,plNe);
+ else
+ AliError("No calorimeter selected");
+
+ //Search highest energy prompt gamma in calorimeter
+ GetPromptGamma(plCalo, plCTS, pGamma, iIsInPHOSorEMCAL) ;
+
+
+ AliDebug(1, Form("Is Gamma in %s? %d",GetCalorimeter().Data(),iIsInPHOSorEMCAL));
+ AliDebug(3,Form("Charged list entries %d, Neutral list entries %d, %s list entries %d",
+ plCTS->GetEntries(),plNe->GetEntries(), GetCalorimeter().Data(),plCalo->GetEntries()));
+
+ //If there is any prompt photon in fCalorimeter,
+ //search jet leading particle
+
+ if(iIsInPHOSorEMCAL){
+
+ if (GetPrintInfo())
+ AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+
+ AliDebug(2, "Make correlation");
+
+ //Search correlation
+ MakeGammaChargedCorrelation(plCTS, pGamma);
+ MakeGammaNeutralCorrelation(plNe, pGamma);
+
+ }//Gamma in Calo
+
+ AliDebug(2, "End of analysis, delete pointers");
+
+ particleList->Delete() ;
+ plCTS->Delete() ;
+ plNe->Delete() ;
+ plCalo->Delete() ;
+ pGamma->Delete();
+
+ delete plNe ;
+ delete plCalo ;
+ delete plCTS ;
+ delete particleList ;
+ // delete pGamma;
+
+ PostData(0, fOutputContainer);
+}
+
+
+//____________________________________________________________________________
+void AliAnaGammaHadron::MakeGammaChargedCorrelation(TClonesArray * pl, TParticle * pGamma) const
+{
+ //Search for the charged particle with highest with
+ //Phi=Phi_gamma-Pi and pT=0.1E_gamma
+ Double_t ptg = pGamma->Pt();
+ Double_t phig = pGamma->Phi();
+ Double_t pt = -100.;
+ Double_t rat = -100.;
+ Double_t phi = -100. ;
+
+ for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
+
+ TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
+
+ pt = particle->Pt();
+ rat = pt/ptg ;
+ phi = particle->Phi() ;
+
+ AliDebug(3,Form("pt %f, phi %f, phi gamma %f. Cuts: delta phi min %f, max%f, pT min %f",pt,phi,phig,fPhiMinCut,fPhiMaxCut,fMinPtPion));
+
+ fhEtaCharged->Fill(ptg,particle->Eta());
+ fhPhiCharged->Fill(ptg,phi);
+ fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-particle->Eta());
+ fhDeltaPhiGammaCharged->Fill(ptg,phig-phi);
+ //Selection within angular and energy limits
+ if(((phig-phi)> fPhiMinCut) && ((phig-phi)<fPhiMaxCut) && pt > fMinPtPion){
+ AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
+ fhCorrelationGammaCharged->Fill(ptg,rat);
+ }
+ }//particle loop
+}
+
+//____________________________________________________________________________
+void AliAnaGammaHadron::MakeGammaNeutralCorrelation(TClonesArray * pl, TParticle * pGamma)
+{
+
+ //Search for the neutral pion with highest with
+ //Phi=Phi_gamma-Pi and pT=0.1E_gamma
+ Double_t pt = -100.;
+ Double_t rat = -100.;
+ Double_t phi = -100. ;
+ Double_t ptg = pGamma->Pt();
+ Double_t phig = pGamma->Phi();
+
+ TIter next(pl);
+ TParticle * particle = 0;
+
+ Int_t iPrimary = -1;
+ TLorentzVector gammai,gammaj;
+ Double_t angle = 0., e = 0., invmass = 0.;
+ Int_t ksPdg = 0;
+ Int_t jPrimary=-1;
+
+ while ( (particle = (TParticle*)next()) ) {
+ iPrimary++;
+ ksPdg = particle->GetPdgCode();
+
+ //2 gamma overlapped, found with PID
+ if(ksPdg == 111){
+ pt = particle->Pt();
+ rat = pt/ptg ;
+ phi = particle->Phi() ;
+ fhEtaCharged->Fill(ptg,particle->Eta());
+ fhPhiCharged->Fill(ptg,phi);
+ fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-particle->Eta());
+ fhDeltaPhiGammaCharged->Fill(ptg,phig-phi);
+ //AliDebug(3,Form("pt %f, phi %f",pt,phi));
+ if (GetPrintInfo())
+ AliInfo(Form("pt %f, phi %f",pt,phi));
+ //Selection within angular and energy limits
+ if((pt> ptg)&& ((phig-phi)>fPhiMinCut)&&((phig-phi)<fPhiMaxCut)){
+ fhCorrelationGammaNeutral ->Fill(ptg,rat);
+ //AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
+ if (GetPrintInfo())
+ AliInfo(Form("Selected: pt %f, phi %f",pt,phi));
+ }// cuts
+ }// pdg = 111
+
+ //Make invariant mass analysis
+ else if(ksPdg == 22){
+ //Search the photon companion in case it comes from a Pi0 decay
+ //Apply several cuts to select the good pair
+ particle->Momentum(gammai);
+ jPrimary=-1;
+ TIter next2(pl);
+ while ( (particle = (TParticle*)next2()) ) {
+ jPrimary++;
+ if(jPrimary>iPrimary){
+ ksPdg = particle->GetPdgCode();
+
+ if(ksPdg == 22){
+ particle->Momentum(gammaj);
+ //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
+ //gammai.Pt(),gammaj.Pt());
+ pt = (gammai+gammaj).Pt();
+ phi = (gammai+gammaj).Phi();
+ if(phi < 0)
+ phi+=TMath::TwoPi();
+ rat = pt/ptg ;
+ invmass = (gammai+gammaj).M();
+ angle = gammaj.Angle(gammai.Vect());
+ e = (gammai+gammaj).E();
+ fhEtaNeutral->Fill(ptg,(gammai+gammaj).Eta());
+ fhPhiNeutral->Fill(ptg,phi);
+ fhDeltaEtaGammaNeutral->Fill(ptg,pGamma->Eta()-(gammai+gammaj).Eta());
+ fhDeltaPhiGammaNeutral->Fill(ptg,phig-phi);
+ // AliDebug(3,Form("pt %f, phi %f",pt,phi));
+ if (GetPrintInfo())
+ AliInfo(Form("pt %f, phi %f",pt,phi));
+
+ //Fill histograms with no cuts applied.
+ fhAnglePairNoCut->Fill(e,angle);
+ fhInvMassPairNoCut->Fill(ptg,invmass);
+
+ //First cut on the energy and azimuth of the pair
+ if((phig-phi) > fPhiMinCut && (phig-phi) < fPhiMaxCut
+ && pt > fMinPtPion){
+ fhAnglePairAzimuthCut ->Fill(e,angle);
+ fhInvMassPairAzimuthCut->Fill(ptg,invmass);
+ AliDebug(3,Form("1st cut: pt %f, phi %f",pt,phi));
+
+ //Second cut on the aperture of the pair
+ if(IsAngleInWindow(angle,e)){
+ fhAnglePairOpeningAngleCut ->Fill(e,angle);
+ fhInvMassPairOpeningAngleCut->Fill(ptg,invmass);
+ AliDebug(3,Form("2nd cut: pt %f, phi %f",pt,phi));
+
+ //Third cut on the invariant mass of the pair
+ if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
+ fhInvMassPairAllCut ->Fill(ptg,invmass);
+ fhAnglePairAllCut ->Fill(e,angle);
+ //Fill correlation histogram
+ fhCorrelationGammaNeutral ->Fill(ptg,rat);
+ //AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
+ if (GetPrintInfo())
+ AliInfo(Form("Selected: pt %f, phi %f",pt,phi));
+ }//(invmass>0.125) && (invmass<0.145)
+ }//Opening angle cut
+ }//Azimuth and pt cut.
+ }//if pdg = 22
+ }//iprimary<jprimary
+ }//while
+ }// if pdg = 22
+ }//while
+}
+
+ //____________________________________________________________________________
+void AliAnaGammaHadron::Init(const Option_t * )
+{
+ // Initialisation of branch container
+ AliAnaGammaDirect::Init();
+
+ //Initialize the parameters of the analysis.
+ //fCalorimeter="PHOS";
+ fAngleMaxParam.Set(4) ;
+ fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
+ fAngleMaxParam.AddAt(-0.25,1) ;
+ fAngleMaxParam.AddAt(0.025,2) ;
+ fAngleMaxParam.AddAt(-2e-4,3) ;
+
+ //fPrintInfo = kTRUE;
+ fInvMassMaxCut = 0.16 ;
+ fInvMassMinCut = 0.11 ;
+ fPhiMaxCut = 4.5;
+ fPhiMinCut = 1.5 ;
+
+ fMinPtPion = 0. ;
+
+ //Fill particle lists when PID is ok
+ // fEMCALPID = kFALSE;
+ // fPHOSPID = kFALSE;
+
+ //Initialization of histograms
+
+ MakeHistos() ;
+
+}
+
+//__________________________________________________________________________-
+Bool_t AliAnaGammaHadron::IsAngleInWindow(const Float_t angle,const Float_t e) {
+ //Check if the opening angle of the candidate pairs is inside
+ //our selection windowd
+
+ Bool_t result = kFALSE;
+ Double_t mpi0 = 0.1349766;
+ Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
+ +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
+ Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
+ Double_t min = 100. ;
+ if(arg>0.)
+ min = TMath::ACos(arg);
+
+ if((angle<max)&&(angle>=min))
+ result = kTRUE;
+
+ return result;
+}
+
+//____________________________________________________________________________
+void AliAnaGammaHadron::MakeHistos()
+{
+ // Create histograms to be saved in output file and
+ // stores them in fOutputContainer
+
+ fOutputContainer = new TObjArray(10000) ;
+
+ //Use histograms in AliAnaGammaDirect
+ TObjArray * outputContainer =GetOutputContainer();
+ for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
+ fOutputContainer->Add(outputContainer->At(i)) ;
+
+ fhPhiCharged = new TH2F
+ ("PhiCharged","#phi_{#pi^{#pm}} vs p_{T #gamma}",
+ 120,0,120,120,0,7);
+ fhPhiCharged->SetYTitle("#phi_{#pi^{#pm}} (rad)");
+ fhPhiCharged->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhPhiCharged) ;
+
+ fhPhiNeutral = new TH2F
+ ("PhiNeutral","#phi_{#pi^{0}} vs p_{T #gamma}",
+ 120,0,120,120,0,7);
+ fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
+ fhPhiNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhPhiNeutral) ;
+
+ fhEtaCharged = new TH2F
+ ("EtaCharged","#eta_{#pi^{#pm}} vs p_{T #gamma}",
+ 120,0,120,120,-1,1);
+ fhEtaCharged->SetYTitle("#eta_{#pi^{#pm}} (rad)");
+ fhEtaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhEtaCharged) ;
+
+ fhEtaNeutral = new TH2F
+ ("EtaNeutral","#eta_{#pi^{0}} vs p_{T #gamma}",
+ 120,0,120,120,-1,1);
+ fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
+ fhEtaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhEtaNeutral) ;
+
+ fhDeltaPhiGammaCharged = new TH2F
+ ("DeltaPhiGammaCharged","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #gamma}",
+ 200,0,120,200,0,6.4);
+ fhDeltaPhiGammaCharged->SetYTitle("#Delta #phi");
+ fhDeltaPhiGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhDeltaPhiGammaCharged) ;
+
+ fhDeltaEtaGammaCharged = new TH2F
+ ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
+ 200,0,120,200,-2,2);
+ fhDeltaEtaGammaCharged->SetYTitle("#Delta #eta");
+ fhDeltaEtaGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhDeltaEtaGammaCharged) ;
+
+ fhDeltaPhiGammaNeutral = new TH2F
+ ("DeltaPhiGammaNeutral","#phi_{#gamma} - #phi_{#pi^{0}} vs p_{T #gamma}",
+ 200,0,120,200,0,6.4);
+ fhDeltaPhiGammaNeutral->SetYTitle("#Delta #phi");
+ fhDeltaPhiGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhDeltaPhiGammaNeutral) ;
+
+ fhDeltaEtaGammaNeutral = new TH2F
+ ("DeltaEtaGammaNeutral","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
+ 200,0,120,200,-2,2);
+ fhDeltaEtaGammaNeutral->SetYTitle("#Delta #eta");
+ fhDeltaEtaGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhDeltaEtaGammaNeutral) ;
+
+ //
+ fhAnglePairAccepted = new TH2F
+ ("AnglePairAccepted",
+ "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window",
+ 200,0,50,200,0,0.2);
+ fhAnglePairAccepted->SetYTitle("Angle (rad)");
+ fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairAccepted) ;
+
+ fhAnglePairNoCut = new TH2F
+ ("AnglePairNoCut",
+ "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2);
+ fhAnglePairNoCut->SetYTitle("Angle (rad)");
+ fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairNoCut) ;
+
+ fhAnglePairAzimuthCut = new TH2F
+ ("AnglePairAzimuthCut",
+ "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}",
+ 200,0,50,200,0,0.2);
+ fhAnglePairAzimuthCut->SetYTitle("Angle (rad)");
+ fhAnglePairAzimuthCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairAzimuthCut) ;
+
+ fhAnglePairOpeningAngleCut = new TH2F
+ ("AnglePairOpeningAngleCut",
+ "Angle between all #gamma pair (opening angle + azimuth cut) vs p_{T #pi^{0}}"
+ ,200,0,50,200,0,0.2);
+ fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
+ fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairOpeningAngleCut) ;
+
+ fhAnglePairAllCut = new TH2F
+ ("AnglePairAllCut",
+ "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs p_{T #pi^{0}}"
+ ,200,0,50,200,0,0.2);
+ fhAnglePairAllCut->SetYTitle("Angle (rad)");
+ fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairAllCut) ;
+
+
+ //
+ fhInvMassPairNoCut = new TH2F
+ ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
+ 120,0,120,360,0,0.5);
+ fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+ fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhInvMassPairNoCut) ;
+
+ fhInvMassPairAzimuthCut = new TH2F
+ ("InvMassPairAzimuthCut",
+ "Invariant Mass of #gamma pair (azimuth cuts) vs p_{T #gamma}",
+ 120,0,120,360,0,0.5);
+ fhInvMassPairAzimuthCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+ fhInvMassPairAzimuthCut->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhInvMassPairAzimuthCut) ;
+
+ fhInvMassPairOpeningAngleCut = new TH2F
+ ("InvMassPairOpeningAngleCut",
+ "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
+ 120,0,120,360,0,0.5);
+ fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+ fhInvMassPairOpeningAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhInvMassPairOpeningAngleCut) ;
+
+ fhInvMassPairAllCut = new TH2F
+ ("InvMassPairAllCut",
+ "Invariant Mass of #gamma pair (opening angle+invmass cut+azimuth) vs p_{T #gamma}",
+ 120,0,120,360,0,0.5);
+ fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+ fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhInvMassPairAllCut) ;
+
+ //
+ fhCorrelationGammaCharged =
+ new TH2F("CorrelationGammaCharged","z_{#gamma #pi} = p_{T #pi^{#pm}} / p_{T #gamma}",
+ 240,0.,120.,1000,0.,1.2);
+ fhCorrelationGammaCharged->SetYTitle("z_{#gamma #pi}");
+ fhCorrelationGammaCharged->SetXTitle("p_{T #gamma}");
+ fOutputContainer->Add(fhCorrelationGammaCharged) ;
+
+ fhCorrelationGammaNeutral =
+ new TH2F("CorrelationGammaNeutral","z_{#gamma #pi} = p_{T #pi^{0}} / p_{T #gamma}",
+ 240,0.,120.,1000,0.,1.2);
+ fhCorrelationGammaNeutral->SetYTitle("z_{#gamma #pi}");
+ fhCorrelationGammaNeutral->SetXTitle("p_{T #gamma}");
+ fOutputContainer->Add(fhCorrelationGammaNeutral) ;
+
+}
+
+//____________________________________________________________________________
+void AliAnaGammaHadron::Print(const Option_t * opt) const
+{
+
+ //Print some relevant parameters set for the analysis
+ if(! opt)
+ return;
+
+ Info("Print", "%s %s", GetName(), GetTitle() ) ;
+ printf("pT Pion > %f\n", fMinPtPion) ;
+ printf("Phi Pion < %f\n", fPhiMaxCut) ;
+ printf("Phi Pion > %f\n", fPhiMinCut) ;
+ printf("M_pair < %f\n", fInvMassMaxCut) ;
+ printf("M_pair > %f\n", fInvMassMinCut) ;
+
+}
+
+//__________________________________________
+void AliAnaGammaHadron::Terminate(Option_t *)
+{
+ // The Terminate() function is the last function to be called during
+ // a query. It always runs on the client, it can be used to present
+ // the results graphically or save the results to file.
+
+
+}
--- /dev/null
+#ifndef ALIANAGAMMAHADRON_H
+#define ALIANAGAMMAHADRON_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+// Class for the analysis of gamma-jet correlations.
+// Basically it seaches for a prompt photon in the Calorimeters acceptance,
+// if so we construct a jet around the highest pt particle in the opposite
+// side in azimuth. This jet has to fullfill several conditions to be
+// accepted. Then the fragmentation function of this jet is constructed
+// Class created from old AliPHOSGammaPion
+
+//*-- Author: Gustavo Conesa (INFN-LNF)
+
+// --- ROOT system ---
+#include <TROOT.h>
+#include <TChain.h>
+#include "TTask.h"
+#include "TArrayD.h"
+#include "TChain.h"
+#include <TH2F.h>
+#include <TTree.h>
+#include "AliAnaGammaDirect.h"
+
+class AliESD ;
+
+class AliAnaGammaHadron : public AliAnaGammaDirect {
+
+public:
+
+ AliAnaGammaHadron(const char *name) ; // default ctor
+ AliAnaGammaHadron(const AliAnaGammaHadron & gj) ; // cpy ctor
+ virtual ~AliAnaGammaHadron() ; //virtual dtor
+ virtual void Exec(Option_t * opt = "") ;
+ virtual void Init(Option_t * opt = "");
+ virtual void Terminate(Option_t * opt = "");
+
+ Double_t GetAngleMaxParam(Int_t i) const {return fAngleMaxParam.At(i) ; }
+ Double_t GetInvMassMaxCut() const {return fInvMassMaxCut ; }
+ Double_t GetInvMassMinCut() const {return fInvMassMinCut ; }
+ Double_t GetPhiMaxCut() const {return fPhiMaxCut ; }
+ Double_t GetPhiMinCut() const {return fPhiMinCut ; }
+ Float_t GetMinPtPion() const {return fMinPtPion ; }
+
+ void Print(const Option_t * opt)const;
+
+ void SetAngleMaxParam(Int_t i, Double_t par)
+ {fAngleMaxParam.AddAt(par,i) ; }
+ void SetMinPtPion(Float_t pt){fMinPtPion = pt; };
+ void SetInvMassCutRange(Double_t invmassmin, Double_t invmassmax)
+ {fInvMassMaxCut =invmassmax; fInvMassMinCut =invmassmin;}
+ void SetPhiCutRange(Double_t phimin, Double_t phimax)
+ {fPhiMaxCut =phimax; fPhiMinCut =phimin;}
+
+ private:
+
+ Bool_t IsAngleInWindow(const Float_t angle, const Float_t e);
+ void MakeGammaChargedCorrelation(TClonesArray * pl, TParticle *pGamma) const ;
+ void MakeGammaNeutralCorrelation(TClonesArray * pl, TParticle *pGamma) ;
+ void MakeHistos() ;
+
+ private:
+
+ // TTree *fChain ; //!pointer to the analyzed TTree or TChain
+ //AliESD *fESD ; //! Declaration of leave types
+
+ Double_t fPhiMaxCut ; //
+ Double_t fPhiMinCut ; //
+ // Double_t fGammaPtCut ; // Min pt in Calorimeter
+ Double_t fInvMassMaxCut ; // Invariant Mass cut maximum
+ Double_t fInvMassMinCut ; // Invariant Masscut minimun
+
+ Double_t fMinPtPion; // Minimum pt of pion
+ TObjArray *fOutputContainer ; //! output data container
+ TString fNamePtThres[10]; // String name of pt th to append to histos
+ TArrayD fAngleMaxParam ; //Max opening angle selection parameters
+ //TString fCalorimeter ; //PHOS or EMCAL detects Gamma
+ //Bool_t fEMCALPID ; //Fill EMCAL particle lists with particles with corresponding pid
+ //Bool_t fPHOSPID; //Fill PHOS particle lists with particles with corresponding pid
+
+ //Histograms
+
+ TH2F * fhPhiCharged ;
+ TH2F * fhPhiNeutral ;
+ TH2F * fhEtaCharged ;
+ TH2F * fhEtaNeutral ;
+ TH2F * fhDeltaPhiGammaCharged ;
+ TH2F * fhDeltaPhiGammaNeutral ;
+ TH2F * fhDeltaEtaGammaCharged ;
+ TH2F * fhDeltaEtaGammaNeutral ;
+
+ TH2F * fhCorrelationGammaNeutral ;
+ TH2F * fhCorrelationGammaCharged ;
+
+ TH2F * fhAnglePairAccepted ;
+ TH2F * fhAnglePairNoCut ;
+ TH2F * fhAnglePairAzimuthCut ;
+ TH2F * fhAnglePairOpeningAngleCut ;
+ TH2F * fhAnglePairAllCut ;
+ TH2F * fhInvMassPairNoCut ;
+ TH2F * fhInvMassPairAzimuthCut ;
+ TH2F * fhInvMassPairOpeningAngleCut ;
+ TH2F * fhInvMassPairAllCut ;
+
+ ClassDef(AliAnaGammaHadron,0)
+} ;
+
+
+#endif //ALIANAGAMMAHADRON_H
+
+
+
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+// Class for the analysis of gamma-jet correlations
+// Basically it seaches for a prompt photon in the Calorimeters acceptance,
+// if so we construct a jet around the highest pt particle in the opposite
+// side in azimuth, inside the Central Tracking System (ITS+TPC) and
+// EMCAL acceptances (only when PHOS detects the gamma). First the leading
+// particle and then the jet have to fullfill several conditions
+// (energy, direction ..) to be accepted. Then the fragmentation function
+// of this jet is constructed
+// Class created from old AliPHOSGammaJet
+// (see AliRoot versions previous Release 4-09)
+//
+//*-- Author: Gustavo Conesa (LNF-INFN)
+//////////////////////////////////////////////////////////////////////////////
+
+
+// --- ROOT system ---
+
+#include <TFile.h>
+#include <TParticle.h>
+#include <TH2.h>
+
+#include "AliAnaGammaJet.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliESDCaloCluster.h"
+#include "Riostream.h"
+#include "AliLog.h"
+
+ClassImp(AliAnaGammaJet)
+
+//____________________________________________________________________________
+AliAnaGammaJet::AliAnaGammaJet(const char *name) :
+ AliAnaGammaDirect(name),
+ fSeveralConeAndPtCuts(0),
+ fPbPb(kFALSE),
+ fJetsOnlyInCTS(0),
+ fEtaEMCALCut(0.),fPhiMaxCut(0.),
+ fPhiMinCut(0.),
+ fInvMassMaxCut(0.), fInvMassMinCut(0.),
+ fJetCTSRatioMaxCut(0.),
+ fJetCTSRatioMinCut(0.), fJetRatioMaxCut(0.),
+ fJetRatioMinCut(0.), fNCone(0),
+ fNPt(0), fCone(0), fPtThreshold(0),
+ fPtJetSelectionCut(0.0),
+ fAngleMaxParam(), fSelect(0)
+{
+
+ // ctor
+ fAngleMaxParam.Set(4) ;
+ fAngleMaxParam.Reset(0.);
+
+ for(Int_t i = 0; i<10; i++){
+ fCones[i] = 0.0 ;
+ fNameCones[i] = "" ;
+ fPtThres[i] = 0.0 ;
+ fNamePtThres[i] = "" ;
+ if( i < 6 ){
+ fJetXMin1[i] = 0.0 ;
+ fJetXMin2[i] = 0.0 ;
+ fJetXMax1[i] = 0.0 ;
+ fJetXMax2[i] = 0.0 ;
+ fBkgMean[i] = 0.0 ;
+ fBkgRMS[i] = 0.0 ;
+ if( i < 2 ){
+ fJetE1[i] = 0.0 ;
+ fJetE2[i] = 0.0 ;
+ fJetSigma1[i] = 0.0 ;
+ fJetSigma2[i] = 0.0 ;
+ fPhiEMCALCut[i] = 0.0 ;
+ }
+ }
+ }
+
+ TList * list = gDirectory->GetListOfKeys() ;
+ TIter next(list) ;
+ TH2F * h = 0 ;
+ Int_t index ;
+ for (index = 0 ; index < list->GetSize()-1 ; index++) {
+ //-1 to avoid GammaJet Task
+ h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
+ fOutputContainer->Add(h) ;
+ }
+
+ // Input slot #0 works with an Ntuple
+ DefineInput(0, TChain::Class());
+ // Output slot #0 writes into a TH1 container
+ DefineOutput(0, TObjArray::Class()) ;
+
+}
+
+
+//____________________________________________________________________________
+AliAnaGammaJet::AliAnaGammaJet(const AliAnaGammaJet & gj) :
+ AliAnaGammaDirect(gj),
+ fSeveralConeAndPtCuts(gj.fSeveralConeAndPtCuts),
+ fPbPb(gj.fPbPb), fJetsOnlyInCTS(gj.fJetsOnlyInCTS),
+ fEtaEMCALCut(gj.fEtaEMCALCut),
+ fPhiMaxCut(gj.fPhiMaxCut), fPhiMinCut(gj.fPhiMinCut),
+ fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut),
+ fRatioMinCut(gj.fRatioMinCut),
+ fJetCTSRatioMaxCut(gj.fJetCTSRatioMaxCut),
+ fJetCTSRatioMinCut(gj.fJetCTSRatioMinCut), fJetRatioMaxCut(gj.fJetRatioMaxCut),
+ fJetRatioMinCut(gj.fJetRatioMinCut), fNCone(gj.fNCone),
+ fNPt(gj.fNPt), fCone(gj.fCone), fPtThreshold(gj.fPtThreshold),
+ fPtJetSelectionCut(gj.fPtJetSelectionCut),
+ fOutputContainer(0), fAngleMaxParam(gj.fAngleMaxParam),
+ fSelect(gj.fSelect)
+{
+ // cpy ctor
+ SetName (gj.GetName()) ;
+ SetTitle(gj.GetTitle()) ;
+
+ for(Int_t i = 0; i<10; i++){
+ fCones[i] = gj.fCones[i] ;
+ fNameCones[i] = gj.fNameCones[i] ;
+ fPtThres[i] = gj.fPtThres[i] ;
+ fNamePtThres[i] = gj.fNamePtThres[i] ;
+ if( i < 6 ){
+ fJetXMin1[i] = gj.fJetXMin1[i] ;
+ fJetXMin2[i] = gj.fJetXMin2[i] ;
+ fJetXMax1[i] = gj.fJetXMax1[i] ;
+ fJetXMax2[i] = gj.fJetXMax2[i] ;
+ fBkgMean[i] = gj.fBkgMean[i] ;
+ fBkgRMS[i] = gj.fBkgRMS[i] ;
+ if( i < 2 ){
+ fJetE1[i] = gj.fJetE1[i] ;
+ fJetE2[i] = gj.fJetE2[i] ;
+ fJetSigma1[i] = gj.fJetSigma1[i] ;
+ fJetSigma2[i] = gj.fJetSigma2[i] ;
+ fPhiEMCALCut[i] = gj.fPhiEMCALCut[i] ;
+ }
+ }
+ }
+}
+
+//____________________________________________________________________________
+AliAnaGammaJet::~AliAnaGammaJet()
+{
+ // Remove all pointers
+ fOutputContainer->Clear() ;
+ delete fOutputContainer ;
+
+ delete fhChargeRatio ;
+ delete fhPi0Ratio ;
+ delete fhDeltaPhiCharge ;
+ delete fhDeltaPhiPi0 ;
+ delete fhDeltaEtaCharge ;
+ delete fhDeltaEtaPi0 ;
+ delete fhAnglePair ;
+ delete fhAnglePairAccepted ;
+ delete fhAnglePairNoCut ;
+ delete fhAnglePairLeadingCut ;
+ delete fhAnglePairAngleCut ;
+ delete fhAnglePairAllCut ;
+ delete fhAnglePairLeading ;
+ delete fhInvMassPairNoCut ;
+ delete fhInvMassPairLeadingCut ;
+ delete fhInvMassPairAngleCut ;
+ delete fhInvMassPairAllCut ;
+ delete fhInvMassPairLeading ;
+ delete fhNBkg ;
+ delete fhNLeading ;
+ delete fhNJet ;
+ delete fhJetRatio ;
+ delete fhJetPt ;
+ delete fhBkgRatio ;
+ delete fhBkgPt ;
+ delete fhJetFragment ;
+ delete fhBkgFragment ;
+ delete fhJetPtDist ;
+ delete fhBkgPtDist ;
+
+ delete [] fhJetRatios;
+ delete [] fhJetPts;
+ delete [] fhBkgRatios;
+ delete [] fhBkgPts;
+
+ delete [] fhNLeadings;
+ delete [] fhNJets;
+ delete [] fhNBkgs;
+
+ delete [] fhJetFragments;
+ delete [] fhBkgFragments;
+ delete [] fhJetPtDists;
+ delete [] fhBkgPtDists;
+
+}
+
+//____________________________________________________________________________
+Double_t AliAnaGammaJet::CalculateJetRatioLimit(const Double_t ptg,
+ const Double_t *par,
+ const Double_t *x) {
+
+ //Parametrized cut for the energy of the jet.
+
+ Double_t epp = par[0] + par[1] * ptg ;
+ Double_t spp = par[2] + par[3] * ptg ;
+ Double_t f = x[0] + x[1] * ptg ;
+ Double_t epb = epp + par[4] ;
+ Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
+ Double_t rat = (epb - spb * f) / ptg ;
+ //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
+ //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
+ return rat ;
+}
+
+
+//____________________________________________________________________________
+void AliAnaGammaJet::Exec(Option_t *)
+{
+
+ // Processing of one event
+
+ //Get ESDs
+ Long64_t entry = GetChain()->GetReadEntry() ;
+
+ if (!GetESD()) {
+ AliError("fESD is not connected to the input!") ;
+ return ;
+ }
+
+ if (GetPrintInfo())
+ AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ;
+
+ //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
+
+ TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
+ TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
+ TClonesArray * plNe = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter
+ TClonesArray * plCalo = new TClonesArray("TParticle",1000); // All particles measured in Prompt Gamma calorimeter
+
+
+ TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
+ TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
+
+ Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
+
+ //Fill lists with photons, neutral particles and charged particles
+ //look for the highest energy photon in the event inside fCalorimeter
+ AliDebug(2, "Fill particle lists, get prompt gamma");
+
+ //Fill particle lists
+
+ if(GetCalorimeter() == "PHOS")
+ CreateParticleList(particleList, plCTS,plNe,plCalo);
+ else if(GetCalorimeter() == "EMCAL")
+ CreateParticleList(particleList, plCTS,plCalo,plNe);
+ else
+ AliError("No calorimeter selected");
+
+ //Search highest energy prompt gamma in calorimeter
+ GetPromptGamma(plCalo, plCTS, pGamma, iIsInPHOSorEMCAL) ;
+
+ AliDebug(1, Form("Is Gamma in %s? %d",GetCalorimeter().Data(),iIsInPHOSorEMCAL));
+ AliDebug(3,Form("Charged list entries %d, Neutral list entries %d, %s list entries %d",
+ plCTS->GetEntries(),plNe->GetEntries(), GetCalorimeter().Data(), plCalo->GetEntries()));
+
+ //If there is any prompt photon in fCalorimeter,
+ //search jet leading particle
+
+ if(iIsInPHOSorEMCAL){
+ if (GetPrintInfo())
+ AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+
+ AliDebug(2, "Get Leading Particles, Make jet");
+
+ //Search leading particles in CTS and EMCAL
+ if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){
+ if (GetPrintInfo())
+ AliInfo(Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ;
+
+ //Search Jet
+ if(!fSeveralConeAndPtCuts)
+ MakeJet(particleList,pGamma,pLeading,"");
+ else{
+ for(Int_t icone = 0; icone<fNCone; icone++) {
+ for(Int_t ipt = 0; ipt<fNPt;ipt++) {
+ TString lastname ="Cone"+ fNameCones[icone]+"Pt"+ fNamePtThres[ipt];
+ MakeJet(particleList, pGamma, pLeading,lastname);
+ }//icone
+ }//ipt
+ }//fSeveralConeAndPtCuts
+ }//Leading
+ }//Gamma in Calo
+
+ AliDebug(2, "End of analysis, delete pointers");
+
+ particleList->Delete() ;
+ plCTS->Delete() ;
+ plNe->Delete() ;
+ plCalo->Delete() ;
+ pLeading->Delete();
+ pGamma->Delete();
+
+ delete plNe ;
+ delete plCalo ;
+ delete plCTS ;
+ delete particleList ;
+ // delete pLeading;
+ // delete pGamma;
+
+ PostData(0, fOutputContainer);
+}
+
+//____________________________________________________________________________
+void AliAnaGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname)
+{
+ //Fill histograms wth jet fragmentation
+ //and number of selected jets and leading particles
+ //and the background multiplicity
+ TParticle * particle = 0 ;
+ Int_t ipr = 0;
+ Float_t charge = 0;
+
+ TIter next(pl) ;
+ while ( (particle = dynamic_cast<TParticle*>(next())) ) {
+ ipr++ ;
+ Double_t pt = particle->Pt();
+
+ charge = TDatabasePDG::Instance()
+ ->GetParticle(particle->GetPdgCode())->Charge();
+ if(charge != 0){//Only jet Charged particles
+ dynamic_cast<TH2F*>
+ (fOutputContainer->FindObject(type+"Fragment"+lastname))
+ ->Fill(ptg,pt/ptg);
+ dynamic_cast<TH2F*>
+ (fOutputContainer->FindObject(type+"PtDist"+lastname))
+ ->Fill(ptg,pt);
+ }//charged
+
+ }//while
+
+ if(type == "Bkg")
+ dynamic_cast<TH1F*>
+ (fOutputContainer->FindObject("NBkg"+lastname))
+ ->Fill(ipr);
+ else{
+ dynamic_cast<TH1F*>
+ (fOutputContainer->FindObject("NJet"+lastname))->
+ Fill(ptg);
+ dynamic_cast<TH2F*>
+ (fOutputContainer->FindObject("NLeading"+lastname))
+ ->Fill(ptg,ptl);
+ }
+
+}
+
+//____________________________________________________________________________
+void AliAnaGammaJet::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const
+{
+ //Search for the charged particle with highest with
+ //Phi=Phi_gamma-Pi and pT=0.1E_gamma
+ Double_t pt = -100.;
+ Double_t phi = -100.;
+
+ for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
+
+ TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
+
+ Double_t ptl = particle->Pt();
+ Double_t rat = ptl/pGamma->Pt() ;
+ Double_t phil = particle->Phi() ;
+ Double_t phig = pGamma->Phi();
+
+ //Selection within angular and energy limits
+ if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
+ (rat > fRatioMinCut) && (rat < fRatioMaxCut) && (ptl > pt)) {
+ phi = phil ;
+ pt = ptl ;
+ pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
+ AliDebug(4,Form("Charge in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, particle->Eta(), phi,rat)) ;
+ }
+ }
+
+ AliDebug(3,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ;
+
+}
+
+
+//____________________________________________________________________________
+void AliAnaGammaJet::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading)
+{
+
+ //Search for the neutral pion with highest with
+ //Phi=Phi_gamma-Pi and pT=0.1E_gamma
+ Double_t pt = -100.;
+ Double_t phi = -100.;
+ Double_t ptl = -100.;
+ Double_t rat = -100.;
+ Double_t phil = -100. ;
+ Double_t ptg = pGamma->Pt();
+ Double_t phig = pGamma->Phi();
+
+ TIter next(pl);
+ TParticle * particle = 0;
+
+ Int_t iPrimary = -1;
+ TLorentzVector gammai,gammaj;
+ Double_t angle = 0., e = 0., invmass = 0.;
+ Double_t anglef = 0., ef = 0., invmassf = 0.;
+ Int_t ksPdg = 0;
+ Int_t jPrimary=-1;
+
+ while ( (particle = (TParticle*)next()) ) {
+ iPrimary++;
+
+ ksPdg = particle->GetPdgCode();
+ ptl = particle->Pt();
+ if(ksPdg == 111){ //2 gamma overlapped, found with PID
+ rat = ptl/ptg ;
+ phil = particle->Phi() ;
+ //Selection within angular and energy limits
+ if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
+ ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
+ phi = phil ;
+ pt = ptl ;
+ pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
+ AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ;
+ }// cuts
+ }// pdg = 111
+ else if(ksPdg == 22){//1 gamma
+ //Search the photon companion in case it comes from a Pi0 decay
+ //Apply several cuts to select the good pair
+ particle->Momentum(gammai);
+ jPrimary=-1;
+ TIter next2(pl);
+ while ( (particle = (TParticle*)next2()) ) {
+ jPrimary++;
+ if(jPrimary>iPrimary){
+ ksPdg = particle->GetPdgCode();
+
+ if(ksPdg == 22){
+ particle->Momentum(gammaj);
+ //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
+ //gammai.Pt(),gammaj.Pt());
+ ptl = (gammai+gammaj).Pt();
+ phil = (gammai+gammaj).Phi();
+ if(phil < 0)
+ phil+=TMath::TwoPi();
+ rat = ptl/ptg ;
+ invmass = (gammai+gammaj).M();
+ angle = gammaj.Angle(gammai.Vect());
+ //Info("GetLeadingPi0","Angle %f", angle);
+ e = (gammai+gammaj).E();
+ //Fill histograms with no cuts applied.
+ fhAnglePairNoCut->Fill(e,angle);
+ fhInvMassPairNoCut->Fill(ptg,invmass);
+ //First cut on the energy and azimuth of the pair
+ if((rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
+ ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
+
+ fhAnglePairLeadingCut->Fill(e,angle);
+ fhInvMassPairLeadingCut->Fill(ptg,invmass);
+ //Second cut on the aperture of the pair
+ if(IsAngleInWindow(angle,e)){
+ fhAnglePairAngleCut->Fill(e,angle);
+ fhInvMassPairAngleCut->Fill(ptg,invmass);
+
+ //Info("GetLeadingPi0","InvMass %f", invmass);
+ //Third cut on the invariant mass of the pair
+ if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
+ fhInvMassPairAllCut->Fill(ptg,invmass);
+ fhAnglePairAllCut->Fill(e,angle);
+ if(ptl > pt ){
+ pt = ptl;
+ phi = phil ;
+ ef = e ;
+ anglef = angle ;
+ invmassf = invmass ;
+ pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
+ AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ;
+ }
+ }//cuts
+ }//(invmass>0.125) && (invmass<0.145)
+ }//gammaj.Angle(gammai.Vect())<0.04
+ }//if pdg = 22
+ }//iprimary<jprimary
+ }//while
+ }// if pdg = 22
+ // }
+ }//while
+
+ if(ef > 0.){//Final pi0 found, highest pair energy, fill histograms
+ fhInvMassPairLeading->Fill(ptg,invmassf);
+ fhAnglePairLeading->Fill(ef,anglef);
+ }
+
+ AliDebug(3,Form("Leading EMCAL: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ;
+}
+
+//____________________________________________________________________________
+Bool_t AliAnaGammaJet::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe,
+ TParticle * pGamma, TParticle * pLeading)
+{
+ //Search Charged or Neutral leading particle, select the highest one.
+
+ TParticle * pLeadingCh = new TParticle();
+ TParticle * pLeadingPi0 = new TParticle();
+
+ Double_t ptg = pGamma->Pt();
+ Double_t phig = pGamma->Phi();
+ Double_t etag = pGamma->Eta();
+
+ if(GetCalorimeter() == "PHOS" && !fJetsOnlyInCTS)
+ {
+ AliDebug(3, "GetLeadingPi0");
+ GetLeadingPi0 (plNe, pGamma, pLeadingPi0) ;
+ AliDebug(3, "GetLeadingCharge");
+ GetLeadingCharge(plCTS, pGamma, pLeadingCh) ;
+
+ Double_t ptch = pLeadingCh->Pt();
+ Double_t phich = pLeadingCh->Phi();
+ Double_t etach = pLeadingCh->Eta();
+ Double_t ptpi = pLeadingPi0->Pt();
+ Double_t phipi = pLeadingPi0->Phi();
+ Double_t etapi = pLeadingPi0->Eta();
+
+ //Is leading cone inside EMCAL acceptance?
+
+ Bool_t insidech = kFALSE ;
+ if((phich - fCone) > fPhiEMCALCut[0] &&
+ (phich + fCone) < fPhiEMCALCut[1] &&
+ (etach-fCone) < fEtaEMCALCut )
+ insidech = kTRUE ;
+
+ Bool_t insidepi = kFALSE ;
+ if((phipi - fCone) > fPhiEMCALCut[0] &&
+ (phipi + fCone) < fPhiEMCALCut[1] &&
+ (etapi-fCone) < fEtaEMCALCut )
+ insidepi = kTRUE ;
+
+ AliDebug(2,Form("Leading: charged pt %f, pi0 pt %f",ptch,ptpi)) ;
+
+ if (ptch > 0 || ptpi > 0)
+ {
+ if(insidech && (ptch > ptpi))
+ {
+ if (GetPrintInfo())
+ AliInfo("Leading found in CTS");
+ pLeading->SetMomentum(pLeadingCh->Px(),pLeadingCh->Py(),pLeadingCh->Pz(),pLeadingCh->Energy());
+ AliDebug(3,Form("Final leading found in CTS, pt %f, phi %f, eta %f",ptch,phich,etach)) ;
+ fhChargeRatio->Fill(ptg,ptch/pGamma->Pt());
+ fhDeltaPhiCharge->Fill(ptg,pGamma->Phi()-phich);
+ fhDeltaEtaCharge->Fill(ptg,pGamma->Eta()-etach);
+ return 1;
+ }
+
+ else if((ptpi > ptch) && insidepi)
+ {
+ if (GetPrintInfo())
+ AliInfo("Leading found in EMCAL");
+ pLeading->SetMomentum(pLeadingPi0->Px(),pLeadingPi0->Py(),pLeadingPi0->Pz(),pLeadingPi0->Energy());
+ AliDebug(3,Form("Final leading found in EMCAL, pt %f, phi %f, eta %f",ptpi,phipi,etapi)) ;
+ fhPi0Ratio ->Fill(ptg,ptpi/ptg);
+ fhDeltaPhiPi0->Fill(ptg,phig-phipi);
+ fhDeltaEtaPi0->Fill(ptg,etag-etapi);
+ return 1;
+ }
+
+ else{
+ AliDebug(3,"NO LEADING PARTICLE FOUND");}
+ return 0;
+ }
+ else{
+ AliDebug(3,"NO LEADING PARTICLE FOUND");
+ return 0;
+ }
+ }
+
+ else
+ {
+ //No calorimeter present for Leading particle detection
+ GetLeadingCharge(plCTS, pGamma, pLeading) ;
+ Double_t ptch = pLeading->Pt();
+ Double_t phich = pLeading->Phi();
+ Double_t etach = pLeading->Eta();
+ if(ptch > 0){
+ fhChargeRatio->Fill(ptg,ptch/ptg);
+ fhDeltaPhiCharge->Fill(ptg,phig-phich);
+ fhDeltaEtaCharge->Fill(ptg,etag-etach);
+ AliDebug(3,Form("Leading found : pt %f, phi %f, eta %f",ptch,phich,etach)) ;
+ return 1;
+ }
+ else
+ {
+ AliDebug(3,"NO LEADING PARTICLE FOUND");
+ return 0;
+ }
+ }
+}
+
+ //____________________________________________________________________________
+void AliAnaGammaJet::Init(const Option_t * )
+{
+// // Initialisation of branch container
+ AliAnaGammaDirect::Init();
+
+ //Initialize the parameters of the analysis.
+ //fCalorimeter="PHOS";
+ fAngleMaxParam.Set(4) ;
+ fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
+ fAngleMaxParam.AddAt(-0.25,1) ;
+ fAngleMaxParam.AddAt(0.025,2) ;
+ fAngleMaxParam.AddAt(-2e-4,3) ;
+ fSeveralConeAndPtCuts = kFALSE ;
+ //fPrintInfo = kTRUE;
+ fPbPb = kFALSE ;
+ fInvMassMaxCut = 0.16 ;
+ fInvMassMinCut = 0.11 ;
+ //fJetsOnlyInCTS = kTRUE ;
+ fEtaEMCALCut = 0.7 ;
+ fPhiEMCALCut[0] = 80. *TMath::Pi()/180.;
+ fPhiEMCALCut[1] = 190.*TMath::Pi()/180.;
+ fPhiMaxCut = 3.4 ;
+ fPhiMinCut = 2.9 ;
+
+ //Jet selection parameters
+ //Fixed cut (old)
+ fRatioMaxCut = 1.0 ;
+ fRatioMinCut = 0.1 ;
+ fJetRatioMaxCut = 1.2 ;
+ fJetRatioMinCut = 0.3 ;
+ fJetsOnlyInCTS = kFALSE ;
+ fJetCTSRatioMaxCut = 1.2 ;
+ fJetCTSRatioMinCut = 0.3 ;
+ fSelect = 0 ;
+
+ //Cut depending on gamma energy
+
+ fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied
+ //Reconstructed jet energy dependence parameters
+ //e_jet = a1+e_gamma b2.
+ //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
+ fJetE1[0] = -5.75; fJetE1[1] = -4.1;
+ fJetE2[0] = 1.005; fJetE2[1] = 1.05;
+
+ //Reconstructed sigma of jet energy dependence parameters
+ //s_jet = a1+e_gamma b2.
+ //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
+ fJetSigma1[0] = 2.65; fJetSigma1[1] = 2.75;
+ fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
+
+ //Background mean energy and RMS
+ //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
+ //Index 2-> (low pt jets)BKG > 0.5 GeV;
+ //Index > 2, same for CTS conf
+ fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
+ fBkgMean[3] = 0.; fBkgMean[4] = 6.4; fBkgMean[5] = 48.6;
+ fBkgRMS[0] = 0.; fBkgRMS[1] = 7.5; fBkgRMS[2] = 22.0;
+ fBkgRMS[3] = 0.; fBkgRMS[4] = 5.4; fBkgRMS[5] = 13.2;
+
+ //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
+ //limits for monoenergetic jets.
+ //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
+ //Index 2-> (low pt jets) BKG > 0.5 GeV;
+ //Index > 2, same for CTS conf
+
+ fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
+ fJetXMin1[3] =-2.0 ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1 ;
+ fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
+ fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
+ fJetXMax1[0] =-3.8 ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6 ;
+ fJetXMax1[3] =-2.7 ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7 ;
+ fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
+ fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
+
+
+ //Different cones and pt thresholds to construct the jet
+
+ fCone = 0.3 ;
+ fPtThreshold = 0. ;
+ fNCone = 3 ;
+ fNPt = 3 ;
+ fCones[1] = 0.2 ; fNameCones[1] = "02" ;
+ fPtThres[0] = 0.0 ; fNamePtThres[0] = "00" ;
+ fCones[0] = 0.3 ; fNameCones[0] = "03" ;
+ fPtThres[1] = 0.5 ; fNamePtThres[1] = "05" ;
+ fCones[2] = 0.4 ; fNameCones[2] = "04" ;
+ fPtThres[2] = 1.0 ; fNamePtThres[2] = "10" ;
+ //Fill particle lists when PID is ok
+ // fEMCALPID = kFALSE;
+ // fPHOSPID = kFALSE;
+
+ //Initialization of histograms
+
+ MakeHistos() ;
+
+}
+
+//__________________________________________________________________________-
+Bool_t AliAnaGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
+ //Check if the opening angle of the candidate pairs is inside
+ //our selection windowd
+
+ Bool_t result = kFALSE;
+ Double_t mpi0 = 0.1349766;
+ Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
+ +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
+ Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
+ Double_t min = 100. ;
+ if(arg>0.)
+ min = TMath::ACos(arg);
+
+ if((angle<max)&&(angle>=min))
+ result = kTRUE;
+
+ return result;
+}
+
+//__________________________________________________________________________-
+Bool_t AliAnaGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj){
+ //Check if the energy of the reconstructed jet is within an energy window
+
+ Double_t par[6];
+ Double_t xmax[2];
+ Double_t xmin[2];
+
+ Int_t iCTS = 0;
+ if(fJetsOnlyInCTS)
+ iCTS = 3 ;
+
+ if(!fPbPb){
+ //Phythia alone, jets with pt_th > 0.2, r = 0.3
+ par[0] = fJetE1[0]; par[1] = fJetE2[0];
+ //Energy of the jet peak
+ //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
+ par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
+ //Sigma of the jet peak
+ //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
+ par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
+ //Parameters reserved for PbPb bkg.
+ xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
+ xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
+ //Factor that multiplies sigma to obtain the best limits,
+ //by observation, of mono jet ratios (ptjet/ptg)
+ //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
+
+ }
+ else{
+ if(ptg > fPtJetSelectionCut){
+ //Phythia +PbPb with pt_th > 2 GeV/c, r = 0.3
+ par[0] = fJetE1[0]; par[1] = fJetE2[0];
+ //Energy of the jet peak, same as in pp
+ //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
+ par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
+ //Sigma of the jet peak, same as in pp
+ //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
+ par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
+ //Mean value and RMS of PbPb Bkg
+ xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
+ xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
+ //Factor that multiplies sigma to obtain the best limits,
+ //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg,
+ //pt_th > 2 GeV, r = 0.3
+ //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
+
+ }
+ else{
+ //Phythia + PbPb with pt_th > 0.5 GeV/c, r = 0.3
+ par[0] = fJetE1[1]; par[1] = fJetE2[1];
+ //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
+ //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
+ par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
+ //Sigma of the jet peak, pt_th > 2 GeV/c, r = 0.3
+ //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
+ par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
+ //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
+ xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
+ xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
+ //Factor that multiplies sigma to obtain the best limits,
+ //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg,
+ //pt_th > 2 GeV, r = 0.3
+ //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
+
+ }//If low pt jet in bkg
+ }//if Bkg
+
+ //Calculate minimum and maximum limits of the jet ratio.
+ Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
+ Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
+
+ AliDebug(3,Form("Jet selection? : Limits min %f, max %f, pt_jet %f, pt_gamma %f, pt_jet / pt_gamma %f",min,max,ptj,ptg,ptj/ptg));
+
+ if(( min < ptj/ptg ) && ( max > ptj/ptg))
+ return kTRUE;
+ else
+ return kFALSE;
+
+}
+
+//____________________________________________________________________________
+void AliAnaGammaJet::MakeHistos()
+{
+ // Create histograms to be saved in output file and
+ // stores them in fOutputContainer
+
+ fOutputContainer = new TObjArray(10000) ;
+
+ TObjArray * outputContainer =GetOutputContainer();
+ for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
+ fOutputContainer->Add(outputContainer->At(i)) ;
+
+ //
+ fhChargeRatio = new TH2F
+ ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
+ 120,0,120,120,0,1);
+ fhChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
+ fhChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhChargeRatio) ;
+
+ fhDeltaPhiCharge = new TH2F
+ ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
+ 200,0,120,200,0,6.4);
+ fhDeltaPhiCharge->SetYTitle("#Delta #phi");
+ fhDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhDeltaPhiCharge) ;
+
+ fhDeltaEtaCharge = new TH2F
+ ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
+ 200,0,120,200,-2,2);
+ fhDeltaEtaCharge->SetYTitle("#Delta #eta");
+ fhDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhDeltaEtaCharge) ;
+
+ //
+ if(!fJetsOnlyInCTS){
+ fhPi0Ratio = new TH2F
+ ("Pi0Ratio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
+ 120,0,120,120,0,1);
+ fhPi0Ratio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}");
+ fhPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhPi0Ratio) ;
+
+ fhDeltaPhiPi0 = new TH2F
+ ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
+ 200,0,120,200,0,6.4);
+ fhDeltaPhiPi0->SetYTitle("#Delta #phi");
+ fhDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhDeltaPhiPi0) ;
+
+ fhDeltaEtaPi0 = new TH2F
+ ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
+ 200,0,120,200,-2,2);
+ fhDeltaEtaPi0->SetYTitle("#Delta #eta");
+ fhDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhDeltaEtaPi0) ;
+
+ //
+ fhAnglePair = new TH2F
+ ("AnglePair",
+ "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}",
+ 200,0,50,200,0,0.2);
+ fhAnglePair->SetYTitle("Angle (rad)");
+ fhAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePair) ;
+
+ fhAnglePairAccepted = new TH2F
+ ("AnglePairAccepted",
+ "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window",
+ 200,0,50,200,0,0.2);
+ fhAnglePairAccepted->SetYTitle("Angle (rad)");
+ fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairAccepted) ;
+
+ fhAnglePairNoCut = new TH2F
+ ("AnglePairNoCut",
+ "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2);
+ fhAnglePairNoCut->SetYTitle("Angle (rad)");
+ fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairNoCut) ;
+
+ fhAnglePairLeadingCut = new TH2F
+ ("AnglePairLeadingCut",
+ "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}",
+ 200,0,50,200,0,0.2);
+ fhAnglePairLeadingCut->SetYTitle("Angle (rad)");
+ fhAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairLeadingCut) ;
+
+ fhAnglePairAngleCut = new TH2F
+ ("AnglePairAngleCut",
+ "Angle between all #gamma pair (angle + leading cut) vs p_{T #pi^{0}}"
+ ,200,0,50,200,0,0.2);
+ fhAnglePairAngleCut->SetYTitle("Angle (rad)");
+ fhAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairAngleCut) ;
+
+ fhAnglePairAllCut = new TH2F
+ ("AnglePairAllCut",
+ "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T #pi^{0}}"
+ ,200,0,50,200,0,0.2);
+ fhAnglePairAllCut->SetYTitle("Angle (rad)");
+ fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairAllCut) ;
+
+ fhAnglePairLeading = new TH2F
+ ("AnglePairLeading",
+ "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}",
+ 200,0,50,200,0,0.2);
+ fhAnglePairLeading->SetYTitle("Angle (rad)");
+ fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+ fOutputContainer->Add(fhAnglePairLeading) ;
+
+ //
+ fhInvMassPairNoCut = new TH2F
+ ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
+ 120,0,120,360,0,0.5);
+ fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+ fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhInvMassPairNoCut) ;
+
+ fhInvMassPairLeadingCut = new TH2F
+ ("InvMassPairLeadingCut",
+ "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
+ 120,0,120,360,0,0.5);
+ fhInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+ fhInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhInvMassPairLeadingCut) ;
+
+ fhInvMassPairAngleCut = new TH2F
+ ("InvMassPairAngleCut",
+ "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
+ 120,0,120,360,0,0.5);
+ fhInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+ fhInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhInvMassPairAngleCut) ;
+
+ fhInvMassPairAllCut = new TH2F
+ ("InvMassPairAllCut",
+ "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
+ 120,0,120,360,0,0.5);
+ fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+ fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhInvMassPairAllCut) ;
+
+ fhInvMassPairLeading = new TH2F
+ ("InvMassPairLeading",
+ "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
+ 120,0,120,360,0,0.5);
+ fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
+ fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhInvMassPairLeading) ;
+ }
+
+
+ if(!fSeveralConeAndPtCuts){
+
+ //Count
+ fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
+ fhNBkg->SetYTitle("counts");
+ fhNBkg->SetXTitle("N");
+ fOutputContainer->Add(fhNBkg) ;
+
+ fhNLeading = new TH2F
+ ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120);
+ fhNLeading->SetYTitle("p_{T charge} (GeV/c)");
+ fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
+ fOutputContainer->Add(fhNLeading) ;
+
+ fhNJet = new TH1F("NJet","Accepted jets",240,0,120);
+ fhNJet->SetYTitle("N");
+ fhNJet->SetXTitle("p_{T #gamma}(GeV/c)");
+ fOutputContainer->Add(fhNJet) ;
+
+ //Ratios and Pt dist of reconstructed (not selected) jets
+ //Jet
+ fhJetRatio = new TH2F
+ ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
+ 240,0,120,200,0,10);
+ fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
+ fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhJetRatio) ;
+
+ fhJetPt = new TH2F
+ ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
+ fhJetPt->SetYTitle("p_{T jet}");
+ fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhJetPt) ;
+
+ //Bkg
+
+ fhBkgRatio = new TH2F
+ ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
+ 240,0,120,200,0,10);
+ fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
+ fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhBkgRatio) ;
+
+ fhBkgPt = new TH2F
+ ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
+ fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
+ fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhBkgPt) ;
+
+ //Jet Distributions
+
+ fhJetFragment =
+ new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
+ 240,0.,120.,1000,0.,1.2);
+ fhJetFragment->SetYTitle("x_{T}");
+ fhJetFragment->SetXTitle("p_{T #gamma}");
+ fOutputContainer->Add(fhJetFragment) ;
+
+ fhBkgFragment = new TH2F
+ ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
+ 240,0.,120.,1000,0.,1.2);
+ fhBkgFragment->SetYTitle("x_{T}");
+ fhBkgFragment->SetXTitle("p_{T #gamma}");
+ fOutputContainer->Add(fhBkgFragment) ;
+
+ fhJetPtDist =
+ new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
+ fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhJetPtDist) ;
+
+ fhBkgPtDist = new TH2F
+ ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
+ fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhBkgPtDist) ;
+
+ }
+ else{
+ //If we want to study the jet for different cones and pt
+
+ for(Int_t icone = 0; icone<fNCone; icone++){
+ for(Int_t ipt = 0; ipt<fNPt;ipt++){
+
+ //Jet
+
+ fhJetRatios[icone][ipt] = new TH2F
+ ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
+ +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
+ 240,0,120,200,0,10);
+ fhJetRatios[icone][ipt]->
+ SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
+ fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhJetRatios[icone][ipt]) ;
+
+
+ fhJetPts[icone][ipt] = new TH2F
+ ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
+ +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
+ 240,0,120,400,0,200);
+ fhJetPts[icone][ipt]->
+ SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
+ fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhJetPts[icone][ipt]) ;
+
+ //Bkg
+ fhBkgRatios[icone][ipt] = new TH2F
+ ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
+ +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
+ 240,0,120,200,0,10);
+ fhBkgRatios[icone][ipt]->
+ SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
+ fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhBkgRatios[icone][ipt]) ;
+
+ fhBkgPts[icone][ipt] = new TH2F
+ ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
+ +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
+ 240,0,120,400,0,200);
+ fhBkgPts[icone][ipt]->
+ SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
+ fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhBkgPts[icone][ipt]) ;
+
+ //Counts
+ fhNBkgs[icone][ipt] = new TH1F
+ ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
+ +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
+ fhNBkgs[icone][ipt]->SetYTitle("counts");
+ fhNBkgs[icone][ipt]->SetXTitle("N");
+ fOutputContainer->Add(fhNBkgs[icone][ipt]) ;
+
+ fhNLeadings[icone][ipt] = new TH2F
+ ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>"
+ +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
+ fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
+ fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
+ fOutputContainer->Add(fhNLeadings[icone][ipt]) ;
+
+ fhNJets[icone][ipt] = new TH1F
+ ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "Number of neutral jets, cone ="+fNameCones[icone]+", pt>"
+ +fNamePtThres[ipt]+" GeV/c",120,0,120);
+ fhNJets[icone][ipt]->SetYTitle("N");
+ fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
+ fOutputContainer->Add(fhNJets[icone][ipt]) ;
+
+ //Fragmentation Function
+ fhJetFragments[icone][ipt] = new TH2F
+ ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
+ +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
+ fhJetFragments[icone][ipt]->SetYTitle("x_{T}");
+ fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
+ fOutputContainer->Add(fhJetFragments[icone][ipt]) ;
+
+ fhBkgFragments[icone][ipt] = new TH2F
+ ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
+ +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
+ fhBkgFragments[icone][ipt]->SetYTitle("x_{T}");
+ fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
+ fOutputContainer->Add(fhBkgFragments[icone][ipt]) ;
+
+ //Jet particle distribution
+
+ fhJetPtDists[icone][ipt] = new TH2F
+ ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
+ " GeV/c",120,0.,120.,120,0.,120.);
+ fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhJetPtDists[icone][ipt]) ;
+
+ fhBkgPtDists[icone][ipt] = new TH2F
+ ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+ "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
+ " GeV/c",120,0.,120.,120,0.,120.);
+ fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fhBkgPtDists[icone][ipt]) ;
+
+ }//ipt
+ } //icone
+ }//If we want to study any cone or pt threshold
+}
+
+//____________________________________________________________________________
+void AliAnaGammaJet::MakeJet(TClonesArray * pl, TParticle * pGamma, TParticle* pLeading,TString lastname)
+{
+ //Fill the jet with the particles around the leading particle with
+ //R=fCone and pt_th = fPtThres. Caclulate the energy of the jet and
+ //check if we select it. Fill jet histograms
+
+ TClonesArray * jetList = new TClonesArray("TParticle",1000);
+ TClonesArray * bkgList = new TClonesArray("TParticle",1000);
+
+ TLorentzVector jet (0,0,0,0);
+ TLorentzVector bkg(0,0,0,0);
+ TLorentzVector lv (0,0,0,0);
+
+ Double_t ptjet = 0.0;
+ Double_t ptbkg = 0.0;
+ Int_t n0 = 0;
+ Int_t n1 = 0;
+ Bool_t b1 = kFALSE;
+ Bool_t b0 = kFALSE;
+
+ Double_t ptg = pGamma->Pt();
+ Double_t phig = pGamma->Phi();
+ Double_t ptl = pLeading->Pt();
+ Double_t phil = pLeading->Phi();
+ Double_t etal = pLeading->Eta();
+
+ Float_t ptcut = 0. ;
+ if(fPbPb){
+ if(ptg > fPtJetSelectionCut) ptcut = 2. ;
+ else ptcut = 0.5;
+ }
+
+ TIter next(pl) ;
+ TParticle * particle = 0 ;
+ while ( (particle = dynamic_cast<TParticle*>(next())) ) {
+
+ b0 = kFALSE;
+ b1 = kFALSE;
+
+ //Particles in jet
+ SetJet(particle, b0, fCone, etal, phil) ;
+
+ if(b0){
+ new((*jetList)[n0++]) TParticle(*particle) ;
+ particle->Momentum(lv);
+ if(particle->Pt() > ptcut ){
+ jet+=lv;
+ ptjet+=particle->Pt();
+ }
+ }
+
+ //Background around (phi_gamma-pi, eta_leading)
+ SetJet(particle, b1, fCone,etal, phig) ;
+
+ if(b1) {
+ new((*bkgList)[n1++]) TParticle(*particle) ;
+ particle->Momentum(lv);
+ if(particle->Pt() > ptcut ){
+ bkg+=lv;
+ ptbkg+=particle->Pt();
+ }
+ }
+ }
+
+ ptjet = jet.Pt();
+ ptbkg = bkg.Pt();
+
+ if(ptjet > 0.) {
+
+ AliDebug(2,Form("Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg));
+
+ //Fill histograms
+
+ Double_t ratjet = ptjet/ptg ;
+ Double_t ratbkg = ptbkg/ptg ;
+
+ dynamic_cast<TH2F*>
+ (fOutputContainer->FindObject("JetRatio"+lastname))
+ ->Fill(ptg,ratjet);
+ dynamic_cast<TH2F*>
+ (fOutputContainer->FindObject("JetPt"+lastname))
+ ->Fill(ptg,ptjet);
+
+ dynamic_cast<TH2F*>
+ (fOutputContainer->FindObject("BkgRatio"+lastname))
+ ->Fill(ptg,ratbkg);
+
+ dynamic_cast<TH2F*>
+ (fOutputContainer->FindObject("BkgPt"+lastname))
+ ->Fill(ptg,ptbkg);
+
+
+ //Jet selection
+ Bool_t kSelect = kFALSE;
+ if(fSelect == 0)
+ kSelect = kTRUE; //Accept all jets, no restriction
+ else if(fSelect == 1){
+ //Selection with parametrized cuts
+ if(IsJetSelected(ptg,ptjet)) kSelect = kTRUE;
+ }
+ else if(fSelect == 2){
+ //Simple selection
+ if(!fJetsOnlyInCTS){
+ if((ratjet < fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE;
+ }
+ else{
+ if((ratjet < fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE;
+ }
+ }
+ else
+ AliError("Jet selection option larger than 2, DONT SELECT JETS");
+
+
+ if(kSelect){
+ if (GetPrintInfo())
+ AliInfo(Form("Jet Selected: pt %f ", ptjet)) ;
+
+ FillJetHistos(jetList, ptg, ptl,"Jet",lastname);
+ FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname);
+ }
+ } //ptjet > 0
+
+ jetList ->Delete();
+ bkgList ->Delete();
+
+}
+
+//____________________________________________________________________________
+void AliAnaGammaJet::Print(const Option_t * opt) const
+{
+
+ //Print some relevant parameters set for the analysis
+ if(! opt)
+ return;
+
+ Info("Print", "%s %s", GetName(), GetTitle() ) ;
+ if(!fJetsOnlyInCTS)
+ printf("EMCAL Acceptance cut : phi min %f, phi max %f, eta %f \n", fPhiEMCALCut[0], fPhiEMCALCut[1], fEtaEMCALCut) ;
+
+ printf("Phi_Leading < %f\n", fPhiMaxCut) ;
+ printf("Phi_Leading > %f\n", fPhiMinCut) ;
+ printf("pT Leading / pT Gamma < %f\n", fRatioMaxCut) ;
+ printf("pT Leading / pT Gamma > %f\n", fRatioMinCut) ;
+ printf("M_pair < %f\n", fInvMassMaxCut) ;
+ printf("M_pair > %f\n", fInvMassMinCut) ;
+ if(fSelect == 2){
+ printf("pT Jet / pT Gamma < %f\n", fJetRatioMaxCut) ;
+ printf("pT Jet / pT Gamma > %f\n", fJetRatioMinCut) ;
+ printf("pT Jet (Only CTS)/ pT Gamma < %f\n", fJetCTSRatioMaxCut) ;
+ printf("pT Jet (Only CTS)/ pT Gamma > %f\n", fJetCTSRatioMinCut) ;
+ }
+
+
+}
+
+//___________________________________________________________________
+void AliAnaGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone,
+ Double_t eta, Double_t phi)
+{
+
+ //Check if the particle is inside the cone defined by the leading particle
+ b = kFALSE;
+
+ if(phi > TMath::TwoPi())
+ phi-=TMath::TwoPi();
+ if(phi < 0.)
+ phi+=TMath::TwoPi();
+
+ Double_t rad = 10000 + cone;
+
+ if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
+ rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
+ TMath::Power(part->Phi()-phi,2));
+ else{
+ if(part->Phi()-phi > TMath::TwoPi() - cone)
+ rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
+ TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
+ if(part->Phi()-phi < -(TMath::TwoPi() - cone))
+ rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
+ TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
+ }
+
+ if(rad < cone )
+ b = kTRUE;
+
+}
+
+
+void AliAnaGammaJet::Terminate(Option_t *)
+{
+ // The Terminate() function is the last function to be called during
+ // a query. It always runs on the client, it can be used to present
+ // the results graphically or save the results to file.
+
+
+}
--- /dev/null
+#ifndef ALIANAGAMMAJET_H
+#define ALIANAGAMMAJET_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+// Class for the analysis of gamma-jet correlations.
+// Basically it seaches for a prompt photon in the Calorimeters acceptance,
+// if so we construct a jet around the highest pt particle in the opposite
+// side in azimuth. This jet has to fullfill several conditions to be
+// accepted. Then the fragmentation function of this jet is constructed
+// Class created from old AliPHOSGammaJet
+
+//*-- Author: Gustavo Conesa (INFN-LNF)
+
+// --- ROOT system ---
+#include <TROOT.h>
+#include <TChain.h>
+#include "TTask.h"
+#include "TArrayD.h"
+#include "TChain.h"
+#include <TH2F.h>
+#include <TTree.h>
+#include "AliAnaGammaDirect.h"
+
+class AliESD ;
+
+class AliAnaGammaJet : public AliAnaGammaDirect {
+
+public:
+
+ AliAnaGammaJet(const char *name) ; // default ctor
+ AliAnaGammaJet(const AliAnaGammaJet & gj) ; // cpy ctor
+ virtual ~AliAnaGammaJet() ; //virtual dtor
+ virtual void Exec(Option_t * opt = "") ;
+ virtual void Init(Option_t * opt = "");
+ virtual void Terminate(Option_t * opt = "");
+
+ Bool_t AreJetOnlyInCTS() const {return fJetsOnlyInCTS ; }
+ Double_t GetAngleMaxParam(Int_t i) const {return fAngleMaxParam.At(i) ; }
+ Double_t GetEtaEMCALCut() const {return fEtaEMCALCut;}
+ Double_t GetPhiEMCALCut(Int_t i) const {return fPhiEMCALCut[i];}
+ Double_t GetInvMassMaxCut() const {return fInvMassMaxCut ; }
+ Double_t GetInvMassMinCut() const {return fInvMassMinCut ; }
+ Double_t GetPhiMaxCut() const {return fPhiMaxCut ; }
+ Double_t GetPhiMinCut() const {return fPhiMinCut ; }
+ Double_t GetPtJetSelectionCut() const {return fPtJetSelectionCut ; }
+ Double_t GetJetRatioMaxCut() const {return fJetRatioMaxCut ; }
+ Double_t GetJetRatioMinCut() const {return fJetRatioMinCut ; }
+ Double_t GetRatioMaxCut() const {return fRatioMaxCut ; }
+ Double_t GetRatioMinCut() const {return fRatioMinCut ; }
+ Int_t GetNCones() const {return fNCone ; }
+ Int_t GetNPtThres() const {return fNPt ; }
+ Float_t GetCone() const {return fCone ; }
+ Float_t GetPtThreshold() const {return fPtThreshold ; }
+ Float_t GetCones(Int_t i) const {return fCones[i] ; }
+ Float_t GetPtThreshold(Int_t i) const {return fPtThres[i] ; }
+ TString GetConeName(Int_t i) const {return fNameCones[i] ; }
+ TString GetPtThresName(Int_t i) const {return fNamePtThres[i] ; }
+
+ Bool_t AreSeveralConeAndPtCuts() const {return fSeveralConeAndPtCuts ; }
+ Bool_t IsPbPb() const {return fPbPb ; }
+
+
+ void Print(const Option_t * opt)const;
+
+ void SetAngleMaxParam(Int_t i, Double_t par)
+ {fAngleMaxParam.AddAt(par,i) ; }
+ void SetSeveralConeAndPtCuts(Bool_t several){fSeveralConeAndPtCuts = several ;}
+ void SetEtaEMCALCut(Double_t etacut) {fEtaEMCALCut = etacut ; }
+ void SetPhiEMCALCut(Double_t phi, Int_t i){fPhiEMCALCut[i] = phi; }
+ void SetPbPb(Bool_t opt){fPbPb = opt; }
+ void SetNCones(Int_t n){fNCone = n ; }
+ void SetNPtThresholds(Int_t n){fNPt = n ; }
+ void SetCones(Int_t i, Float_t cone, TString sc)
+ {fCones[i] = cone ; fNameCones[i] = sc; };
+ void SetCone(Float_t cone)
+ {fCone = cone; }
+ void SetPtThreshold(Float_t pt){fPtThreshold = pt; };
+ void SetPtThresholds(Int_t i,Float_t pt, TString spt){fPtThres[i] = pt ;
+ fNamePtThres[i] = spt; };
+ void SetInvMassCutRange(Double_t invmassmin, Double_t invmassmax)
+ {fInvMassMaxCut =invmassmax; fInvMassMinCut =invmassmin;}
+ void SetJetRatioCutRange(Double_t ratiomin, Double_t ratiomax)
+ {fJetRatioMaxCut =ratiomax; fJetRatioMinCut = ratiomin ; }
+ void SetJetsOnlyInCTS(Bool_t opt){fJetsOnlyInCTS = opt; }
+ void SetJetCTSRatioCutRange(Double_t ratiomin, Double_t ratiomax)
+ {fJetCTSRatioMaxCut =ratiomax; fJetCTSRatioMinCut = ratiomin ; }
+ void SetPhiCutRange(Double_t phimin, Double_t phimax)
+ {fPhiMaxCut =phimax; fPhiMinCut =phimin;}
+ void SetPtJetSelectionCut(Double_t cut){fPtJetSelectionCut = cut; }
+ void SetJetSelection(UInt_t select){ fSelect= select ; }
+ void SetRatioCutRange(Double_t ratiomin, Double_t ratiomax)
+ {fRatioMaxCut = ratiomax; fRatioMinCut = ratiomin;}
+
+
+ private:
+
+ Double_t CalculateJetRatioLimit(const Double_t ptg, const Double_t *param,
+ const Double_t *x);
+
+ void FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl,TString type, TString lastname);
+
+ Bool_t IsAngleInWindow(const Float_t angle, const Float_t e);
+ Bool_t IsJetSelected(const Double_t ptg, const Double_t ptjet);
+
+ void MakeJet(TClonesArray * particleList, TParticle * pGamma, TParticle* pLeading, TString lastname);
+
+ void GetLeadingCharge(TClonesArray * pl, TParticle *pGamma, TParticle * pLeading) const ;
+ void GetLeadingPi0 (TClonesArray * pl, TParticle *pGamma, TParticle * pLeading) ;
+ Bool_t GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe, TParticle *pGamma, TParticle * pLeading) ;
+ void MakeHistos() ;
+
+ void SetJet(TParticle * part, Bool_t & b, Float_t cone, Double_t eta,
+ Double_t phi);
+
+ private:
+
+ // TTree *fChain ; //!pointer to the analyzed TTree or TChain
+ //AliESD *fESD ; //! Declaration of leave types
+
+ Bool_t fSeveralConeAndPtCuts; // To play with the jet cone size and pt th.
+ //Bool_t fPrintInfo ; //Print most interesting information on screen
+ // and give interesting information
+
+ Bool_t fPbPb; // PbPb event
+ Bool_t fJetsOnlyInCTS ; // Jets measured only in TPC+ITS.
+ Double_t fEtaEMCALCut ; // Eta EMCAL acceptance
+ Double_t fPhiEMCALCut[2] ; // Phi cut maximum
+ Double_t fPhiMaxCut ; // Phi EMCAL maximum acceptance
+ Double_t fPhiMinCut ; // Phi EMCAL minimum acceptance
+ // Double_t fGammaPtCut ; // Min pt in Calorimeter
+ Double_t fInvMassMaxCut ; // Invariant Mass cut maximum
+ Double_t fInvMassMinCut ; // Invariant Masscut minimun
+ Double_t fRatioMaxCut ; // Leading particle/gamma Ratio cut maximum
+ Double_t fRatioMinCut ; // Leading particle/gamma Ratio cut minimum
+
+ //Jet selection parameters
+ //Fixed cuts (old)
+ Double_t fJetCTSRatioMaxCut ; // Leading particle/gamma Ratio cut maximum
+ Double_t fJetCTSRatioMinCut ; // Leading particle/gamma Ratio cut minimum
+ Double_t fJetRatioMaxCut ; // Jet/gamma Ratio cut maximum
+ Double_t fJetRatioMinCut ; // Jet/gamma Ratio cut minimum
+
+ //Cuts depending on jet pt
+ Double_t fJetE1[2]; //Rec. jet energy parameters
+ Double_t fJetE2[2]; //Rec. jet energy parameters
+ Double_t fJetSigma1[2];//Rec. sigma of jet energy parameters
+ Double_t fJetSigma2[2];//Rec. sigma of jet energy parameters
+ Double_t fBkgMean[6]; //Background mean energy
+ Double_t fBkgRMS[6]; //Background RMS
+ Double_t fJetXMin1[6]; //X Factor to set jet min limit for pp
+ Double_t fJetXMin2[6]; //X Factor to set jet min limit for PbPb
+ Double_t fJetXMax1[6]; //X Factor to set jet max limit for pp
+ Double_t fJetXMax2[6]; //X Factor to set jet max limit for PbPb
+
+ Int_t fNCone ; // Number of jet cones sizes
+ Int_t fNPt ; // Number of jet particle pT threshold
+ Double_t fCone ; // Jet cone sizes under study (!fSeveralConeAndPtCuts)
+ Double_t fCones[10]; // Jet cone sizes under study (fSeveralConeAndPtCuts)
+ TString fNameCones[10]; // String name of cone to append to histos
+ Double_t fPtThreshold; // Jet pT threshold under study(!fSeveralConeAndPtCuts)
+ Double_t fPtThres[10]; // Jet pT threshold under study(fSeveralConeAndPtCuts)
+ Double_t fPtJetSelectionCut; // Jet pt to change to low pt jets analysis
+ TObjArray *fOutputContainer ; //! output data container
+ TString fNamePtThres[10]; // String name of pt th to append to histos
+ TArrayD fAngleMaxParam ; //Max opening angle selection parameters
+ UInt_t fSelect ; //kTRUE: Selects all jets, no limits.
+ //TString fCalorimeter ; //PHOS or EMCAL detects Gamma
+ //Bool_t fEMCALPID ; //Fill EMCAL particle lists with particles with corresponding pid
+ //Bool_t fPHOSPID; //Fill PHOS particle lists with particles with corresponding pid
+
+ //Histograms
+
+ TH2F * fhChargeRatio ;
+ TH2F * fhPi0Ratio ;
+ TH2F * fhDeltaPhiCharge ;
+ TH2F * fhDeltaPhiPi0 ;
+ TH2F * fhDeltaEtaCharge ;
+ TH2F * fhDeltaEtaPi0 ;
+ TH2F * fhAnglePair ;
+ TH2F * fhAnglePairAccepted ;
+ TH2F * fhAnglePairNoCut ;
+ TH2F * fhAnglePairLeadingCut ;
+ TH2F * fhAnglePairAngleCut ;
+ TH2F * fhAnglePairAllCut ;
+ TH2F * fhAnglePairLeading ;
+ TH2F * fhInvMassPairNoCut ;
+ TH2F * fhInvMassPairLeadingCut ;
+ TH2F * fhInvMassPairAngleCut ;
+ TH2F * fhInvMassPairAllCut ;
+ TH2F * fhInvMassPairLeading ;
+ TH1F * fhNBkg ;
+ TH2F * fhNLeading ;
+ TH1F * fhNJet ;
+ TH2F * fhJetRatio ;
+ TH2F * fhJetPt ;
+ TH2F * fhBkgRatio ;
+ TH2F * fhBkgPt ;
+ TH2F * fhJetFragment ;
+ TH2F * fhBkgFragment ;
+ TH2F * fhJetPtDist ;
+ TH2F * fhBkgPtDist ;
+
+ TH2F * fhJetRatios[5][5];
+ TH2F * fhJetPts[5][5];
+ TH2F * fhBkgRatios[5][5];
+ TH2F * fhBkgPts[5][5];
+
+ TH2F * fhNLeadings[5][5];
+ TH1F * fhNJets[5][5];
+ TH1F * fhNBkgs[5][5];
+
+ TH2F * fhJetFragments[5][5];
+ TH2F * fhBkgFragments[5][5];
+ TH2F * fhJetPtDists[5][5];
+ TH2F * fhBkgPtDists[5][5];
+
+ ClassDef(AliAnaGammaJet,0)
+} ;
+
+
+#endif //ALIANAGAMMAJET_H
+
+
+
--- /dev/null
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliAnaGammaDirect+;
+#pragma link C++ class AliAnaGammaHadron+;
+#pragma link C++ class AliAnaGammaJet+;
+
+#endif
--- /dev/null
+
+include $(ROOTSYS)/test/Makefile.arch
+
+default-target: libGamma.so
+
+ALICEINC = -I.
+
+### define include dir for local case and par case
+ifneq ($(ANALYSIS_NEW_INCLUDE),)
+ ALICEINC += -I../$(ESD_INCLUDE) -I../$(ANALYSIS_NEW_INCLUDE)
+else
+ ifneq ($(ALICE_ROOT),)
+ ALICEINC += -I$(ALICE_ROOT)/include
+ endif
+endif
+
+# for building of Gamma.par
+ifneq ($(GAMMA_INCLUDE),)
+ ALICEINC += -I../$(GAMMA_INCLUDE)
+endif
+
+CXXFLAGS += $(ALICEINC) -g
+
+PACKAGE = Gamma
+include lib$(PACKAGE).pkg
+
+DHDR_Gamma := $(DHDR)
+HDRS_Gamma := $(HDRS)
+SRCS_Gamma := $(SRCS) G__$(PACKAGE).cxx
+OBJS_Gamma := $(SRCS_Gamma:.cxx=.o)
+
+PARFILE = $(PACKAGE).par
+
+
+lib$(PACKAGE).so: $(OBJS_Gamma)
+ @echo "Linking" $@ ...
+ @/bin/rm -f $@
+ifeq ($(PLATFORM),macosx)
+ @$(LD) -bundle -undefined $(UNDEFOPT) $(LDFLAGS) $^ -o $@
+else
+ @$(LD) $(SOFLAGS) $(LDFLAGS) $^ -o $@
+endif
+ @chmod a+x $@
+ @echo "done"
+
+%.o: %.cxx %.h
+ $(CXX) $(CXXFLAGS) -c $< -o $@
+
+clean:
+ @rm -f $(OBJS_Gamma) *.so G__$(PACKAGE).* $(PARFILE)
+
+G__$(PACKAGE).cxx G__$(PACKAGE).h: $(HDRS) $(DHDR)
+ @echo "Generating dictionary ..."
+ rootcint -f $@ -c $(ALICEINC) $^
+
+### CREATE PAR FILE
+
+$(PARFILE): $(patsubst %,$(PACKAGE)/%,$(filter-out G__%, $(HDRS_Gamma) $(SRCS_Gamma) $(DHDR_Gamma) Makefile Makefile.arch lib$(PACKAGE).pkg PROOF-INF))
+ @echo "Creating archive" $@ ...
+ @tar cfzh $@ $(PACKAGE)
+ @rm -rf $(PACKAGE)
+ @echo "done"
+
+$(PACKAGE)/Makefile: Makefile #.$(PACKAGE)
+ @echo Copying $< to $@ with transformations
+ @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+ @sed 's/include \$$(ROOTSYS)\/test\/Makefile.arch/include Makefile.arch/' < $^ > $@
+
+$(PACKAGE)/Makefile.arch: $(ROOTSYS)/test/Makefile.arch
+ @echo Copying $< to $@
+ @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+ @cp -a $^ $@
+
+$(PACKAGE)/PROOF-INF: PROOF-INF.$(PACKAGE)
+ @echo Copying $< to $@
+ @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+ @cp -a -r $^ $@
+
+$(PACKAGE)/%: %
+ @echo Copying $< to $@
+ @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+ @cp -a $< $@
+
+test-%.par: %.par
+ @echo "INFO: The file $< is now tested, in case of an error check in par-tmp."
+ @mkdir -p par-tmp
+ @cd par-tmp; tar xfz ../$<; cd $(subst .par,,$<); PROOF-INF/BUILD.sh
+ @rm -rf par-tmp
+ @echo "INFO: Testing succeeded (already cleaned up)"
--- /dev/null
+SRCS = AliAnaGammaDirect.cxx AliAnaGammaHadron.cxx AliAnaGammaJet.cxx
+
+HDRS:= $(SRCS:.cxx=.h)
+
+DHDR:= GammaLinkDef.h
+
+EXPORT:=$(SRCS:.cxx=.h)
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff