From f9cea31c748bdc44068bda92a902021cb420e707 Mon Sep 17 00:00:00 2001 From: schutz Date: Tue, 23 Jan 2007 17:17:29 +0000 Subject: [PATCH] New Gamma package --- ESDCheck/PROOF-INF.AnalysisCheck/SETUP.C | 8 +- PWG4/AliAnaGammaDirect.cxx | 642 ++++++++++ PWG4/AliAnaGammaDirect.h | 132 +++ PWG4/AliAnaGammaHadron.cxx | 600 ++++++++++ PWG4/AliAnaGammaHadron.h | 120 ++ PWG4/AliAnaGammaJet.cxx | 1354 ++++++++++++++++++++++ PWG4/AliAnaGammaJet.h | 233 ++++ PWG4/GammaLinkDef.h | 11 + PWG4/Makefile | 89 ++ PWG4/libGamma.pkg | 9 + 10 files changed, 3194 insertions(+), 4 deletions(-) create mode 100644 PWG4/AliAnaGammaDirect.cxx create mode 100644 PWG4/AliAnaGammaDirect.h create mode 100644 PWG4/AliAnaGammaHadron.cxx create mode 100644 PWG4/AliAnaGammaHadron.h create mode 100644 PWG4/AliAnaGammaJet.cxx create mode 100644 PWG4/AliAnaGammaJet.h create mode 100644 PWG4/GammaLinkDef.h create mode 100644 PWG4/Makefile create mode 100644 PWG4/libGamma.pkg diff --git a/ESDCheck/PROOF-INF.AnalysisCheck/SETUP.C b/ESDCheck/PROOF-INF.AnalysisCheck/SETUP.C index d63773cc17a..028db33d994 100644 --- a/ESDCheck/PROOF-INF.AnalysisCheck/SETUP.C +++ b/ESDCheck/PROOF-INF.AnalysisCheck/SETUP.C @@ -2,12 +2,12 @@ void SETUP() { // 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"); } diff --git a/PWG4/AliAnaGammaDirect.cxx b/PWG4/AliAnaGammaDirect.cxx new file mode 100644 index 00000000000..7c6ee3bf5b2 --- /dev/null +++ b/PWG4/AliAnaGammaDirect.cxx @@ -0,0 +1,642 @@ +/************************************************************************** + * 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 +#include +#include + +#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(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 "< 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(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(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; iconePt(), 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(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(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(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 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(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 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; iconeSetYTitle("#Sigma p_{T} (GeV/c)"); + fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)"); + fOutputContainer->Add(fhPtSumIsolated[icone]) ; + + for(Int_t ipt = 0; iptSetXTitle("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. + + +} diff --git a/PWG4/AliAnaGammaDirect.h b/PWG4/AliAnaGammaDirect.h new file mode 100644 index 00000000000..65a983ded34 --- /dev/null +++ b/PWG4/AliAnaGammaDirect.h @@ -0,0 +1,132 @@ +#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 +#include +#include "TTask.h" +#include "TArrayD.h" +#include "TChain.h" +#include +#include +#include +#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 + + + diff --git a/PWG4/AliAnaGammaHadron.cxx b/PWG4/AliAnaGammaHadron.cxx new file mode 100644 index 00000000000..8b0049d9fb9 --- /dev/null +++ b/PWG4/AliAnaGammaHadron.cxx @@ -0,0 +1,600 @@ +/************************************************************************** + * 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 +#include +#include + +#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(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(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(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) 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)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) && (invmassFill(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 + }//iprimary0.) + min = TMath::ACos(arg); + + if((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. + + +} diff --git a/PWG4/AliAnaGammaHadron.h b/PWG4/AliAnaGammaHadron.h new file mode 100644 index 00000000000..cf334715a89 --- /dev/null +++ b/PWG4/AliAnaGammaHadron.h @@ -0,0 +1,120 @@ +#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 +#include +#include "TTask.h" +#include "TArrayD.h" +#include "TChain.h" +#include +#include +#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 + + + diff --git a/PWG4/AliAnaGammaJet.cxx b/PWG4/AliAnaGammaJet.cxx new file mode 100644 index 00000000000..66a4c3a184b --- /dev/null +++ b/PWG4/AliAnaGammaJet.cxx @@ -0,0 +1,1354 @@ +/************************************************************************** + * 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 +#include +#include + +#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(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(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; iconeDelete() ; + 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(next())) ) { + ipr++ ; + Double_t pt = particle->Pt(); + + charge = TDatabasePDG::Instance() + ->GetParticle(particle->GetPdgCode())->Charge(); + if(charge != 0){//Only jet Charged particles + dynamic_cast + (fOutputContainer->FindObject(type+"Fragment"+lastname)) + ->Fill(ptg,pt/ptg); + dynamic_cast + (fOutputContainer->FindObject(type+"PtDist"+lastname)) + ->Fill(ptg,pt); + }//charged + + }//while + + if(type == "Bkg") + dynamic_cast + (fOutputContainer->FindObject("NBkg"+lastname)) + ->Fill(ipr); + else{ + dynamic_cast + (fOutputContainer->FindObject("NJet"+lastname))-> + Fill(ptg); + dynamic_cast + (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(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) 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)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)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) && (invmassFill(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 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=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" +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(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 + (fOutputContainer->FindObject("JetRatio"+lastname)) + ->Fill(ptg,ratjet); + dynamic_cast + (fOutputContainer->FindObject("JetPt"+lastname)) + ->Fill(ptg,ptjet); + + dynamic_cast + (fOutputContainer->FindObject("BkgRatio"+lastname)) + ->Fill(ptg,ratbkg); + + dynamic_cast + (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. + + +} diff --git a/PWG4/AliAnaGammaJet.h b/PWG4/AliAnaGammaJet.h new file mode 100644 index 00000000000..b6714c2b669 --- /dev/null +++ b/PWG4/AliAnaGammaJet.h @@ -0,0 +1,233 @@ +#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 +#include +#include "TTask.h" +#include "TArrayD.h" +#include "TChain.h" +#include +#include +#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 + + + diff --git a/PWG4/GammaLinkDef.h b/PWG4/GammaLinkDef.h new file mode 100644 index 00000000000..70fa9ca14c4 --- /dev/null +++ b/PWG4/GammaLinkDef.h @@ -0,0 +1,11 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class AliAnaGammaDirect+; +#pragma link C++ class AliAnaGammaHadron+; +#pragma link C++ class AliAnaGammaJet+; + +#endif diff --git a/PWG4/Makefile b/PWG4/Makefile new file mode 100644 index 00000000000..d8961562b27 --- /dev/null +++ b/PWG4/Makefile @@ -0,0 +1,89 @@ + +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)" diff --git a/PWG4/libGamma.pkg b/PWG4/libGamma.pkg new file mode 100644 index 00000000000..95caba9449a --- /dev/null +++ b/PWG4/libGamma.pkg @@ -0,0 +1,9 @@ +SRCS = AliAnaGammaDirect.cxx AliAnaGammaHadron.cxx AliAnaGammaJet.cxx + +HDRS:= $(SRCS:.cxx=.h) + +DHDR:= GammaLinkDef.h + +EXPORT:=$(SRCS:.cxx=.h) + + -- 2.31.1