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