]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/AliAnaGammaDirect.cxx
coding conventions, compilation warnings, code cleanup
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaDirect.cxx
index 82d4f8d73f527799af8c83fbaec3368456d1e575..6611cff720e5d5e16f931a4fafa81c281537ca72 100644 (file)
 /* History of cvs commits:
  *
  * $Log$
- * Revision 1.3  2007/03/08 10:24:32  schutz
- * Coding convention
- *
- * Revision 1.2  2007/02/09 18:40:40  schutz
- * B\bNew version from Gustavo
- *
- * Revision 1.1  2007/01/23 17:17:29  schutz
- * New Gamma package
+ * Revision 1.4.4.4  2007/07/26 10:32:09  schutz
+ * new analysis classes in the the new analysis framework
  *
  *
  */
 
 //_________________________________________________________________________
-// 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 for the prompt gamma analysis, isolation cut
 //
 //  Class created from old AliPHOSGammaJet 
 //  (see AliRoot versions previous Release 4-09)
 //
 //*-- Author: Gustavo Conesa (LNF-INFN) 
 //////////////////////////////////////////////////////////////////////////////
-
-
-// --- ROOT system ---
-
-#include <TFile.h>
+  
+  
+// --- ROOT system --- 
 #include <TParticle.h>
 #include <TH2.h>
-#include <TChain.h>
+#include <TList.h>
 #include "AliAnaGammaDirect.h" 
-#include "AliESD.h"
-#include "AliESDtrack.h"
-#include "AliESDCaloCluster.h"
 #include "Riostream.h"
 #include "AliLog.h"
-
+  
 ClassImp(AliAnaGammaDirect)
-
+  
 //____________________________________________________________________________
-  AliAnaGammaDirect::AliAnaGammaDirect(const char *name) : 
-    AliAnalysisTask(name,""), fChain(0), fESD(0),
-    fOutputContainer(new TObjArray(100)), 
-    fPrintInfo(0), fMinGammaPt(0.),
-    fCalorimeter(""), fEMCALPID(0),fPHOSPID(0),
-    fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.), fPHOSPhotonWeight(0.),
+  AliAnaGammaDirect::AliAnaGammaDirect() : 
+    TObject(),
+    fMinGammaPt(0.),
     fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
-    fMakeICMethod(0),fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0)
+    fICMethod(0),fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0),  
+    //kSeveralIC
+    fNCones(0),fNPtThres(0), fConeSizes(),  fPtThresholds(), 
+    fhPtThresIsolated(), fhPtSumIsolated()
 {
-  //Ctor
-        
+  //default ctor
+  
   //Initialize parameters
   InitParameters();
 
-  // 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),
-  fEMCALPhotonWeight(g.fEMCALPhotonWeight), 
-  fEMCALPi0Weight(g.fEMCALPi0Weight), 
-  fPHOSPhotonWeight(g.fPHOSPhotonWeight),
+  TObject(g),
+  fMinGammaPt(g.fMinGammaPt), 
   fConeSize(g.fConeSize),
   fPtThreshold(g.fPtThreshold),
   fPtSumThreshold(g.fPtSumThreshold), 
-  fMakeICMethod(g.fMakeICMethod),
-  fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),fhEtaGamma(g.fhEtaGamma)
+  fICMethod(g.fICMethod),
+  fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),fhEtaGamma(g.fhEtaGamma),  
+  //kSeveralIC
+  fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(), 
+  fhPtThresIsolated(), fhPtSumIsolated()
 {
   // cpy ctor
-  SetName (g.GetName()) ; 
-  SetTitle(g.GetTitle()) ; 
+  
+  //kSeveralIC
+  for(Int_t i = 0; i < fNCones ; i++){
+    fConeSizes[i] =  g.fConeSizes[i];
+    fhPtSumIsolated[i] = g.fhPtSumIsolated[i]; 
+    for(Int_t j = 0; j < fNPtThres ; j++)
+      fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j]; 
+  }
+  
+  for(Int_t i = 0; i < fNPtThres ; i++)
+    fPtThresholds[i]=   g.fPtThresholds[i];
 }
 
 //_________________________________________________________________________
 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & source)
 {
   // assignment operator
-
+  
   if(&source == this) return *this;
 
-  fChain = source.fChain ; 
-  fESD = source.fESD ;
-  fOutputContainer = source.fOutputContainer ;
-  fPrintInfo = source.fPrintInfo ;
   fMinGammaPt = source.fMinGammaPt ;   
-  fCalorimeter = source. fCalorimeter ;  
-  fEMCALPID = source.fEMCALPID ;
-  fPHOSPID = source.fPHOSPID ;
-  fEMCALPhotonWeight = source. fEMCALPhotonWeight ;
-  fEMCALPi0Weight = source.fEMCALPi0Weight ;
-  fPHOSPhotonWeight = source.fPHOSPhotonWeight ;
   fConeSize = source.fConeSize ;
   fPtThreshold = source.fPtThreshold ;
   fPtSumThreshold = source.fPtSumThreshold ; 
-  fMakeICMethod = source.fMakeICMethod ;
+  fICMethod = source.fICMethod ;
   fhNGamma = source.fhNGamma ; 
   fhPhiGamma = source.fhPhiGamma ;
   fhEtaGamma = source.fhEtaGamma ;
-
-
+  
+  //kSeveralIC
+  fNCones = source.fNCones ;
+  fNPtThres = source.fNPtThres ; 
+   
+  for(Int_t i = 0; i < fNCones ; i++){
+    fConeSizes[i] =  source.fConeSizes[i];
+    fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
+    for(Int_t j = 0; j < fNPtThres ; j++)
+      fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
+  }
+  
+  for(Int_t i = 0; i < fNPtThres ; i++)
+    fPtThresholds[i]=   source.fPtThresholds[i];
+  
   return *this;
-
+  
 }
 
 //____________________________________________________________________________
 AliAnaGammaDirect::~AliAnaGammaDirect() 
 {
   // Remove all pointers
-  fOutputContainer->Clear() ; 
-  delete fOutputContainer ;
+  
   delete fhNGamma    ;  
   delete fhPhiGamma  ; 
   delete fhEtaGamma   ;  
-
-}
-
-//______________________________________________________________________________
-void AliAnaGammaDirect::ConnectInputData(const Option_t*)
-{
-  // Initialisation of branch container and histograms 
-    
-  AliInfo(Form("*** Initialization of %s", GetName())) ; 
   
-  // Get input data
-  fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
-  if (!fChain) {
-    AliError(Form("Input 0 for %s not found\n", GetName()));
-    return ;
-  }
+  //kSeveralIC
+  delete [] fhPtThresIsolated ;
+  delete [] fhPtSumIsolated ;
   
-  // 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);
-  } else {
-    fESD = new AliESD();
-    SetBranchAddress(0, "ESD", &fESD);
-  }
-
 }
 
 //________________________________________________________________________
-void AliAnaGammaDirect::CreateOutputObjects()
+TList *  AliAnaGammaDirect::GetCreateOutputObjects()
 {  
 
   // Create histograms to be saved in output file and 
-  // store them in fOutputContainer
-
-  fOutputContainer = new TObjArray(3) ;
-
+  // store them in outputContainer
+  TList * outputContainer = new TList() ; 
+  outputContainer->SetName("DirectGammaHistos") ; 
+  
   //Histograms of highest gamma identified in Event
-  fhNGamma  = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120); 
+  fhNGamma  = new TH1F("NGamma","Number of #gamma over calorimeter",240,0,120); 
   fhNGamma->SetYTitle("N");
   fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
-  fOutputContainer->Add(fhNGamma) ; 
+  outputContainer->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) ; 
+  outputContainer->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) ;
-
-}
-
-//____________________________________________________________________________
-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
-       TLorentzVector momentum ;
-       clus->GetMomentum(momentum);
-       TParticle * particle = new TParticle() ;
-       //particle->SetMomentum(px,py,pz,en) ;
-       particle->SetMomentum(momentum) ;
-
-       AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
-       
-       //Select only photons
-       
-       pid=clus->GetPid();
-       //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
-       if( !fPHOSPID)
-         new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
-       else if( pid[AliPID::kPhoton] > fPHOSPhotonWeight)
-         new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
-      }
-    }
+  outputContainer->Add(fhEtaGamma) ;
+
+  if(fICMethod == kSeveralIC){
+    char name[128];
+    char title[128];
+    for(Int_t icone = 0; icone<fNCones; icone++){
+      sprintf(name,"PtSumIsolated_Cone_%d",icone);
+      sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+      fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
+      fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+      fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
+      outputContainer->Add(fhPtSumIsolated[icone]) ; 
+      
+      for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
+       sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
+       sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+       fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
+       fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+       outputContainer->Add(fhPtThresIsolated[icone][ipt]) ; 
+      }//icone loop
+    }//ipt loop
   }
 
-  //########## #####################
-  //Prepare bool array for EMCAL track matching
-
-  // step 1, set the flag in a way that it rejects all not-V1 clusters
-  // but at this point all V1 clusters are accepted
-  
-  Int_t begem = fESD->GetFirstEMCALCluster();  
-  Int_t endem = fESD->GetFirstEMCALCluster() + 
-    fESD->GetNumberOfEMCALClusters() ;  
-//   if(endem < begem+12)
-//     AliError("Number of pseudoclusters smaller than 12");
-  Bool_t *useCluster = new Bool_t[endem+1];
-  
-//   for (npar =  0; npar <  endem; npar++){
-//     if(npar < begem+12)
-//       useCluster[npar] =kFALSE; //EMCAL Pseudoclusters and PHOS clusters
-//     else
-//       useCluster[npar] =kTRUE;   //EMCAL clusters 
-//   }
-  for (npar =  0; npar <  endem; npar++)
-    useCluster[npar] =kFALSE; //EMCAL Pseudoclusters and clusters
-  
-  //########### 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));
-
-  Int_t iemcalMatch  = -1 ;
-  for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
-    AliESDtrack * track = fESD->GetTrack(npar) ; // retrieve track from esd
+  return outputContainer ;
 
-    // step 2 for EMCAL matching, change the flag for all matched clusters found in tracks
-    iemcalMatch = track->GetEMCALcluster(); 
-    if(iemcalMatch > 0) useCluster[iemcalMatch] = kTRUE; // reject matched cluster
-    
-    //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 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 && !useCluster[npar] ){
-       TLorentzVector momentum ;
-       clus->GetMomentum(momentum);
-       TParticle * particle = new TParticle() ;
-       //particle->SetMomentum(px,py,pz,en) ;
-       particle->SetMomentum(momentum) ;
-       cout<<"GOOD EMCAL "<<particle->Pt()<<endl;
-       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] > fEMCALPhotonWeight)
-             new((*plEMCAL)[indexNe++])       TParticle(*particle) ;       
-         }
-       else
-         {
-             Int_t pdg = 0;
-             if(fEMCALPID) 
-               {
-                 if( pid[AliPID::kPhoton] > fEMCALPhotonWeight) 
-                   pdg = 22;
-                 else if( pid[AliPID::kPi0] > fEMCALPi0Weight)
-                   pdg = 111;
-               }
-             else
-               pdg = 22; //No PID, assume all photons
-             
-             TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
-                                                  momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
-             AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
-             
-             new((*plEMCAL)[indexNe++])       TParticle(*particle) ; 
-             new((*pl)[index++])           TParticle(*particle) ;
-         }
-      }
-    }
-    
-    AliDebug(3,"Particle lists filled");
-    
 }
 
-
-
 //____________________________________________________________________________
-void AliAnaGammaDirect::Exec(Option_t *) 
-{
-  
-  // Processing of one event
-
-  //Get ESDs
-  Long64_t entry = fChain->GetReadEntry() ;
-  
-  if (!fESD) {
-    AliError("fESD is not connected to the input!") ; 
-    return ; 
-  }
-  
-  if (fPrintInfo) 
-     AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ; 
-
-  //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
-
-  TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
-  TClonesArray * plCTS         = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
-  TClonesArray * plNe         = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter (EMCAL)
-  TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured in PHOS as Gamma calorimeter
-  TClonesArray * plEMCAL   = new TClonesArray("TParticle",1000);  // All particles measured in EMCAL as Gamma calorimeter
-
-  TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
-  //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;
-  
-  
-  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
-  
-  AliDebug(2, "End of analysis, delete pointers");
-  
-  particleList->Delete() ; 
-  plCTS->Delete() ;
-  plNe->Delete() ;
-  plEMCAL->Delete() ;
-  plPHOS->Delete() ;
-  pGamma->Delete();
-  PostData(0, fOutputContainer);
-
- delete plNe ;
-  delete plCTS ;
-  //delete plPHOS ;
-  //delete plEMCAL ;
-  delete particleList ;
-
-  //  delete pGamma;
-
-}    
-
-
-//____________________________________________________________________________
-void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &Is) const 
+void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &found) const 
 {
   //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
 
@@ -466,30 +202,29 @@ void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS,
       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;
+      found  = kTRUE;
     }
   }
 
   //Do Isolation?
-  if(fMakeICMethod && Is)
+  if( ( fICMethod == kPtIC  ||  fICMethod == kSumPtIC )  && found)
     {
       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(fICMethod == kPtIC) //Pt thres method
+       found = icPtThres ;
+      if(fICMethod == kSumPtIC) //Pt cone sum method
+       found = 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());
+  if(found){
+    AliDebug(1,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
+    AliDebug(1,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());
   }
@@ -502,21 +237,21 @@ void AliAnaGammaDirect::InitParameters()
 {
  
   //Initialize the parameters of the analysis.
-  fCalorimeter="PHOS";
-  fPrintInfo           = kTRUE;
   fMinGammaPt  = 5. ;
 
   //Fill particle lists when PID is ok
-  fEMCALPID = kFALSE;
-  fPHOSPID = kFALSE;
-  fEMCALPhotonWeight = 0.5 ;
-  fEMCALPi0Weight = 0.5 ;
-  fPHOSPhotonWeight = 0.8 ;
   fConeSize             = 0.2 ; 
   fPtThreshold         = 2.0; 
   fPtSumThreshold  = 1.; 
 
-  fMakeICMethod = 1; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
+  fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
+
+ //-----------kSeveralIC-----------------
+  fNCones           = 4 ; 
+  fNPtThres         = 4 ; 
+  fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
+  fPtThresholds[0]=1.; fPtThresholds[1]=2.; fPtThresholds[2]=3.; fPtThresholds[3]=4.;
+
 }
 
 //__________________________________________________________________
@@ -535,7 +270,7 @@ void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
   Float_t phi    = -100.  ;
   Float_t rad   = -100 ;
   Int_t    n        = 0 ;
-  TParticle * particle  = new TParticle();
+  TParticle * particle  = new TParticle;
 
   coneptsum = 0; 
   icmpt = kFALSE;
@@ -584,27 +319,71 @@ void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
 
 }
 
-void AliAnaGammaDirect::Print(const Option_t * opt) const
+//__________________________________________________________________
+void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS) 
 {
+  //Isolation Cut Analysis for both methods and different pt cuts and cones
 
+  if (fICMethod != kSeveralIC)
+    AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
+  
+  for(Int_t ipr = 0; ipr < plCalo->GetEntries() ; ipr ++ ){
+    TParticle * pCandidate = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
+    
+    if(pCandidate->Pt() > fMinGammaPt){
+      
+      Bool_t  icPtThres   = kFALSE;
+      Bool_t  icPtSum     = kFALSE;
+      
+      Float_t ptC      = pCandidate->Pt() ;
+   
+      fhNGamma->Fill(ptC);
+      fhPhiGamma->Fill( ptC,pCandidate->Phi());
+      fhEtaGamma->Fill(ptC,pCandidate->Eta());
+    
+      for(Int_t icone = 0; icone<fNCones; icone++){
+       fConeSize=fConeSizes[icone] ;
+       Float_t coneptsum = 0 ;
+       for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
+         fPtThreshold=fPtThresholds[ipt] ;
+         MakeIsolationCut(plCTS,plCalo, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
+         AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
+                         pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
+
+         fhPtThresIsolated[icone][ipt]->Fill(ptC); 
+       }//pt thresh loop
+       fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
+      }//cone size loop
+    }//min pt candidate
+  }//candidate loop
+}
+
+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("IC method               =     %d\n", fMakeICMethod) ; 
+  
+  printf("Min Gamma pT      =     %f\n",  fMinGammaPt) ;
+  printf("IC method               =     %d\n", fICMethod) ; 
   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.
-    
+   if(fICMethod == kPtIC) printf("pT threshold           =     %f\n", fPtThreshold) ;
+   if(fICMethod == kSumPtIC) printf("pT sum threshold   =     %f\n", fPtSumThreshold) ;
+   
+  if(fICMethod == kSeveralIC){
+    printf("N Cone Sizes               =     %d\n", fNCones) ; 
+    printf("N pT thresholds           =     %d\n", fNPtThres) ;
+    printf("Cone Sizes                  =    \n") ;
+    for(Int_t i = 0; i < fNCones; i++)
+      printf("   %f;",  fConeSizes[i]) ;
+    printf("    \n") ;
+    for(Int_t i = 0; i < fNPtThres; i++)
+      printf("   %f;",  fPtThresholds[i]) ;
+  }
 
-}
+  printf("    \n") ;
+  
+}