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