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