PWGGA/CaloTrackCorrelations/AliAnaParticleIsolation.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes hereby granted      *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 // Class for analysis of particle isolation
  18 // Input is selected particles put in AOD branch (AliAODPWG4ParticleCorrelation)
  19 //
  20 //  Class created from old AliPHOSGammaJet
  21 //  (see AliRoot versions previous Release 4-09)
  22 //
  23 // -- Author: Gustavo Conesa (LNF-INFN)
  24
  25 //-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
  26 //////////////////////////////////////////////////////////////////////////////
  27
  28
  29 // --- ROOT system ---
  30 #include <TClonesArray.h>
  31 #include <TList.h>
  32 #include <TObjString.h>
  33 #include <TH2F.h>
  34 #include <TClass.h>
  35
  36 // --- Analysis system ---
  37 #include "AliAnaParticleIsolation.h"
  38 #include "AliCaloTrackReader.h"
  39 #include "AliIsolationCut.h"
  40 #include "AliNeutralMesonSelection.h"
  41 #include "AliAODPWG4ParticleCorrelation.h"
  42 #include "AliMCAnalysisUtils.h"
  43 #include "AliVTrack.h"
  44 #include "AliVCluster.h"
  45
  46 ClassImp(AliAnaParticleIsolation)
  47
  48 //______________________________________________________________________________
  49 AliAnaParticleIsolation::AliAnaParticleIsolation() :
  50 AliAnaCaloTrackCorrBaseClass(),   fCalorimeter(""),
  51 fReMakeIC(0),                     fMakeSeveralIC(0),
  52 fFillTMHisto(0),                  fFillSSHisto(0),
  53 // Several IC
  54 fNCones(0),                       fNPtThresFrac(0),
  55 fConeSizes(),                     fPtThresholds(),
  56 fPtFractions(),                   fSumPtThresholds(),
  57 // Histograms
  58 fhEIso(0),                        fhPtIso(0),
  59 fhPhiIso(0),                      fhEtaIso(0),                     fhEtaPhiIso(0),
  60 fhEtaPhiNoIso(0),
  61 fhPtNoIso(0),                     fhPtDecayIso(0),                 fhPtDecayNoIso(0),
  62 fhEtaPhiDecayIso(0),              fhEtaPhiDecayNoIso(0),
  63 fhConeSumPt(0),                   fhPtInCone(0),
  64 fhFRConeSumPt(0),                 fhPtInFRCone(0),
  65 // MC histograms
  66 fhPtIsoPrompt(0),                 fhPhiIsoPrompt(0),               fhEtaIsoPrompt(0),
  67 fhPtThresIsolatedPrompt(),        fhPtFracIsolatedPrompt(),        fhPtSumIsolatedPrompt(),
  68 fhPtIsoFragmentation(0),          fhPhiIsoFragmentation(0),        fhEtaIsoFragmentation(0),
  69 fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
  70 fhPtIsoPi0Decay(0),               fhPhiIsoPi0Decay(0),             fhEtaIsoPi0Decay(0),
  71 fhPtThresIsolatedPi0Decay(),      fhPtFracIsolatedPi0Decay(),      fhPtSumIsolatedPi0Decay(),
  72 fhPtIsoEtaDecay(0),               fhPhiIsoEtaDecay(0),             fhEtaIsoEtaDecay(0),
  73 fhPtThresIsolatedEtaDecay(),      fhPtFracIsolatedEtaDecay(),      fhPtSumIsolatedEtaDecay(),
  74 fhPtIsoOtherDecay(0),             fhPhiIsoOtherDecay(0),           fhEtaIsoOtherDecay(0),
  75 fhPtThresIsolatedOtherDecay(),    fhPtFracIsolatedOtherDecay(),    fhPtSumIsolatedOtherDecay(),
  76 fhPtIsoConversion(0),             fhPhiIsoConversion(0),           fhEtaIsoConversion(0),
  77 fhPtThresIsolatedConversion(),    fhPtFracIsolatedConversion(),    fhPtSumIsolatedConversion(),
  78 fhPtIsoUnknown(0),                fhPhiIsoUnknown(0),              fhEtaIsoUnknown(0),
  79 fhPtThresIsolatedUnknown(),       fhPtFracIsolatedUnknown(),       fhPtSumIsolatedUnknown(),
  80 fhPtNoIsoPi0Decay(0),             fhPtNoIsoEtaDecay(0),            fhPtNoIsoOtherDecay(0),
  81 fhPtNoIsoPrompt(0),               fhPtIsoMCPhoton(0),              fhPtNoIsoMCPhoton(0),
  82 fhPtNoIsoConversion(0),           fhPtNoIsoFragmentation(0),       fhPtNoIsoUnknown(0),
  83 // Hist several IC
  84 fhSumPtLeadingPt(),fhPtLeadingPt(),fhFRSumPtLeadingPt(),fhFRPtLeadingPt(),
  85 fhPtThresIsolated(),              fhPtFracIsolated(),              fhPtSumIsolated(),
  86 fhEtaPhiPtThresIso(),             fhEtaPhiPtThresDecayIso(),       fhPtPtThresDecayIso(),
  87 fhEtaPhiPtFracIso(),              fhEtaPhiPtFracDecayIso(),        fhPtPtFracDecayIso(),
  88 fhPtPtSumDecayIso(),              fhEtaPhiSumDensityIso(),         fhEtaPhiSumDensityDecayIso(),
  89 fhPtSumDensityIso(),              fhPtSumDensityDecayIso(),
  90 fhPtFracPtSumIso(),               fhPtFracPtSumDecayIso(),
  91 fhEtaPhiFracPtSumIso(),           fhEtaPhiFracPtSumDecayIso(),
  92 // Cluster control histograms
  93 fhTrackMatchedDEta(0x0),          fhTrackMatchedDPhi(0x0),         fhTrackMatchedDEtaDPhi(0x0),
  94 fhdEdx(0),                        fhEOverP(0),                     fhTrackMatchedMCParticle(0),
  95 fhELambda0(0),                    fhELambda1(0),
  96 fhELambda0TRD(0),                 fhELambda1TRD(0),
  97 // Number of local maxima in cluster
  98 fhNLocMax(0),
  99 fhELambda0LocMax1(0),             fhELambda1LocMax1(0),
 100 fhELambda0LocMax2(0),             fhELambda1LocMax2(0),
 101 fhELambda0LocMaxN(0),             fhELambda1LocMaxN(0),
 102 // Histograms settings
 103 fHistoNPtSumBins(0),              fHistoPtSumMax(0.),              fHistoPtSumMin(0.),
 104 fHistoNPtInConeBins(0),           fHistoPtInConeMax(0.),           fHistoPtInConeMin(0.)
 105 {
 106   //default ctor
 107
 108   //Initialize parameters
 109   InitParameters();
 110
 111   for(Int_t i = 0; i < 5 ; i++)
 112   {
 113     fConeSizes[i]      = 0 ;
 114
 115     fhPtSumIsolatedPrompt       [i] = 0 ;
 116     fhPtSumIsolatedFragmentation[i] = 0 ;
 117     fhPtSumIsolatedPi0Decay     [i] = 0 ;
 118     fhPtSumIsolatedEtaDecay     [i] = 0 ;
 119     fhPtSumIsolatedOtherDecay   [i] = 0 ;
 120     fhPtSumIsolatedConversion   [i] = 0 ;
 121     fhPtSumIsolatedUnknown      [i] = 0 ;
 122
 123     for(Int_t j = 0; j < 5 ; j++)
 124     {
 125       fhPtThresIsolated      [i][j] = 0 ;
 126       fhPtFracIsolated       [i][j] = 0 ;
 127       fhPtSumIsolated        [i][j] = 0 ;
 128
 129       fhEtaPhiPtThresIso     [i][j] = 0 ;
 130       fhEtaPhiPtThresDecayIso[i][j] = 0 ;
 131       fhPtPtThresDecayIso    [i][j] = 0 ;
 132
 133       fhEtaPhiPtFracIso      [i][j] = 0 ;
 134       fhEtaPhiPtFracDecayIso [i][j] = 0 ;
 135       fhPtPtFracDecayIso     [i][j] = 0 ;
 136       fhPtPtSumDecayIso      [i][j] = 0 ;
 137       fhPtSumDensityIso      [i][j] = 0 ;
 138       fhPtSumDensityDecayIso [i][j] = 0 ;
 139       fhEtaPhiSumDensityIso      [i][j] = 0 ;
 140       fhEtaPhiSumDensityDecayIso [i][j] = 0 ;
 141       fhPtFracPtSumIso       [i][j] = 0 ;
 142       fhPtFracPtSumDecayIso  [i][j] = 0 ;
 143       fhEtaPhiFracPtSumIso       [i][j] = 0 ;
 144       fhEtaPhiFracPtSumDecayIso  [i][j] = 0 ;
 145
 146       fhPtThresIsolatedPrompt       [i][j] = 0 ;
 147       fhPtThresIsolatedFragmentation[i][j] = 0 ;
 148       fhPtThresIsolatedPi0Decay     [i][j] = 0 ;
 149       fhPtThresIsolatedEtaDecay     [i][j] = 0 ;
 150       fhPtThresIsolatedOtherDecay   [i][j] = 0 ;
 151       fhPtThresIsolatedConversion   [i][j] = 0 ;
 152       fhPtThresIsolatedUnknown      [i][j] = 0 ;
 153
 154       fhPtFracIsolatedPrompt        [i][j] = 0 ;
 155       fhPtFracIsolatedFragmentation [i][j] = 0 ;
 156       fhPtFracIsolatedPi0Decay      [i][j] = 0 ;
 157       fhPtFracIsolatedEtaDecay      [i][j] = 0 ;
 158       fhPtFracIsolatedOtherDecay    [i][j] = 0 ;
 159       fhPtFracIsolatedConversion    [i][j] = 0 ;
 160       fhPtFracIsolatedUnknown       [i][j] = 0 ;
 161
 162     }
 163   }
 164
 165   for(Int_t i = 0; i < 5 ; i++)
 166   {
 167     fPtFractions    [i] = 0 ;
 168     fPtThresholds   [i] = 0 ;
 169     fSumPtThresholds[i] = 0 ;
 170   }
 171
 172 }
 173
 174 //________________________________________________________________________________________________
 175 void AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(
 176                                                                             const Int_t clusterID,
 177                                                                             const Int_t nMaxima,
 178                                                                             const Int_t mcTag
 179                                                                             )
 180 {
 181   // Fill Track matching and Shower Shape control histograms
 182
 183   if(!fFillTMHisto &&  !fFillSSHisto) return;
 184
 185   if(clusterID < 0 )
 186   {
 187     printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(), ID of cluster = %d, not possible! ", clusterID);
 188     return;
 189   }
 190
 191   Int_t iclus = -1;
 192   TObjArray* clusters = 0x0;
 193   if     (fCalorimeter == "EMCAL") clusters = GetEMCALClusters();
 194   else if(fCalorimeter == "PHOS" ) clusters = GetPHOSClusters();
 195
 196   if(clusters)
 197   {
 198
 199     AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
 200     Float_t energy = cluster->E();
 201
 202     if(fFillSSHisto)
 203     {
 204       fhELambda0   ->Fill(energy, cluster->GetM02() );
 205       fhELambda1   ->Fill(energy, cluster->GetM20() );
 206
 207       if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
 208       {
 209         fhELambda0TRD   ->Fill(energy, cluster->GetM02() );
 210         fhELambda1TRD   ->Fill(energy, cluster->GetM20() );
 211       }
 212
 213       fhNLocMax->Fill(energy,nMaxima);
 214       if     (nMaxima==1) { fhELambda0LocMax1->Fill(energy,cluster->GetM02()); fhELambda1LocMax1->Fill(energy,cluster->GetM20()); }
 215       else if(nMaxima==2) { fhELambda0LocMax2->Fill(energy,cluster->GetM02()); fhELambda1LocMax2->Fill(energy,cluster->GetM20()); }
 216       else                { fhELambda0LocMaxN->Fill(energy,cluster->GetM02()); fhELambda1LocMaxN->Fill(energy,cluster->GetM20()); }
 217
 218     } // SS histo fill
 219
 220
 221     if(fFillTMHisto)
 222     {
 223       Float_t dZ  = cluster->GetTrackDz();
 224       Float_t dR  = cluster->GetTrackDx();
 225
 226       if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
 227       {
 228         dR = 2000., dZ = 2000.;
 229         GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 230       }
 231
 232       //printf("ParticleIsolation: dPhi %f, dEta %f\n",dR,dZ);
 233       if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
 234       {
 235         fhTrackMatchedDEta->Fill(energy,dZ);
 236         fhTrackMatchedDPhi->Fill(energy,dR);
 237         if(energy > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
 238       }
 239
 240       // Check dEdx and E/p of matched clusters
 241
 242       if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
 243       {
 244
 245         AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
 246
 247         if(track)
 248         {
 249           Float_t dEdx = track->GetTPCsignal();
 250           fhdEdx->Fill(cluster->E(), dEdx);
 251
 252           Float_t eOverp = cluster->E()/track->P();
 253           fhEOverP->Fill(cluster->E(),  eOverp);
 254         }
 255         //else
 256         //  printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
 257
 258
 259         if(IsDataMC())
 260         {
 261           if ( !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)  )
 262           {
 263             if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
 264                       GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle->Fill(energy, 2.5 );
 265             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle->Fill(energy, 0.5 );
 266             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle->Fill(energy, 1.5 );
 267             else                                                                                   fhTrackMatchedMCParticle->Fill(energy, 3.5 );
 268
 269           }
 270           else
 271           {
 272             if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
 273                       GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle->Fill(energy, 6.5 );
 274             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle->Fill(energy, 4.5 );
 275             else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle->Fill(energy, 5.5 );
 276             else                                                                                   fhTrackMatchedMCParticle->Fill(energy, 7.5 );
 277           }
 278
 279         }  // MC
 280
 281       } // match window
 282
 283     }// TM histos fill
 284
 285   } // clusters array available
 286
 287 }
 288
 289 //______________________________________________________
 290 TObjString *  AliAnaParticleIsolation::GetAnalysisCuts()
 291 {
 292   //Save parameters used for analysis
 293   TString parList ; //this will be list of parameters used for this analysis.
 294   const Int_t buffersize = 255;
 295   char onePar[buffersize] ;
 296
 297   snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;
 298   parList+=onePar ;
 299   snprintf(onePar, buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 300   parList+=onePar ;
 301   snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
 302   parList+=onePar ;
 303   snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
 304   parList+=onePar ;
 305   snprintf(onePar, buffersize,"fFillTMHisto=%d (Flag for track matching histograms) \n",fFillTMHisto) ;
 306   parList+=onePar ;
 307   snprintf(onePar, buffersize,"fFillSSHisto=%d (Flag for shower shape histograms) \n",fFillSSHisto) ;
 308   parList+=onePar ;
 309
 310   if(fMakeSeveralIC)
 311   {
 312     snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
 313     parList+=onePar ;
 314     snprintf(onePar, buffersize,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
 315     parList+=onePar ;
 316
 317     for(Int_t icone = 0; icone < fNCones ; icone++)
 318     {
 319       snprintf(onePar, buffersize,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
 320       parList+=onePar ;
 321     }
 322     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
 323     {
 324       snprintf(onePar, buffersize,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
 325       parList+=onePar ;
 326     }
 327     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
 328     {
 329       snprintf(onePar, buffersize,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;
 330       parList+=onePar ;
 331     }
 332     for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
 333     {
 334       snprintf(onePar, buffersize,"fSumPtThresholds[%d]=%1.2f (isolation sum pt threshold) \n",ipt, fSumPtThresholds[ipt]) ;
 335       parList+=onePar ;
 336     }
 337   }
 338
 339   //Get parameters set in base class.
 340   parList += GetBaseParametersList() ;
 341
 342   //Get parameters set in IC class.
 343   if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;
 344
 345   return new TObjString(parList) ;
 346
 347 }
 348
 349 //________________________________________________________
 350 TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
 351 {
 352   // Create histograms to be saved in output file and
 353   // store them in outputContainer
 354   TList * outputContainer = new TList() ;
 355   outputContainer->SetName("IsolatedParticleHistos") ;
 356
 357   Int_t   nptbins  = GetHistogramRanges()->GetHistoPtBins();
 358   Int_t   nphibins = GetHistogramRanges()->GetHistoPhiBins();
 359   Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
 360   Float_t ptmax    = GetHistogramRanges()->GetHistoPtMax();
 361   Float_t phimax   = GetHistogramRanges()->GetHistoPhiMax();
 362   Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
 363   Float_t ptmin    = GetHistogramRanges()->GetHistoPtMin();
 364   Float_t phimin   = GetHistogramRanges()->GetHistoPhiMin();
 365   Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();
 366   Int_t   ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();
 367   Float_t ssmax    = GetHistogramRanges()->GetHistoShowerShapeMax();
 368   Float_t ssmin    = GetHistogramRanges()->GetHistoShowerShapeMin();
 369
 370   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
 371   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
 372   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
 373   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
 374   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
 375   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
 376
 377   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
 378   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
 379   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
 380   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
 381   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
 382   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
 383
 384
 385   Int_t   nptsumbins    = fHistoNPtSumBins;
 386   Float_t ptsummax      = fHistoPtSumMax;
 387   Float_t ptsummin      = fHistoPtSumMin;
 388   Int_t   nptinconebins = fHistoNPtInConeBins;
 389   Float_t ptinconemax   = fHistoPtInConeMax;
 390   Float_t ptinconemin   = fHistoPtInConeMin;
 391
 392   Float_t ptthre = GetIsolationCut()->GetPtThreshold();
 393   Float_t ptfrac = GetIsolationCut()->GetPtFraction();
 394   Float_t r      = GetIsolationCut()->GetConeSize();
 395
 396   if(!fMakeSeveralIC)
 397   {
 398     if(fFillTMHisto)
 399     {
 400       fhTrackMatchedDEta  = new TH2F
 401       ("hTrackMatchedDEta",
 402        Form("d#eta of cluster-track vs cluster energy for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 403        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
 404       fhTrackMatchedDEta->SetYTitle("d#eta");
 405       fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
 406
 407       fhTrackMatchedDPhi  = new TH2F
 408       ("hTrackMatchedDPhi",
 409        Form("d#phi of cluster-track vs cluster energy for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 410        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
 411       fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
 412       fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
 413
 414       fhTrackMatchedDEtaDPhi  = new TH2F
 415       ("hTrackMatchedDEtaDPhi",
 416        Form("d#eta vs d#phi of cluster-track for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 417        nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
 418       fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
 419       fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
 420
 421       outputContainer->Add(fhTrackMatchedDEta) ;
 422       outputContainer->Add(fhTrackMatchedDPhi) ;
 423       outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
 424
 425       fhdEdx  = new TH2F ("hdEdx",
 426                           Form("Matched track <dE/dx> vs cluster E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 427                           nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 428       fhdEdx->SetXTitle("E (GeV)");
 429       fhdEdx->SetYTitle("<dE/dx>");
 430       outputContainer->Add(fhdEdx);
 431
 432       fhEOverP  = new TH2F ("hEOverP",
 433                             Form("Matched track E/p vs cluster E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 434                             nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 435       fhEOverP->SetXTitle("E (GeV)");
 436       fhEOverP->SetYTitle("E/p");
 437       outputContainer->Add(fhEOverP);
 438
 439       if(IsDataMC())
 440       {
 441         fhTrackMatchedMCParticle  = new TH2F
 442         ("hTrackMatchedMCParticle",
 443          Form("Origin of particle vs energy vs cluster E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 444          nptbins,ptmin,ptmax,8,0,8);
 445         fhTrackMatchedMCParticle->SetXTitle("E (GeV)");
 446         //fhTrackMatchedMCParticle->SetYTitle("Particle type");
 447
 448         fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(1 ,"Photon");
 449         fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(2 ,"Electron");
 450         fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
 451         fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(4 ,"Rest");
 452         fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
 453         fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
 454         fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
 455         fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
 456
 457         outputContainer->Add(fhTrackMatchedMCParticle);
 458       }
 459     }
 460
 461     if(fFillSSHisto)
 462     {
 463       fhELambda0  = new TH2F
 464       ("hELambda0","Selected #pi^{0} pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 465       fhELambda0->SetYTitle("#lambda_{0}^{2}");
 466       fhELambda0->SetXTitle("E (GeV)");
 467       outputContainer->Add(fhELambda0) ;
 468
 469       fhELambda1  = new TH2F
 470       ("hELambda1","Selected #pi^{0} pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 471       fhELambda1->SetYTitle("#lambda_{1}^{2}");
 472       fhELambda1->SetXTitle("E (GeV)");
 473       outputContainer->Add(fhELambda1) ;
 474
 475       if(fCalorimeter=="EMCAL")
 476       {
 477         fhELambda0TRD  = new TH2F
 478         ("hELambda0TRD","Selected #pi^{0} pairs: E vs #lambda_{0}, SM behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 479         fhELambda0TRD->SetYTitle("#lambda_{0}^{2}");
 480         fhELambda0TRD->SetXTitle("E (GeV)");
 481         outputContainer->Add(fhELambda0TRD) ;
 482
 483         fhELambda1TRD  = new TH2F
 484         ("hELambda1TRD","Selected #pi^{0} pairs: E vs #lambda_{1}, SM behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 485         fhELambda1TRD->SetYTitle("#lambda_{1}^{2}");
 486         fhELambda1TRD->SetXTitle("E (GeV)");
 487         outputContainer->Add(fhELambda1TRD) ;
 488       }
 489
 490       fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
 491                            nptbins,ptmin,ptmax,10,0,10);
 492       fhNLocMax ->SetYTitle("N maxima");
 493       fhNLocMax ->SetXTitle("E (GeV)");
 494       outputContainer->Add(fhNLocMax) ;
 495
 496       fhELambda0LocMax1  = new TH2F
 497       ("hELambda0LocMax1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, 1 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 498       fhELambda0LocMax1->SetYTitle("#lambda_{0}^{2}");
 499       fhELambda0LocMax1->SetXTitle("E (GeV)");
 500       outputContainer->Add(fhELambda0LocMax1) ;
 501
 502       fhELambda1LocMax1  = new TH2F
 503       ("hELambda1LocMax1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, 1 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 504       fhELambda1LocMax1->SetYTitle("#lambda_{1}^{2}");
 505       fhELambda1LocMax1->SetXTitle("E (GeV)");
 506       outputContainer->Add(fhELambda1LocMax1) ;
 507
 508       fhELambda0LocMax2  = new TH2F
 509       ("hELambda0LocMax2","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, 2 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 510       fhELambda0LocMax2->SetYTitle("#lambda_{0}^{2}");
 511       fhELambda0LocMax2->SetXTitle("E (GeV)");
 512       outputContainer->Add(fhELambda0LocMax2) ;
 513
 514       fhELambda1LocMax2  = new TH2F
 515       ("hELambda1LocMax2","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, 2 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 516       fhELambda1LocMax2->SetYTitle("#lambda_{1}^{2}");
 517       fhELambda1LocMax2->SetXTitle("E (GeV)");
 518       outputContainer->Add(fhELambda1LocMax2) ;
 519
 520       fhELambda0LocMaxN  = new TH2F
 521       ("hELambda0LocMaxN","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, N>2 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 522       fhELambda0LocMaxN->SetYTitle("#lambda_{0}^{2}");
 523       fhELambda0LocMaxN->SetXTitle("E (GeV)");
 524       outputContainer->Add(fhELambda0LocMaxN) ;
 525
 526       fhELambda1LocMaxN  = new TH2F
 527       ("hELambda1LocMaxN","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, N>2 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 528       fhELambda1LocMaxN->SetYTitle("#lambda_{1}^{2}");
 529       fhELambda1LocMaxN->SetXTitle("E (GeV)");
 530       outputContainer->Add(fhELambda1LocMaxN) ;
 531
 532     }
 533
 534     fhConeSumPt  = new TH2F("hConePtSum",
 535                             Form("#Sigma p_{T} in isolation cone for R = %2.2f",r),
 536                             nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 537     fhConeSumPt->SetYTitle("#Sigma p_{T}");
 538     fhConeSumPt->SetXTitle("p_{T} (GeV/c)");
 539     outputContainer->Add(fhConeSumPt) ;
 540
 541     fhPtInCone  = new TH2F("hPtInCone",
 542                            Form("p_{T} in isolation cone for R = %2.2f",r),
 543                            nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 544     fhPtInCone->SetYTitle("p_{T in cone} (GeV/c)");
 545     fhPtInCone->SetXTitle("p_{T} (GeV/c)");
 546     outputContainer->Add(fhPtInCone) ;
 547
 548     fhFRConeSumPt  = new TH2F("hFRConePtSum",
 549                               Form("#Sigma p_{T} in the forward region isolation cone for R = %2.2f",r),
 550                               nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 551     fhFRConeSumPt->SetYTitle("#Sigma p_{T}");
 552     fhFRConeSumPt->SetXTitle("p_{T} (GeV/c)");
 553     outputContainer->Add(fhFRConeSumPt) ;
 554
 555     fhPtInFRCone  = new TH2F("hPtInFRCone",
 556                              Form("p_{T} in forward region isolation cone for R = %2.2f",r),
 557                              nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 558     fhPtInFRCone->SetYTitle("p_{T in cone} (GeV/c)");
 559     fhPtInFRCone->SetXTitle("p_{T} (GeV/c)");
 560     outputContainer->Add(fhPtInFRCone) ;
 561
 562     fhEIso   = new TH1F("hE",
 563                         Form("Number of isolated particles vs E for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 564                         nptbins,ptmin,ptmax);
 565     fhEIso->SetYTitle("dN / dE");
 566     fhEIso->SetXTitle("E (GeV/c)");
 567     outputContainer->Add(fhEIso) ;
 568
 569     fhPtIso  = new TH1F("hPt",
 570                         Form("Number of isolated particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 571                         nptbins,ptmin,ptmax);
 572     fhPtIso->SetYTitle("dN / p_{T}");
 573     fhPtIso->SetXTitle("p_{T} (GeV/c)");
 574     outputContainer->Add(fhPtIso) ;
 575
 576     fhPhiIso  = new TH2F("hPhi",
 577                          Form("Number of isolated particles vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 578                          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 579     fhPhiIso->SetYTitle("#phi");
 580     fhPhiIso->SetXTitle("p_{T} (GeV/c)");
 581     outputContainer->Add(fhPhiIso) ;
 582
 583     fhEtaIso  = new TH2F("hEta",
 584                          Form("Number of isolated particles vs #eta for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 585                          nptbins,ptmin,ptmax,netabins,etamin,etamax);
 586     fhEtaIso->SetYTitle("#eta");
 587     fhEtaIso->SetXTitle("p_{T} (GeV/c)");
 588     outputContainer->Add(fhEtaIso) ;
 589
 590     fhEtaPhiIso  = new TH2F("hEtaPhiIso",
 591                             Form("Number of isolated particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 592                             netabins,etamin,etamax,nphibins,phimin,phimax);
 593     fhEtaPhiIso->SetXTitle("#eta");
 594     fhEtaPhiIso->SetYTitle("#phi");
 595     outputContainer->Add(fhEtaPhiIso) ;
 596
 597     fhPtDecayIso  = new TH1F("hPtDecayIso",
 598                              Form("Number of isolated #pi^{0} decay particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 599                              nptbins,ptmin,ptmax);
 600     fhPtDecayIso->SetYTitle("N");
 601     fhPtDecayIso->SetXTitle("p_{T}(GeV/c)");
 602     outputContainer->Add(fhPtDecayIso) ;
 603
 604     fhEtaPhiDecayIso  = new TH2F("hEtaPhiDecayIso",
 605                                  Form("Number of isolated Pi0 decay particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 606                                  netabins,etamin,etamax,nphibins,phimin,phimax);
 607     fhEtaPhiDecayIso->SetXTitle("#eta");
 608     fhEtaPhiDecayIso->SetYTitle("#phi");
 609     outputContainer->Add(fhEtaPhiDecayIso) ;
 610
 611     if(IsDataMC())
 612     {
 613       fhPtIsoPrompt  = new TH1F("hPtMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax);
 614       fhPtIsoPrompt->SetYTitle("N");
 615       fhPtIsoPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
 616       outputContainer->Add(fhPtIsoPrompt) ;
 617
 618       fhPhiIsoPrompt  = new TH2F
 619       ("hPhiMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 620       fhPhiIsoPrompt->SetYTitle("#phi");
 621       fhPhiIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 622       outputContainer->Add(fhPhiIsoPrompt) ;
 623
 624       fhEtaIsoPrompt  = new TH2F
 625       ("hEtaMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 626       fhEtaIsoPrompt->SetYTitle("#eta");
 627       fhEtaIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 628       outputContainer->Add(fhEtaIsoPrompt) ;
 629
 630       fhPtIsoFragmentation  = new TH1F("hPtMCFragmentation","Number of isolated #gamma",nptbins,ptmin,ptmax);
 631       fhPtIsoFragmentation->SetYTitle("N");
 632       fhPtIsoFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
 633       outputContainer->Add(fhPtIsoFragmentation) ;
 634
 635       fhPhiIsoFragmentation  = new TH2F
 636       ("hPhiMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 637       fhPhiIsoFragmentation->SetYTitle("#phi");
 638       fhPhiIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 639       outputContainer->Add(fhPhiIsoFragmentation) ;
 640
 641       fhEtaIsoFragmentation  = new TH2F
 642       ("hEtaMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 643       fhEtaIsoFragmentation->SetYTitle("#eta");
 644       fhEtaIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 645       outputContainer->Add(fhEtaIsoFragmentation) ;
 646
 647       fhPtIsoPi0Decay  = new TH1F("hPtMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
 648       fhPtIsoPi0Decay->SetYTitle("N");
 649       fhPtIsoPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
 650       outputContainer->Add(fhPtIsoPi0Decay) ;
 651
 652       fhPhiIsoPi0Decay  = new TH2F
 653       ("hPhiMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 654       fhPhiIsoPi0Decay->SetYTitle("#phi");
 655       fhPhiIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 656       outputContainer->Add(fhPhiIsoPi0Decay) ;
 657
 658       fhEtaIsoPi0Decay  = new TH2F
 659       ("hEtaMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 660       fhEtaIsoPi0Decay->SetYTitle("#eta");
 661       fhEtaIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 662       outputContainer->Add(fhEtaIsoPi0Decay) ;
 663
 664       fhPtIsoEtaDecay  = new TH1F("hPtMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax);
 665       fhPtIsoEtaDecay->SetYTitle("N");
 666       fhPtIsoEtaDecay->SetXTitle("p_{T #gamma}(GeV/c)");
 667       outputContainer->Add(fhPtIsoEtaDecay) ;
 668
 669       fhPhiIsoEtaDecay  = new TH2F
 670       ("hPhiMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 671       fhPhiIsoEtaDecay->SetYTitle("#phi");
 672       fhPhiIsoEtaDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 673       outputContainer->Add(fhPhiIsoEtaDecay) ;
 674
 675       fhEtaIsoEtaDecay  = new TH2F
 676       ("hEtaMCEtaDecay","Number of isolated #gamma from #eta decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 677       fhEtaIsoEtaDecay->SetYTitle("#eta");
 678       fhEtaIsoEtaDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 679       outputContainer->Add(fhEtaIsoEtaDecay) ;
 680
 681       fhPtIsoOtherDecay  = new TH1F("hPtMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax);
 682       fhPtIsoOtherDecay->SetYTitle("N");
 683       fhPtIsoOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
 684       outputContainer->Add(fhPtIsoOtherDecay) ;
 685
 686       fhPhiIsoOtherDecay  = new TH2F
 687       ("hPhiMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 688       fhPhiIsoOtherDecay->SetYTitle("#phi");
 689       fhPhiIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 690       outputContainer->Add(fhPhiIsoOtherDecay) ;
 691
 692       fhEtaIsoOtherDecay  = new TH2F
 693       ("hEtaMCOtherDecay","Number of isolated #gamma non-#pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 694       fhEtaIsoOtherDecay->SetYTitle("#eta");
 695       fhEtaIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 696       outputContainer->Add(fhEtaIsoOtherDecay) ;
 697
 698       fhPtIsoConversion  = new TH1F("hPtMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax);
 699       fhPtIsoConversion->SetYTitle("N");
 700       fhPtIsoConversion->SetXTitle("p_{T #gamma}(GeV/c)");
 701       outputContainer->Add(fhPtIsoConversion) ;
 702
 703       fhPhiIsoConversion  = new TH2F
 704       ("hPhiMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 705       fhPhiIsoConversion->SetYTitle("#phi");
 706       fhPhiIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 707       outputContainer->Add(fhPhiIsoConversion) ;
 708
 709       fhEtaIsoConversion  = new TH2F
 710       ("hEtaMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 711       fhEtaIsoConversion->SetYTitle("#eta");
 712       fhEtaIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 713       outputContainer->Add(fhEtaIsoConversion) ;
 714
 715       fhPtIsoUnknown  = new TH1F("hPtMCUnknown","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax);
 716       fhPtIsoUnknown->SetYTitle("N");
 717       fhPtIsoUnknown->SetXTitle("p_{T}(GeV/c)");
 718       outputContainer->Add(fhPtIsoUnknown) ;
 719
 720       fhPhiIsoUnknown  = new TH2F
 721       ("hPhiMCUnknown","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 722       fhPhiIsoUnknown->SetYTitle("#phi");
 723       fhPhiIsoUnknown->SetXTitle("p_{T} (GeV/c)");
 724       outputContainer->Add(fhPhiIsoUnknown) ;
 725
 726       fhEtaIsoUnknown  = new TH2F
 727       ("hEtaMCUnknown","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 728       fhEtaIsoUnknown->SetYTitle("#eta");
 729       fhEtaIsoUnknown->SetXTitle("p_{T} (GeV/c)");
 730       outputContainer->Add(fhEtaIsoUnknown) ;
 731
 732     }//Histos with MC
 733
 734   }
 735
 736   // Not Isolated histograms, reference histograms
 737
 738   fhPtNoIso  = new TH1F("hPtNoIso",
 739                         Form("Number of not isolated leading particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 740                         nptbins,ptmin,ptmax);
 741   fhPtNoIso->SetYTitle("N");
 742   fhPtNoIso->SetXTitle("p_{T}(GeV/c)");
 743   outputContainer->Add(fhPtNoIso) ;
 744
 745
 746   fhEtaPhiNoIso  = new TH2F("hEtaPhiNoIso",
 747                             Form("Number of not isolated leading particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 748                             netabins,etamin,etamax,nphibins,phimin,phimax);
 749   fhEtaPhiNoIso->SetXTitle("#eta");
 750   fhEtaPhiNoIso->SetYTitle("#phi");
 751   outputContainer->Add(fhEtaPhiNoIso) ;
 752
 753   fhPtDecayNoIso  = new TH1F("hPtDecayNoIso",
 754                              Form("Number of not isolated leading pi0 decay particles vs p_{T} for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 755                              nptbins,ptmin,ptmax);
 756   fhPtDecayNoIso->SetYTitle("N");
 757   fhPtDecayNoIso->SetXTitle("p_{T}(GeV/c)");
 758   outputContainer->Add(fhPtDecayNoIso) ;
 759
 760   fhEtaPhiDecayNoIso  = new TH2F("hEtaPhiDecayNoIso",
 761                                  Form("Number of not isolated leading Pi0 decay particles #eta vs #phi for R = %2.2f, p_{T}^{th} = %2.2f, p_{T}^{fr} = %2.2f",r,ptthre,ptfrac),
 762                                  netabins,etamin,etamax,nphibins,phimin,phimax);
 763   fhEtaPhiDecayNoIso->SetXTitle("#eta");
 764   fhEtaPhiDecayNoIso->SetYTitle("#phi");
 765   outputContainer->Add(fhEtaPhiDecayNoIso) ;
 766
 767
 768
 769   if(IsDataMC())
 770   {
 771     fhPtNoIsoPi0Decay  = new TH1F
 772     ("hPtNoIsoPi0Decay","Number of not isolated leading #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
 773     fhPtNoIsoPi0Decay->SetYTitle("N");
 774     fhPtNoIsoPi0Decay->SetXTitle("p_{T} (GeV/c)");
 775     outputContainer->Add(fhPtNoIsoPi0Decay) ;
 776
 777     fhPtNoIsoEtaDecay  = new TH1F
 778     ("hPtNoIsoEtaDecay","Number of not isolated leading #gamma from eta decay",nptbins,ptmin,ptmax);
 779     fhPtNoIsoEtaDecay->SetYTitle("N");
 780     fhPtNoIsoEtaDecay->SetXTitle("p_{T} (GeV/c)");
 781     outputContainer->Add(fhPtNoIsoEtaDecay) ;
 782
 783     fhPtNoIsoOtherDecay  = new TH1F
 784     ("hPtNoIsoOtherDecay","Number of not isolated leading #gamma from other decay",nptbins,ptmin,ptmax);
 785     fhPtNoIsoOtherDecay->SetYTitle("N");
 786     fhPtNoIsoOtherDecay->SetXTitle("p_{T} (GeV/c)");
 787     outputContainer->Add(fhPtNoIsoOtherDecay) ;
 788
 789     fhPtNoIsoPrompt  = new TH1F
 790     ("hPtNoIsoPrompt","Number of not isolated leading prompt #gamma",nptbins,ptmin,ptmax);
 791     fhPtNoIsoPrompt->SetYTitle("N");
 792     fhPtNoIsoPrompt->SetXTitle("p_{T} (GeV/c)");
 793     outputContainer->Add(fhPtNoIsoPrompt) ;
 794
 795     fhPtIsoMCPhoton  = new TH1F
 796     ("hPtIsoMCPhoton","Number of isolated leading  #gamma",nptbins,ptmin,ptmax);
 797     fhPtIsoMCPhoton->SetYTitle("N");
 798     fhPtIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
 799     outputContainer->Add(fhPtIsoMCPhoton) ;
 800
 801     fhPtNoIsoMCPhoton  = new TH1F
 802     ("hPtNoIsoMCPhoton","Number of not isolated leading #gamma",nptbins,ptmin,ptmax);
 803     fhPtNoIsoMCPhoton->SetYTitle("N");
 804     fhPtNoIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
 805     outputContainer->Add(fhPtNoIsoMCPhoton) ;
 806
 807     fhPtNoIsoConversion  = new TH1F
 808     ("hPtNoIsoConversion","Number of not isolated leading conversion #gamma",nptbins,ptmin,ptmax);
 809     fhPtNoIsoConversion->SetYTitle("N");
 810     fhPtNoIsoConversion->SetXTitle("p_{T} (GeV/c)");
 811     outputContainer->Add(fhPtNoIsoConversion) ;
 812
 813     fhPtNoIsoFragmentation  = new TH1F
 814     ("hPtNoIsoFragmentation","Number of not isolated leading fragmentation #gamma",nptbins,ptmin,ptmax);
 815     fhPtNoIsoFragmentation->SetYTitle("N");
 816     fhPtNoIsoFragmentation->SetXTitle("p_{T} (GeV/c)");
 817     outputContainer->Add(fhPtNoIsoFragmentation) ;
 818
 819     fhPtNoIsoUnknown  = new TH1F
 820     ("hPtNoIsoUnknown","Number of not isolated leading hadrons",nptbins,ptmin,ptmax);
 821     fhPtNoIsoUnknown->SetYTitle("N");
 822     fhPtNoIsoUnknown->SetXTitle("p_{T} (GeV/c)");
 823     outputContainer->Add(fhPtNoIsoUnknown) ;
 824
 825   }//Histos with MC
 826
 827
 828   if(fMakeSeveralIC)
 829   {
 830     const Int_t buffersize = 255;
 831     char name[buffersize];
 832     char title[buffersize];
 833     for(Int_t icone = 0; icone<fNCones; icone++)
 834     {
 835       // sum pt in cone vs. pt leading
 836       snprintf(name, buffersize,"hSumPtLeadingPt_Cone_%d",icone);
 837       snprintf(title, buffersize,"#Sigma p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
 838       fhSumPtLeadingPt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 839       fhSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}p_{T} (GeV/c)");//#Sigma p_{T}
 840       fhSumPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
 841       outputContainer->Add(fhSumPtLeadingPt[icone]) ;
 842
 843       // pt in cone vs. pt leading
 844       snprintf(name, buffersize,"hPtLeadingPt_Cone_%d",icone);
 845       snprintf(title, buffersize,"p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
 846       fhPtLeadingPt[icone]  = new TH2F(name, title,  nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 847       fhPtLeadingPt[icone] ->SetYTitle("p_{T}^{cone} (GeV/c)");
 848       fhPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
 849       outputContainer->Add(fhPtLeadingPt[icone]) ;
 850
 851        // sum pt in cone vs. pt leading in the forward region (for background subtraction studies)
 852         snprintf(name, buffersize,"hFRSumPtLeadingPt_Cone_%d",icone);
 853       snprintf(title, buffersize,"#Sigma p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
 854       fhFRSumPtLeadingPt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 855       fhFRSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}p_{T} (GeV/c)");//#Sigma p_{T}
 856       fhFRSumPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
 857       outputContainer->Add(fhFRSumPtLeadingPt[icone]) ;
 858
 859       // pt in cone vs. pt leading in the forward region (for background subtraction studies)
 860         snprintf(name, buffersize,"hFRPtLeadingPt_Cone_%d",icone);
 861       snprintf(title, buffersize,"p_{T} in isolation cone for R = %2.2f",fConeSizes[icone]);
 862       fhFRPtLeadingPt[icone]  = new TH2F(name, title,  nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
 863       fhFRPtLeadingPt[icone] ->SetYTitle("p_{T}^{cone} (GeV/c)");
 864       fhFRPtLeadingPt[icone] ->SetXTitle("p_{T}^{leading} (GeV/c)");
 865       outputContainer->Add(fhFRPtLeadingPt[icone]) ;
 866
 867
 868       if(IsDataMC())
 869       {
 870         snprintf(name, buffersize,"hPtSumPrompt_Cone_%d",icone);
 871         snprintf(title, buffersize,"Candidate Prompt cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
 872         fhPtSumIsolatedPrompt[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 873         fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 874         fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
 875         outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ;
 876
 877         snprintf(name, buffersize,"hPtSumFragmentation_Cone_%d",icone);
 878         snprintf(title, buffersize,"Candidate Fragmentation cone sum p_{T} for R = %2.2fvs candidate p_{T}",fConeSizes[icone]);
 879         fhPtSumIsolatedFragmentation[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 880         fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 881         fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
 882         outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ;
 883
 884         snprintf(name, buffersize,"hPtSumPi0Decay_Cone_%d",icone);
 885         snprintf(title, buffersize,"Candidate Pi0Decay cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
 886         fhPtSumIsolatedPi0Decay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 887         fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 888         fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
 889         outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ;
 890
 891         snprintf(name, buffersize,"hPtSumEtaDecay_Cone_%d",icone);
 892         snprintf(title, buffersize,"Candidate EtaDecay cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
 893         fhPtSumIsolatedEtaDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 894         fhPtSumIsolatedEtaDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 895         fhPtSumIsolatedEtaDecay[icone]->SetXTitle("p_{T} (GeV/c)");
 896         outputContainer->Add(fhPtSumIsolatedEtaDecay[icone]) ;
 897
 898         snprintf(name, buffersize,"hPtSumOtherDecay_Cone_%d",icone);
 899         snprintf(title, buffersize,"Candidate OtherDecay cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
 900         fhPtSumIsolatedOtherDecay[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 901         fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 902         fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
 903         outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ;
 904
 905         snprintf(name, buffersize,"hPtSumConversion_Cone_%d",icone);
 906         snprintf(title, buffersize,"Candidate Conversion cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
 907         fhPtSumIsolatedConversion[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 908         fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 909         fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
 910         outputContainer->Add(fhPtSumIsolatedConversion[icone]) ;
 911
 912         snprintf(name, buffersize,"hPtSumUnknown_Cone_%d",icone);
 913         snprintf(title, buffersize,"Candidate Unknown cone sum p_{T} for R = %2.2f vs candidate p_{T}",fConeSizes[icone]);
 914         fhPtSumIsolatedUnknown[icone]  = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 915         fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 916         fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");
 917         outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ;
 918
 919       }//Histos with MC
 920
 921       for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++)
 922       {
 923
 924         snprintf(name, buffersize,"hPtThres_Cone_%d_Pt%d",icone,ipt);
 925         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
 926         fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 927         fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 928         outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
 929
 930         snprintf(name, buffersize,"hPtFrac_Cone_%d_Pt%d",icone,ipt);
 931         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
 932         fhPtFracIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 933         fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 934         outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
 935
 936
 937         snprintf(name, buffersize,"hPtSum_Cone_%d_Pt%d",icone,ipt);
 938         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
 939         fhPtSumIsolated[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 940         // fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma p_{T} (GeV/c)");
 941         fhPtSumIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 942         outputContainer->Add(fhPtSumIsolated[icone][ipt]) ;
 943
 944         snprintf(name, buffersize,"hPtSumDensity_Cone_%d_Pt%d",icone,ipt);
 945         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for density in R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
 946         fhPtSumDensityIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
 947         //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma p_{T} (GeV/c)");
 948         fhPtSumDensityIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 949         outputContainer->Add(fhPtSumDensityIso[icone][ipt]) ;
 950
 951         snprintf(name, buffersize,"hPtFracPtSum_Cone_%d_Pt%d",icone,ipt);
 952         snprintf(title, buffersize,"Isolated candidate p_{T} distribution for PtFracPtSum in R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
 953         fhPtFracPtSumIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
 954         //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma p_{T} (GeV/c)");
 955         fhPtFracPtSumIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 956         outputContainer->Add(fhPtFracPtSumIso[icone][ipt]) ;
 957
 958         // pt decays isolated
 959         snprintf(name, buffersize,"hPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
 960         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
 961         fhPtPtThresDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 962         fhPtPtThresDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 963         outputContainer->Add(fhPtPtThresDecayIso[icone][ipt]) ;
 964
 965         snprintf(name, buffersize,"hPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
 966         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
 967         fhPtPtFracDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
 968         fhPtPtFracDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 969         outputContainer->Add(fhPtPtFracDecayIso[icone][ipt]) ;
 970
 971         snprintf(name, buffersize,"hPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
 972         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
 973         fhPtPtSumDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
 974         //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 975         fhPtPtSumDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 976         outputContainer->Add(fhPtPtSumDecayIso[icone][ipt]) ;
 977
 978         snprintf(name, buffersize,"hPtSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
 979         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for density in R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
 980         fhPtSumDensityDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
 981         //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 982         fhPtSumDensityDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 983         outputContainer->Add(fhPtSumDensityDecayIso[icone][ipt]) ;
 984
 985         snprintf(name, buffersize,"hPtFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
 986         snprintf(title, buffersize,"Isolated decay candidate p_{T} distribution for PtFracPtSum in R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
 987         fhPtFracPtSumDecayIso[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
 988         //  fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 989         fhPtFracPtSumDecayIso[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 990         outputContainer->Add(fhPtFracPtSumDecayIso[icone][ipt]) ;
 991
 992
 993         // eta:phi
 994         snprintf(name, buffersize,"hEtaPhiPtThres_Cone_%d_Pt%d",icone,ipt);
 995         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
 996         fhEtaPhiPtThresIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
 997         fhEtaPhiPtThresIso[icone][ipt]->SetXTitle("#eta");
 998         fhEtaPhiPtThresIso[icone][ipt]->SetYTitle("#phi");
 999         outputContainer->Add(fhEtaPhiPtThresIso[icone][ipt]) ;
1000
1001         snprintf(name, buffersize,"hEtaPhiPtFrac_Cone_%d_Pt%d",icone,ipt);
1002         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1003         fhEtaPhiPtFracIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1004         fhEtaPhiPtFracIso[icone][ipt]->SetXTitle("#eta");
1005         fhEtaPhiPtFracIso[icone][ipt]->SetYTitle("#phi");
1006         outputContainer->Add(fhEtaPhiPtFracIso[icone][ipt]) ;
1007
1008         snprintf(name, buffersize,"hEtaPhiPtSum_Cone_%d_Pt%d",icone,ipt);
1009         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1010         fhEtaPhiPtSumIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1011         fhEtaPhiPtSumIso[icone][ipt]->SetXTitle("#eta");
1012         fhEtaPhiPtSumIso[icone][ipt]->SetYTitle("#phi");
1013         outputContainer->Add(fhEtaPhiPtSumIso[icone][ipt]) ;
1014
1015         snprintf(name, buffersize,"hEtaPhiSumDensity_Cone_%d_Pt%d",icone,ipt);
1016         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for density R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1017         fhEtaPhiSumDensityIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1018         fhEtaPhiSumDensityIso[icone][ipt]->SetXTitle("#eta");
1019         fhEtaPhiSumDensityIso[icone][ipt]->SetYTitle("#phi");
1020         outputContainer->Add(fhEtaPhiSumDensityIso[icone][ipt]) ;
1021
1022         snprintf(name, buffersize,"hEtaPhiFracPtSum_Cone_%d_Pt%d",icone,ipt);
1023         snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for FracPtSum R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1024         fhEtaPhiFracPtSumIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1025         fhEtaPhiFracPtSumIso[icone][ipt]->SetXTitle("#eta");
1026         fhEtaPhiFracPtSumIso[icone][ipt]->SetYTitle("#phi");
1027         outputContainer->Add(fhEtaPhiFracPtSumIso[icone][ipt]) ;
1028
1029         // eta:phi decays
1030         snprintf(name, buffersize,"hEtaPhiPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
1031         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for R = %2.2f and p_{T}^{th} = %2.2f GeV/c",fConeSizes[icone],fPtThresholds[ipt]);
1032         fhEtaPhiPtThresDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1033         fhEtaPhiPtThresDecayIso[icone][ipt]->SetXTitle("#eta");
1034         fhEtaPhiPtThresDecayIso[icone][ipt]->SetYTitle("#phi");
1035         outputContainer->Add(fhEtaPhiPtThresDecayIso[icone][ipt]) ;
1036
1037         snprintf(name, buffersize,"hEtaPhiPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
1038         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1039         fhEtaPhiPtFracDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1040         fhEtaPhiPtFracDecayIso[icone][ipt]->SetXTitle("#eta");
1041         fhEtaPhiPtFracDecayIso[icone][ipt]->SetYTitle("#phi");
1042         outputContainer->Add(fhEtaPhiPtFracDecayIso[icone][ipt]) ;
1043
1044
1045         snprintf(name, buffersize,"hEtaPhiPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
1046         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1047         fhEtaPhiPtSumDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1048         fhEtaPhiPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
1049         fhEtaPhiPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
1050         outputContainer->Add(fhEtaPhiPtSumDecayIso[icone][ipt]) ;
1051
1052         snprintf(name, buffersize,"hEtaPhiSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
1053         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for density R = %2.2f and p_{T}^{sum} = %2.2f GeV/c",fConeSizes[icone],fSumPtThresholds[ipt]);
1054         fhEtaPhiSumDensityDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1055         fhEtaPhiSumDensityDecayIso[icone][ipt]->SetXTitle("#eta");
1056         fhEtaPhiSumDensityDecayIso[icone][ipt]->SetYTitle("#phi");
1057         outputContainer->Add(fhEtaPhiSumDensityDecayIso[icone][ipt]) ;
1058
1059         snprintf(name, buffersize,"hEtaPhiFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
1060         snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for FracPtSum R = %2.2f and p_{T}^{fr} = %2.2f GeV/c",fConeSizes[icone],fPtFractions[ipt]);
1061         fhEtaPhiFracPtSumDecayIso[icone][ipt]  = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
1062         fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
1063         fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
1064         outputContainer->Add(fhEtaPhiFracPtSumDecayIso[icone][ipt]) ;
1065
1066
1067         if(IsDataMC())
1068         {
1069           snprintf(name, buffersize,"hPtThresMCPrompt_Cone_%d_Pt%d",icone,ipt);
1070           snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1071           fhPtThresIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1072           fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1073           outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ;
1074
1075           snprintf(name, buffersize,"hPtFracMCPrompt_Cone_%d_Pt%d",icone,ipt);
1076           snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1077           fhPtFracIsolatedPrompt[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1078           fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1079           outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ;
1080
1081           snprintf(name, buffersize,"hPtThresMCFragmentation_Cone_%d_Pt%d",icone,ipt);
1082           snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1083           fhPtThresIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1084           fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1085           outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ;
1086
1087           snprintf(name, buffersize,"hPtFracMCFragmentation_Cone_%d_Pt%d",icone,ipt);
1088           snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1089           fhPtFracIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1090           fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1091           outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ;
1092
1093           snprintf(name, buffersize,"hPtThresMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
1094           snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1095           fhPtThresIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1096           fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1097           outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ;
1098
1099           snprintf(name, buffersize,"hPtFracMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
1100           snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1101           fhPtFracIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1102           fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1103           outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ;
1104
1105           snprintf(name, buffersize,"hPtThresMCEtaDecay_Cone_%d_Pt%d",icone,ipt);
1106           snprintf(title, buffersize,"Isolated candidate EtaDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1107           fhPtThresIsolatedEtaDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1108           fhPtThresIsolatedEtaDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1109           outputContainer->Add(fhPtThresIsolatedEtaDecay[icone][ipt]) ;
1110
1111           snprintf(name, buffersize,"hPtFracMCEtaDecay_Cone_%d_Pt%d",icone,ipt);
1112           snprintf(title, buffersize,"Isolated candidate EtaDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1113           fhPtFracIsolatedEtaDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1114           fhPtFracIsolatedEtaDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1115           outputContainer->Add(fhPtFracIsolatedEtaDecay[icone][ipt]) ;
1116
1117
1118           snprintf(name, buffersize,"hPtThresMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
1119           snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1120           fhPtThresIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1121           fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1122           outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ;
1123
1124           snprintf(name, buffersize,"hPtFracMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
1125           snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1126           fhPtFracIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1127           fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1128           outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
1129
1130           snprintf(name, buffersize,"hPtThresMCConversion_Cone_%d_Pt%d",icone,ipt);
1131           snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1132           fhPtThresIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1133           fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1134           outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ;
1135
1136           snprintf(name, buffersize,"hPtFracMCConversion_Cone_%d_Pt%d",icone,ipt);
1137           snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1138           fhPtFracIsolatedConversion[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1139           fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1140           outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
1141
1142           snprintf(name, buffersize,"hPtThresMCUnknown_Cone_%d_Pt%d",icone,ipt);
1143           snprintf(title, buffersize,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1144           fhPtThresIsolatedUnknown[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1145           fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1146           outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ;
1147
1148           snprintf(name, buffersize,"hPtFracMCUnknown_Cone_%d_Pt%d",icone,ipt);
1149           snprintf(title, buffersize,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
1150           fhPtFracIsolatedUnknown[icone][ipt]  = new TH1F(name, title,nptbins,ptmin,ptmax);
1151           fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
1152           outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;
1153
1154         }//Histos with MC
1155       }//icone loop
1156     }//ipt loop
1157   }
1158
1159   return outputContainer ;
1160
1161 }
1162
1163 //__________________________________
1164 void AliAnaParticleIsolation::Init()
1165 {
1166   // Do some checks and init stuff
1167
1168   // In case of several cone and thresholds analysis, open the cuts for the filling of the
1169   // track and cluster reference arrays in cone when done in the MakeAnalysisFillAOD().
1170   // The different cones, thresholds are tested for this list of tracks, clusters.
1171   if(fMakeSeveralIC)
1172   {
1173     printf("AliAnaParticleIsolation::Init() - Open default isolation cuts for multiple Isolation analysis\n");
1174     GetIsolationCut()->SetPtThreshold(100);
1175     GetIsolationCut()->SetPtFraction(100);
1176     GetIsolationCut()->SetConeSize(1);
1177   }
1178 }
1179
1180 //____________________________________________
1181 void AliAnaParticleIsolation::InitParameters()
1182 {
1183
1184   //Initialize the parameters of the analysis.
1185   SetInputAODName("PWG4Particle");
1186   SetAODObjArrayName("IsolationCone");
1187   AddToHistogramsName("AnaIsolation_");
1188
1189   fCalorimeter = "PHOS" ;
1190   fReMakeIC = kFALSE ;
1191   fMakeSeveralIC = kFALSE ;
1192
1193   //----------- Several IC-----------------
1194   fNCones             = 5 ;
1195   fNPtThresFrac       = 5 ;
1196   fConeSizes      [0] = 0.1;     fConeSizes      [1] = 0.2;   fConeSizes      [2] = 0.3; fConeSizes      [3] = 0.4;  fConeSizes      [4] = 0.5;
1197   fPtThresholds   [0] = 1.;      fPtThresholds   [1] = 2.;    fPtThresholds   [2] = 3.;  fPtThresholds   [3] = 4.;   fPtThresholds   [4] = 5.;
1198   fPtFractions    [0] = 0.05;    fPtFractions    [1] = 0.075; fPtFractions    [2] = 0.1; fPtFractions    [3] = 1.25; fPtFractions    [4] = 1.5;
1199   fSumPtThresholds[0] = 1.;      fSumPtThresholds[1] = 2.;    fSumPtThresholds[2] = 3.;  fSumPtThresholds[3] = 4.;   fSumPtThresholds[4] = 5.;
1200
1201   //------------- Histograms settings -------
1202   fHistoNPtSumBins = 100 ;
1203   fHistoPtSumMax   = 50 ;
1204   fHistoPtSumMin   = 0.  ;
1205
1206   fHistoNPtInConeBins = 100 ;
1207   fHistoPtInConeMax   = 50 ;
1208   fHistoPtInConeMin   = 0.  ;
1209
1210 }
1211
1212 //__________________________________________________
1213 void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
1214 {
1215   //Do analysis and fill aods
1216   //Search for the isolated photon in fCalorimeter with pt > GetMinPt()
1217
1218   if(!GetInputAODBranch())
1219   {
1220     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1221     abort();
1222   }
1223
1224   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1225   {
1226     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
1227     abort();
1228   }
1229
1230   Int_t n = 0, nfrac = 0;
1231   Bool_t  isolated  = kFALSE ;
1232   Float_t coneptsum = 0 ;
1233   TObjArray * pl    = 0x0; ;
1234
1235   //Select the calorimeter for candidate isolation with neutral particles
1236   if      (fCalorimeter == "PHOS" )
1237     pl = GetPHOSClusters();
1238   else if (fCalorimeter == "EMCAL")
1239     pl = GetEMCALClusters();
1240
1241   //Loop on AOD branch, filled previously in AliAnaPhoton, find leading particle to do isolation only with it
1242   Double_t ptLeading = 0. ;
1243   Int_t    idLeading = -1 ;
1244   TLorentzVector mom ;
1245   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1246   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Input aod branch entries %d\n", naod);
1247
1248   for(Int_t iaod = 0; iaod < naod; iaod++)
1249   {
1250     AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1251
1252     //If too small or too large pt, skip
1253     if(aodinput->Pt() < GetMinPt() || aodinput->Pt() > GetMaxPt() ) continue ;
1254
1255     //check if it is low pt trigger particle
1256     if((aodinput->Pt() < GetIsolationCut()->GetPtThreshold() ||
1257         aodinput->Pt() < GetIsolationCut()->GetSumPtThreshold()) &&
1258        !fMakeSeveralIC)
1259     {
1260       continue ; //trigger should not come from underlying event
1261     }
1262
1263     //vertex cut in case of mixing
1264     Int_t check = CheckMixedEventVertex(aodinput->GetCaloLabel(0), aodinput->GetTrackLabel(0));
1265     if(check ==  0) continue;
1266     if(check == -1) return;
1267
1268     //find the leading particles with highest momentum
1269     if ( aodinput->Pt() > ptLeading )
1270     {
1271       ptLeading = aodinput->Pt() ;
1272       idLeading = iaod ;
1273     }
1274
1275     aodinput->SetLeadingParticle(kFALSE);
1276
1277   }//finish searching for leading trigger particle
1278
1279   // Check isolation of leading particle
1280   if(idLeading < 0) return;
1281
1282   AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(idLeading));
1283   aodinput->SetLeadingParticle(kTRUE);
1284
1285   // Check isolation only of clusters in fiducial region
1286   if(IsFiducialCutOn())
1287   {
1288     Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodinput->Momentum(),aodinput->GetDetector()) ;
1289     if(! in ) return ;
1290   }
1291
1292   //After cuts, study isolation
1293   n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
1294   GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,
1295                                       GetReader(), GetCaloPID(),
1296                                       kTRUE, aodinput, GetAODObjArrayName(),
1297                                       n,nfrac,coneptsum, isolated);
1298
1299   if(!fMakeSeveralIC) aodinput->SetIsolated(isolated);
1300
1301   if(GetDebug() > 1)
1302   {
1303     if(isolated)printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",idLeading);
1304     printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");
1305   }
1306
1307 }
1308
1309 //_________________________________________________________
1310 void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
1311 {
1312   //Do analysis and fill histograms
1313
1314   Int_t   n = 0, nfrac = 0;
1315   Bool_t  isolated  = kFALSE ;
1316   Float_t coneptsum = 0 ;
1317
1318   //Loop on stored AOD
1319   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1320   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
1321
1322   //Get vertex for photon momentum calculation
1323   Double_t vertex[] = {0,0,0} ; //vertex ;
1324   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1325   {
1326     GetReader()->GetVertex(vertex);
1327   }
1328
1329   for(Int_t iaod = 0; iaod < naod ; iaod++)
1330   {
1331     AliAODPWG4ParticleCorrelation* aod =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1332
1333     if(!aod->IsLeadingParticle()) continue; // Try to isolate only leading cluster or track
1334
1335     // Check isolation only of clusters in fiducial region
1336     if(IsFiducialCutOn())
1337     {
1338       Bool_t in = GetFiducialCut()->IsInFiducialCut(*aod->Momentum(),aod->GetDetector()) ;
1339       if(! in ) continue ;
1340     }
1341
1342     Bool_t  isolation  = aod->IsIsolated();
1343     Bool_t  decay      = aod->IsTagged();
1344     Float_t energy     = aod->E();
1345     Float_t pt         = aod->Pt();
1346     Float_t phi        = aod->Phi();
1347     Float_t eta        = aod->Eta();
1348     Float_t conesize   = GetIsolationCut()->GetConeSize();
1349
1350     //Recover reference arrays with clusters and tracks
1351     TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
1352     TObjArray * reftracks   = aod->GetObjArray(GetAODObjArrayName()+"Tracks");
1353
1354     //If too small or too large pt, skip
1355     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
1356
1357     // --- In case of redoing isolation from delta AOD ----
1358
1359     if(fMakeSeveralIC)
1360     {
1361       //Analysis of multiple IC at same time
1362       MakeSeveralICAnalysis(aod);
1363       continue;
1364     }
1365     else if(fReMakeIC)
1366     {
1367       //In case a more strict IC is needed in the produced AOD
1368       n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
1369       GetIsolationCut()->MakeIsolationCut(reftracks,   refclusters,
1370                                           GetReader(), GetCaloPID(),
1371                                           kFALSE, aod, "",
1372                                           n,nfrac,coneptsum, isolated);
1373       fhConeSumPt->Fill(pt,coneptsum);
1374       if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);
1375     }
1376
1377     // ---------------------------------------------------
1378
1379     //Fill pt distribution of particles in cone
1380     //Tracks
1381     coneptsum = 0;
1382     Double_t sumptFR = 0. ;
1383     TObjArray * trackList   = GetCTSTracks() ;
1384     for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
1385     {
1386       AliVTrack* track = (AliVTrack *) trackList->At(itrack);
1387       //fill the histograms at forward range
1388       if(!track)
1389       {
1390         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
1391         continue;
1392       }
1393
1394       Double_t dPhi = phi - track->Phi() + TMath::PiOver2();
1395       Double_t dEta = eta - track->Eta();
1396       Double_t arg  = dPhi*dPhi + dEta*dEta;
1397       if(TMath::Sqrt(arg) < conesize)
1398       {
1399         fhPtInFRCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
1400         sumptFR+=track->Pt();
1401       }
1402
1403       dPhi = phi - track->Phi() - TMath::PiOver2();
1404       arg  = dPhi*dPhi + dEta*dEta;
1405       if(TMath::Sqrt(arg) < conesize)
1406       {
1407         fhPtInFRCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
1408         sumptFR+=track->Pt();
1409       }
1410     }
1411
1412     fhFRConeSumPt->Fill(pt,sumptFR);
1413     if(reftracks)
1414     {
1415       for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
1416       {
1417         AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
1418         fhPtInCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
1419         coneptsum+=track->Pt();
1420       }
1421     }
1422
1423     //CaloClusters
1424     if(refclusters)
1425     {
1426       TLorentzVector mom ;
1427       for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
1428       {
1429         AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
1430         calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
1431
1432         fhPtInCone->Fill(pt, mom.Pt());
1433         coneptsum+=mom.Pt();
1434       }
1435     }
1436
1437     if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f\n", iaod, coneptsum);
1438
1439     if(!fReMakeIC) fhConeSumPt->Fill(pt,coneptsum);
1440
1441     Int_t mcTag = aod->GetTag() ;
1442     Int_t clID  = aod->GetCaloLabel(0) ;
1443
1444     if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeAnalysisFillHistograms() - pt %1.1f, eta %1.1f, phi %1.1f\n",pt, eta, phi);
1445
1446     if(isolation)
1447     {
1448       if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED: fill histograms\n", iaod);
1449
1450       FillTrackMatchingShowerShapeControlHistograms(clID,aod->GetFiducialArea(),mcTag);
1451
1452       fhEIso      ->Fill(energy);
1453       fhPtIso     ->Fill(pt);
1454       fhPhiIso    ->Fill(pt,phi);
1455       fhEtaIso    ->Fill(pt,eta);
1456       fhEtaPhiIso ->Fill(eta,phi);
1457
1458       if(decay)
1459       {
1460         fhPtDecayIso->Fill(pt);
1461         fhEtaPhiDecayIso->Fill(eta,phi);
1462       }
1463
1464       if(IsDataMC())
1465       {
1466
1467         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
1468         {
1469           fhPtIsoMCPhoton  ->Fill(pt);
1470         }
1471
1472         if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt))
1473         {
1474           fhPtIsoPrompt  ->Fill(pt);
1475           fhPhiIsoPrompt ->Fill(pt,phi);
1476           fhEtaIsoPrompt ->Fill(pt,eta);
1477         }
1478         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation))
1479         {
1480           fhPtIsoFragmentation  ->Fill(pt);
1481           fhPhiIsoFragmentation ->Fill(pt,phi);
1482           fhEtaIsoFragmentation ->Fill(pt,eta);
1483         }
1484         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))
1485         {
1486           fhPtIsoPi0Decay  ->Fill(pt);
1487           fhPhiIsoPi0Decay ->Fill(pt,phi);
1488           fhEtaIsoPi0Decay ->Fill(pt,eta);
1489         }
1490         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))
1491         {
1492           fhPtIsoEtaDecay  ->Fill(pt);
1493           fhPhiIsoEtaDecay ->Fill(pt,phi);
1494           fhEtaIsoEtaDecay ->Fill(pt,eta);
1495         }
1496         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))
1497         {
1498           fhPtIsoOtherDecay  ->Fill(pt);
1499           fhPhiIsoOtherDecay ->Fill(pt,phi);
1500           fhEtaIsoOtherDecay ->Fill(pt,eta);
1501         }
1502         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))
1503         {
1504           fhPtIsoConversion  ->Fill(pt);
1505           fhPhiIsoConversion ->Fill(pt,phi);
1506           fhEtaIsoConversion ->Fill(pt,eta);
1507         }
1508         else // anything else
1509         {
1510           fhPtIsoUnknown  ->Fill(pt);
1511           fhPhiIsoUnknown ->Fill(pt,phi);
1512           fhEtaIsoUnknown ->Fill(pt,eta);
1513         }
1514       }//Histograms with MC
1515
1516     }//Isolated histograms
1517
1518     if(!isolation)
1519     {
1520       if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d NOT ISOLATED, fill histograms\n", iaod);
1521
1522       fhPtNoIso    ->Fill(pt);
1523       fhEtaPhiNoIso->Fill(eta,phi);
1524
1525       if(decay)
1526       {
1527         fhPtDecayNoIso    ->Fill(pt);
1528         fhEtaPhiDecayNoIso->Fill(eta,phi);
1529       }
1530
1531       if(IsDataMC())
1532       {
1533         if     (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))        fhPtNoIsoMCPhoton     ->Fill(pt);
1534         if     (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtNoIsoPi0Decay     ->Fill(pt);
1535         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtNoIsoEtaDecay     ->Fill(pt);
1536         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtNoIsoOtherDecay   ->Fill(pt);
1537         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt))        fhPtNoIsoPrompt       ->Fill(pt);
1538         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation)) fhPtNoIsoFragmentation->Fill(pt);
1539         else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))    fhPtNoIsoConversion   ->Fill(pt);
1540         else                                                                                   fhPtNoIsoUnknown      ->Fill(pt);
1541       }
1542     }
1543
1544   }// aod loop
1545
1546 }
1547
1548
1549 //_____________________________________________________________________________________
1550 void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation* ph)
1551 {
1552
1553   //Isolation Cut Analysis for both methods and different pt cuts and cones
1554   Float_t ptC   = ph->Pt();
1555   Float_t etaC  = ph->Eta();
1556   Float_t phiC  = ph->Phi();
1557   Int_t   tag   = ph->GetTag();
1558   Bool_t  decay = ph->IsTagged();
1559
1560   if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeSeveralICAnalysis() - Isolate pT %2.2f\n",ptC);
1561
1562   //Keep original setting used when filling AODs, reset at end of analysis
1563   Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
1564   Float_t ptfracorg  = GetIsolationCut()->GetPtFraction();
1565   Float_t rorg       = GetIsolationCut()->GetConeSize();
1566
1567   Float_t coneptsum = 0 ;
1568   Int_t   n    [10][10];//[fNCones][fNPtThresFrac];
1569   Int_t   nfrac[10][10];//[fNCones][fNPtThresFrac];
1570   Bool_t  isolated   = kFALSE;
1571   Int_t   n_cone;
1572   Int_t   nfrac_cone;
1573
1574
1575   // fill hist with all particles before isolation criteria
1576   fhPtNoIso    ->Fill(ptC);
1577   fhEtaPhiNoIso->Fill(etaC,phiC);
1578
1579   if(IsDataMC())
1580   {
1581     if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))        fhPtNoIsoMCPhoton     ->Fill(ptC);
1582     if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtNoIsoPi0Decay     ->Fill(ptC);
1583     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtNoIsoEtaDecay     ->Fill(ptC);
1584     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtNoIsoOtherDecay   ->Fill(ptC);
1585     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtNoIsoPrompt       ->Fill(ptC);
1586     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtNoIsoFragmentation->Fill(ptC);
1587     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtNoIsoConversion   ->Fill(ptC);
1588     else                                                                                 fhPtNoIsoUnknown      ->Fill(ptC);
1589   }
1590
1591   if(decay)
1592   {
1593     fhPtDecayNoIso    ->Fill(ptC);
1594     fhEtaPhiDecayNoIso->Fill(etaC,phiC);
1595   }
1596   //Get vertex for photon momentum calculation
1597   Double_t vertex[] = {0,0,0} ; //vertex ;
1598   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1599     {
1600       GetReader()->GetVertex(vertex);
1601     }
1602
1603   //Loop on cone sizes
1604   for(Int_t icone = 0; icone<fNCones; icone++)
1605   {
1606     //Recover reference arrays with clusters and tracks
1607     TObjArray * refclusters = ph->GetObjArray(GetAODObjArrayName()+"Clusters");
1608     TObjArray * reftracks   = ph->GetObjArray(GetAODObjArrayName()+"Tracks");
1609
1610     //If too small or too large pt, skip
1611     if(ptC < GetMinPt() || ptC > GetMaxPt() ) continue ;
1612
1613    //In case a more strict IC is needed in the produced AOD
1614
1615   n_cone=0; nfrac_cone = 0; isolated = kFALSE; coneptsum = 0;
1616
1617   GetIsolationCut()->SetSumPtThreshold(100);
1618   GetIsolationCut()->SetPtThreshold(100);
1619   GetIsolationCut()->SetPtFraction(100);
1620   GetIsolationCut()->SetConeSize(fConeSizes[icone]);
1621   GetIsolationCut()->MakeIsolationCut(reftracks,   refclusters,
1622                                           GetReader(), GetCaloPID(),
1623                                           kFALSE, ph, "",
1624                                           n_cone,nfrac_cone,coneptsum, isolated);
1625
1626
1627     fhSumPtLeadingPt[icone]->Fill(ptC,coneptsum);
1628
1629     // retreive pt tracks to fill histo vs. pt leading
1630     //Fill pt distribution of particles in cone
1631     //fhPtLeadingPt(),fhFRSumPtLeadingPt(),fhFRPtLeadingPt(),
1632
1633     //Tracks
1634     coneptsum = 0;
1635     Double_t sumptFR = 0. ;
1636     TObjArray * trackList   = GetCTSTracks() ;
1637     for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
1638     {
1639       AliVTrack* track = (AliVTrack *) trackList->At(itrack);
1640       //fill the histograms at forward range
1641       if(!track)
1642       {
1643         printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
1644         continue;
1645       }
1646
1647       Double_t dPhi = phiC - track->Phi() + TMath::PiOver2();
1648       Double_t dEta = etaC - track->Eta();
1649       Double_t arg  = dPhi*dPhi + dEta*dEta;
1650       if(TMath::Sqrt(arg) < fConeSizes[icone])
1651       {
1652         fhFRPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
1653         sumptFR+=track->Pt();
1654       }
1655
1656       dPhi = phiC - track->Phi() - TMath::PiOver2();
1657       arg  = dPhi*dPhi + dEta*dEta;
1658       if(TMath::Sqrt(arg) < fConeSizes[icone])
1659       {
1660         fhFRPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
1661         sumptFR+=track->Pt();
1662       }
1663     }
1664     fhFRSumPtLeadingPt[icone]->Fill(ptC,sumptFR);
1665     if(reftracks)
1666     {
1667       for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
1668       {
1669         AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
1670         fhPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
1671         coneptsum+=track->Pt();
1672       }
1673     }
1674     //CaloClusters
1675     if(refclusters)
1676     {
1677       TLorentzVector mom ;
1678       for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
1679       {
1680         AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
1681         calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
1682
1683         fhPtLeadingPt[icone]->Fill(ptC, mom.Pt());
1684         coneptsum+=mom.Pt();
1685       }
1686     }
1687     ///////////////////
1688
1689
1690     //Loop on ptthresholds
1691     for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++)
1692     {
1693       n    [icone][ipt]=0;
1694       nfrac[icone][ipt]=0;
1695       GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
1696       GetIsolationCut()->SetPtFraction(fPtFractions[ipt]) ;
1697       GetIsolationCut()->SetSumPtThreshold(fSumPtThresholds[ipt]);
1698
1699       GetIsolationCut()->MakeIsolationCut(reftracks, refclusters,
1700                                           GetReader(), GetCaloPID(),
1701                                           kFALSE, ph, "",
1702                                           n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
1703
1704       if(!isolated) continue;
1705       //Normal ptThreshold cut
1706
1707       if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - cone size %1.1f, ptThres  %1.1f, sumptThresh  %1.1f, n %d, nfrac %d, coneptsum %2.2f, isolated %d\n",
1708                                 fConeSizes[icone],fPtThresholds[ipt],fSumPtThresholds[ipt],n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
1709       if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - pt %1.1f, eta %1.1f, phi %1.1f\n",ptC, etaC, phiC);
1710
1711       if(n[icone][ipt] == 0)
1712       {
1713         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling pt threshold loop\n");
1714         fhPtThresIsolated[icone][ipt]->Fill(ptC);
1715         fhEtaPhiPtThresIso[icone][ipt]->Fill(etaC,phiC);
1716
1717         if(decay)
1718         {
1719           fhPtPtThresDecayIso[icone][ipt]->Fill(ptC);
1720           //      fhEtaPhiPtThresDecayIso[icone][ipt]->Fill(etaC,phiC);
1721         }
1722
1723         if(IsDataMC())
1724         {
1725           if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtThresIsolatedPrompt[icone][ipt]       ->Fill(ptC) ;
1726           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtThresIsolatedConversion[icone][ipt]   ->Fill(ptC) ;
1727           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
1728           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtThresIsolatedPi0Decay[icone][ipt]     ->Fill(ptC) ;
1729           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtThresIsolatedEtaDecay[icone][ipt]     ->Fill(ptC) ;
1730           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtThresIsolatedOtherDecay[icone][ipt]   ->Fill(ptC) ;
1731           else                                                                                  fhPtThresIsolatedUnknown[icone][ipt]      ->Fill(ptC) ;
1732         }
1733       }
1734
1735       // pt in cone fraction
1736       if(nfrac[icone][ipt] == 0)
1737       {
1738         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling frac loop\n");
1739         fhPtFracIsolated[icone][ipt]->Fill(ptC);
1740         fhEtaPhiPtFracIso[icone][ipt]->Fill(etaC,phiC);
1741
1742         if(decay)
1743         {
1744           fhPtPtFracDecayIso[icone][ipt]->Fill(ptC);
1745           fhEtaPhiPtFracDecayIso[icone][ipt]->Fill(etaC,phiC);
1746         }
1747
1748         if(IsDataMC())
1749         {
1750           if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtFracIsolatedPrompt[icone][ipt]       ->Fill(ptC) ;
1751           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtFracIsolatedConversion[icone][ipt]   ->Fill(ptC) ;
1752           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
1753           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtFracIsolatedPi0Decay[icone][ipt]     ->Fill(ptC) ;
1754           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtFracIsolatedEtaDecay[icone][ipt]     ->Fill(ptC) ;
1755           else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtFracIsolatedOtherDecay[icone][ipt]   ->Fill(ptC) ;
1756           else  fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC) ;
1757         }
1758       }
1759
1760       if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - checking IC method : %i\n",GetIsolationCut()->GetICMethod());
1761
1762       //Pt threshold on pt cand/ sum in cone histograms
1763       if(coneptsum<fSumPtThresholds[ipt])
1764       {//      if((GetIsolationCut()->GetICMethod())==1){//kSumPtIC){
1765         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling sum loop\n");
1766         fhPtSumIsolated[icone][ipt]->Fill(ptC) ;
1767         fhEtaPhiPtSumIso[icone][ipt]->Fill(etaC, phiC) ;
1768         if(decay)
1769         {
1770           fhPtPtSumDecayIso[icone][ipt]->Fill(ptC);
1771           fhEtaPhiPtSumDecayIso[icone][ipt]->Fill(etaC, phiC) ;
1772         }
1773       }
1774
1775
1776
1777
1778     // pt sum pt frac method
1779     if( ((fPtFractions[ipt]*ptC < fSumPtThresholds[ipt]) && (coneptsum < fSumPtThresholds[ipt])) || ((fPtFractions[ipt]*ptC > fSumPtThresholds[ipt]) && (coneptsum < fPtFractions[ipt]*ptC)) )
1780       {
1781         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling PtFrac PtSum loop\n");
1782         fhPtFracPtSumIso[icone][ipt]->Fill(ptC) ;
1783         fhEtaPhiFracPtSumIso[icone][ipt]->Fill(etaC,phiC) ;
1784
1785         if(decay)
1786         {
1787           fhPtFracPtSumDecayIso[icone][ipt]->Fill(ptC);
1788           fhEtaPhiFracPtSumDecayIso[icone][ipt]->Fill(etaC,phiC);
1789         }
1790       }
1791
1792       // density method
1793       Float_t cellDensity = GetIsolationCut()->GetCellDensity( ph, GetReader());
1794       if(coneptsum<fSumPtThresholds[ipt]*cellDensity)
1795       {//(GetIsolationCut()->GetICMethod())==4){//kSumDensityIC) {
1796         if(GetDebug()>0) printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling density loop\n");
1797         fhPtSumDensityIso[icone][ipt]->Fill(ptC) ;
1798         fhEtaPhiSumDensityIso[icone][ipt]->Fill(etaC,phiC) ;
1799
1800         if(decay)
1801         {
1802           fhPtSumDensityDecayIso[icone][ipt]->Fill(ptC);
1803           fhEtaPhiSumDensityDecayIso[icone][ipt]->Fill(etaC, phiC);
1804         }
1805
1806       }
1807     }//pt thresh loop
1808
1809     if(IsDataMC())
1810     {
1811       if     ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))        fhPtSumIsolatedPrompt[icone]       ->Fill(ptC,coneptsum) ;
1812       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))    fhPtSumIsolatedConversion[icone]   ->Fill(ptC,coneptsum) ;
1813       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
1814       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))      fhPtSumIsolatedPi0Decay[icone]     ->Fill(ptC,coneptsum) ;
1815       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))      fhPtSumIsolatedEtaDecay[icone]     ->Fill(ptC,coneptsum) ;
1816       else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))    fhPtSumIsolatedOtherDecay[icone]   ->Fill(ptC,coneptsum) ;
1817       else  fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;
1818     }
1819
1820   }//cone size loop
1821
1822   //Reset original parameters for AOD analysis
1823   GetIsolationCut()->SetPtThreshold(ptthresorg);
1824   GetIsolationCut()->SetPtFraction(ptfracorg);
1825   GetIsolationCut()->SetConeSize(rorg);
1826
1827 }
1828
1829 //_____________________________________________________________
1830 void AliAnaParticleIsolation::Print(const Option_t * opt) const
1831 {
1832
1833   //Print some relevant parameters set for the analysis
1834   if(! opt)
1835     return;
1836
1837   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1838   AliAnaCaloTrackCorrBaseClass::Print(" ");
1839
1840   printf("ReMake Isolation          = %d \n",  fReMakeIC) ;
1841   printf("Make Several Isolation    = %d \n",  fMakeSeveralIC) ;
1842   printf("Calorimeter for isolation = %s \n",  fCalorimeter.Data()) ;
1843
1844   if(fMakeSeveralIC)
1845   {
1846     printf("N Cone Sizes       =     %d\n", fNCones) ;
1847     printf("Cone Sizes          =    \n") ;
1848     for(Int_t i = 0; i < fNCones; i++)
1849       printf("  %1.2f;",  fConeSizes[i]) ;
1850     printf("    \n") ;
1851
1852     printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
1853     printf(" pT thresholds         =    \n") ;
1854     for(Int_t i = 0; i < fNPtThresFrac; i++)
1855       printf("   %2.2f;",  fPtThresholds[i]) ;
1856
1857     printf("    \n") ;
1858
1859     printf(" pT fractions         =    \n") ;
1860     for(Int_t i = 0; i < fNPtThresFrac; i++)
1861       printf("   %2.2f;",  fPtFractions[i]) ;
1862
1863     printf("    \n") ;
1864
1865     printf("sum pT thresholds         =    \n") ;
1866     for(Int_t i = 0; i < fNPtThresFrac; i++)
1867       printf("   %2.2f;",  fSumPtThresholds[i]) ;
1868
1869
1870   }
1871
1872   printf("Histograms: %3.1f < pT sum < %3.1f,  Nbin = %d\n",    fHistoPtSumMin,    fHistoPtSumMax,    fHistoNPtSumBins   );
1873   printf("Histograms: %3.1f < pT in cone < %3.1f, Nbin = %d\n", fHistoPtInConeMin, fHistoPtInConeMax, fHistoNPtInConeBins);
1874
1875   printf("    \n") ;
1876
1877 }
1878