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