PWGGA/CaloTrackCorrelations/AliAnaGeneratorKine.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 is 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 // Do direct photon/decay photon (eta, pi0, other)/pi0/eta isolation
  18 // and correlation with partons/jets/hadrons analysis at the generator level.
  19 // For MC kinematics at ESD and AOD level.
  20 // Jets only considered in the case of Pythia, check what to do with others.
  21 //
  22 // -- Author: Gustavo Conesa (LPSC-CNRS-Grenoble)
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25 // --- ROOT system ---
  26 #include "TH2F.h"
  27 #include "TParticle.h"
  28 #include "TDatabasePDG.h"
  29
  30 //---- ANALYSIS system ----
  31 #include "AliAnaGeneratorKine.h"
  32 #include "AliStack.h"
  33 #include "AliGenPythiaEventHeader.h"
  34 #include "AliAODMCParticle.h"
  35
  36 ClassImp(AliAnaGeneratorKine)
  37
  38
  39 //__________________________________________
  40 AliAnaGeneratorKine::AliAnaGeneratorKine() :
  41 AliAnaCaloTrackCorrBaseClass(),
  42 fTriggerDetector(),  fTriggerDetectorString(),
  43 fFidCutTrigger(0),
  44 fMinChargedPt(0),    fMinNeutralPt(0),
  45 fStack(0),           fAODMCparticles(0),
  46 //fParton2(0),         fParton3(0),
  47 fParton6(),          fParton7(),
  48 fParton6PDG(0),      fParton7PDG(0),
  49 fJet6(),             fJet7(),
  50 fTrigger(),          fLVTmp(),
  51 fNPrimaries(0),      fPtHard(0),
  52 fhPtHard(0),         fhPtParton(0),    fhPtJet(0),
  53 fhPtPartonPtHard(0), fhPtJetPtHard(0), fhPtJetPtParton(0)
  54 {
  55   //Default Ctor
  56
  57   //Initialize parameters
  58   InitParameters();
  59
  60
  61
  62   for(Int_t p = 0; p < fgkNmcPrimTypes; p++)
  63   {
  64     fhPt[p] = 0;
  65
  66     for(Int_t i = 0; i < fgkNIso; i++)
  67     {
  68       fhPtLeading[p][i]          = 0;
  69       fhPtLeadingSumPt[p][i]     = 0;
  70       fhPtLeadingIsolated[p][i]  = 0;
  71       for(Int_t j = 0; j < 2; j++)
  72       {
  73         fhZHard[p][j][i]           = 0;
  74         fhZHardIsolated[p][j][i]   = 0;
  75         fhZParton[p][j][i]         = 0;
  76         fhZPartonIsolated[p][j][i] = 0;
  77         fhZJet[p][j][i]            = 0;
  78         fhZJetIsolated[p][j][i]    = 0;
  79         fhXE[p][j][i]              = 0;
  80         fhXEIsolated[p][j][i]      = 0;
  81         fhXEUE[p][j][i]            = 0;
  82         fhXEUEIsolated[p][j][i]    = 0;
  83
  84         fhPtPartonTypeNear[p][j][i]         = 0;
  85         fhPtPartonTypeNearIsolated[p][j][i] = 0;
  86
  87         fhPtPartonTypeAway[p][j][i]         = 0;
  88         fhPtPartonTypeAwayIsolated[p][j][i] = 0;
  89
  90         if( p == 0 )  fhPtAcceptedGammaJet[j][i] = 0;
  91       }
  92     }
  93   }
  94
  95 }
  96
  97 //___________________________________________________________________________
  98 Bool_t  AliAnaGeneratorKine::CorrelateWithPartonOrJet(Int_t   indexTrig,
  99                                                       Int_t   partType,
 100                                                       Bool_t  leading [fgkNIso],
 101                                                       Bool_t  isolated[fgkNIso],
 102                                                       Int_t & iparton )
 103 {
 104   // Correlate trigger with partons or jets, get z
 105
 106   AliDebug(1,"Start");
 107
 108   if( fNPrimaries < 7 )
 109   {
 110     AliDebug(1,Form("End, not enough partons, n primaries %d",fNPrimaries));
 111     return kFALSE ;
 112   }
 113
 114   //
 115   // Get the index of the mother
 116   //
 117
 118   if(fStack) // ESD
 119   {
 120     iparton =  (fStack->Particle(indexTrig))->GetFirstMother();
 121     TParticle * mother = fStack->Particle(iparton);
 122     while (iparton > 7)
 123     {
 124       iparton   = mother->GetFirstMother();
 125       if(iparton < 0)
 126       {
 127         AliWarning("Negative index, skip event");
 128         return kFALSE;
 129       }
 130       mother = fStack->Particle(iparton);
 131     }
 132   }
 133   else // AOD
 134   {
 135     iparton =  ((AliAODMCParticle*) fAODMCparticles->At(indexTrig))->GetMother();
 136     AliAODMCParticle * mother = (AliAODMCParticle*) fAODMCparticles->At(iparton);
 137     while (iparton > 7)
 138     {
 139       iparton   = mother->GetMother();
 140       if(iparton < 0)
 141       {
 142         AliWarning("Negative index, skip event");
 143         return kFALSE;
 144       }
 145       mother = (AliAODMCParticle*) fAODMCparticles->At(iparton);
 146     }
 147   }
 148
 149   //printf("Mother is parton %d with pdg %d\n",imom,mother->GetPdgCode());
 150
 151   if(iparton < 6)
 152   {
 153     AliDebug(1,Form("This particle is not from hard process - pdg %d, parton index %d\n",partType, iparton));
 154     return kFALSE;
 155   }
 156
 157   //
 158   // Get the kinematics
 159   //
 160
 161   Float_t  ptTrig  = fTrigger.Pt();
 162   Float_t phiTrig  = fTrigger.Phi();
 163   Float_t etaTrig  = fTrigger.Eta();
 164
 165   AliDebug(1,Form("Trigger pdg %d pT %2.2f, phi %2.2f, eta %2.2f", partType,ptTrig,phiTrig*TMath::RadToDeg(),etaTrig));
 166
 167   Float_t jetPt    = fJet6.Pt();
 168   Float_t jetPhi   = fJet6.Phi();
 169   Float_t jetEta   = fJet6.Eta();
 170
 171   AliDebug(1,Form("Jet 6 pT %2.2f, phi %2.2f, eta %2.2f",jetPt,jetPhi*TMath::RadToDeg(),jetEta));
 172
 173   Float_t awayJetPt  = fJet7.Pt();
 174   Float_t awayJetEta = fJet7.Eta();
 175   Float_t awayJetPhi = fJet7.Phi();
 176
 177   AliDebug(1,Form("Jet 7 pT %2.2f, phi %2.2f, eta %2.2f",awayJetPt,awayJetPhi*TMath::RadToDeg(),awayJetEta));
 178
 179   Float_t partonPt = fParton6.Pt();
 180
 181   Int_t nearPDG = fParton6PDG;
 182   Int_t awayPDG = fParton7PDG;
 183
 184   AliDebug(1,Form("Parton6 pT pT %2.2f, pdg %d",fParton6.Pt(),fParton6PDG));
 185   AliDebug(1,Form("Parton7 pT pT %2.2f, pdg %d",fParton7.Pt(),fParton7PDG));
 186
 187   if ( iparton == 7 )
 188   {
 189     partonPt = fParton7.Pt();
 190
 191     jetPt  = fJet7.Pt();
 192     jetPhi = fJet7.Phi();
 193     jetEta = fJet7.Eta();
 194
 195     awayJetPt  = fJet6.Pt();
 196     awayJetEta = fJet6.Eta();
 197     awayJetPhi = fJet6.Phi();
 198
 199     nearPDG = fParton7PDG;
 200     awayPDG = fParton6PDG;
 201   }
 202
 203   Float_t deltaPhi = TMath::Abs(phiTrig-awayJetPhi) *TMath::RadToDeg();
 204   AliDebug(1,Form("Trigger Away jet phi %2.2f\n",deltaPhi));
 205
 206   //
 207   // Get id of parton in near and away side
 208   //
 209   Int_t away = -1;
 210   Int_t near = -1;
 211   if     (nearPDG == 22) near = 0;
 212   else if(nearPDG == 21) near = 1;
 213   else                   near = 2;
 214
 215   if     (awayPDG == 22) away = 0;
 216   else if(awayPDG == 21) away = 1;
 217   else                   away = 2;
 218
 219   // RATIOS
 220
 221   Float_t zHard = -1;
 222   Float_t zPart = -1;
 223   Float_t zJet  = -1;
 224
 225   if( fPtHard  > 0 ) zHard = ptTrig / fPtHard ;
 226   if( partonPt > 0 ) zPart = ptTrig / partonPt;
 227   if( jetPt    > 0 ) zJet  = ptTrig / jetPt   ;
 228
 229   //if(zHard > 1 ) printf("*** Particle energy larger than pT hard z=%f\n",zHard);
 230
 231   //printf("Z: hard %2.2f, parton %2.2f, jet %2.2f\n",zHard,zPart,zJet);
 232
 233   // conditions loop
 234   for( Int_t i = 0; i < fgkNIso ; i++ )
 235   {
 236     fhPtPartonTypeNear[partType][leading[i]][i]->Fill(ptTrig,near);
 237     fhPtPartonTypeAway[partType][leading[i]][i]->Fill(ptTrig,away);
 238
 239     fhZHard  [partType][leading[i]][i]->Fill(ptTrig,zHard);
 240     fhZParton[partType][leading[i]][i]->Fill(ptTrig,zPart);
 241     fhZJet   [partType][leading[i]][i]->Fill(ptTrig,zJet );
 242
 243     if(isolated[i])
 244     {
 245       fhPtPartonTypeNearIsolated[partType][leading[i]][i]->Fill(ptTrig,near);
 246       fhPtPartonTypeAwayIsolated[partType][leading[i]][i]->Fill(ptTrig,away);
 247
 248       fhZHardIsolated  [partType][leading[i]][i]->Fill(ptTrig,zHard);
 249       fhZPartonIsolated[partType][leading[i]][i]->Fill(ptTrig,zPart);
 250       fhZJetIsolated   [partType][leading[i]][i]->Fill(ptTrig,zJet);
 251     }
 252
 253     if(partType == kmcPrimPhoton && deltaPhi < 220 && deltaPhi > 140 && TMath::Abs(awayJetEta) < 0.6)
 254     {
 255       //printf("Accept jet\n");
 256       fhPtAcceptedGammaJet[leading[i]][i]->Fill(ptTrig,away);
 257     }
 258     //else printf("Reject jet!!!\n");
 259
 260   } // conditions loop
 261
 262   AliDebug(1,"End TRUE");
 263
 264   return kTRUE;
 265 }
 266
 267
 268 //____________________________________________________
 269 TList *  AliAnaGeneratorKine::GetCreateOutputObjects()
 270 {
 271   // Create histograms to be saved in output file
 272
 273   TList * outputContainer = new TList() ;
 274   outputContainer->SetName("GenKineHistos") ;
 275
 276   Int_t   nptbins    = GetHistogramRanges()->GetHistoPtBins();
 277   Float_t ptmax      = GetHistogramRanges()->GetHistoPtMax();
 278   Float_t ptmin      = GetHistogramRanges()->GetHistoPtMin();
 279
 280   Int_t   nptsumbins = GetHistogramRanges()->GetHistoNPtSumBins();
 281   Float_t ptsummax   = GetHistogramRanges()->GetHistoPtSumMax();
 282   Float_t ptsummin   = GetHistogramRanges()->GetHistoPtSumMin();
 283
 284
 285   fhPtHard  = new TH1F("hPtHard"," pt hard for selected triggers",nptbins,ptmin,ptmax);
 286   fhPtHard->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
 287   outputContainer->Add(fhPtHard);
 288
 289   fhPtParton  = new TH1F("hPtParton"," pt parton for selected triggers",nptbins,ptmin,ptmax);
 290   fhPtParton->SetXTitle("#it{p}_{T}^{parton} (GeV/#it{c})");
 291   outputContainer->Add(fhPtParton);
 292
 293   fhPtJet  = new TH1F("hPtJet"," pt jet for selected triggers",nptbins,ptmin,ptmax);
 294   fhPtJet->SetXTitle("#it{p}_{T}^{jet} (GeV/#it{c})");
 295   outputContainer->Add(fhPtJet);
 296
 297   fhPtPartonPtHard  = new TH2F("hPtPartonPtHard","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
 298   fhPtPartonPtHard->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
 299   fhPtPartonPtHard->SetYTitle("#it{p}_{T}^{parton}/#it{p}_{T}^{hard}");
 300   outputContainer->Add(fhPtPartonPtHard);
 301
 302   fhPtJetPtHard  = new TH2F("hPtJetPtHard","jet pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
 303   fhPtJetPtHard->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
 304   fhPtJetPtHard->SetYTitle("#it{p}_{T}^{jet}/#it{p}_{T}^{hard}");
 305   outputContainer->Add(fhPtJetPtHard);
 306
 307   fhPtJetPtParton  = new TH2F("hPtJetPtParton","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
 308   fhPtJetPtParton->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
 309   fhPtJetPtParton->SetYTitle("#it{p}_{T}^{jet}/#it{p}_{T}^{parton}");
 310   outputContainer->Add(fhPtJetPtParton);
 311
 312   TString name   [] = {"","_EMC","_Photon","_EMC_Photon"};
 313   TString title  [] = {"",", neutral in EMCal",", neutral only #gamma-like",", neutral in EMCal and only #gamma-like"};
 314   TString leading[] = {"NotLeading","Leading"};
 315
 316   TString partTitl[] = {"#gamma_{direct}","#gamma_{decay}^{#pi}","#gamma_{decay}^{#eta}","#gamma_{decay}^{other}","#pi^{0}","#eta"};
 317   TString particle[] = {"DirectPhoton"   ,"Pi0DecayPhoton"      ,"EtaDecayPhoton"       ,"OtherDecayPhoton"      ,"Pi0"    ,"Eta"};
 318
 319   for(Int_t p = 0; p < fgkNmcPrimTypes; p++)
 320   {
 321     fhPt[p]  = new TH1F(Form("h%sPt",particle[p].Data()),Form("Input %s p_{T}",partTitl[p].Data()),nptbins,ptmin,ptmax);
 322     fhPt[p]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 323     outputContainer->Add(fhPt[p]);
 324
 325     for(Int_t i = 0; i < fgkNIso; i++)
 326     {
 327       // Pt
 328
 329       fhPtLeading[p][i]  = new TH1F(Form("h%sPtLeading%s", particle[p].Data(), name[i].Data()),
 330                                     Form("%s: Leading of all particles%s",partTitl[p].Data(),title[i].Data()),
 331                                     nptbins,ptmin,ptmax);
 332       fhPtLeading[p][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 333       outputContainer->Add(fhPtLeading[p][i]);
 334
 335       fhPtLeadingIsolated[p][i]  = new TH1F(Form("h%sPtLeadingIsolated%s", particle[p].Data(), name[i].Data()),
 336                                             Form("%s: Leading of all particles%s, isolated",partTitl[p].Data(),title[i].Data()),
 337                                             nptbins,ptmin,ptmax);
 338       fhPtLeadingIsolated[p][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 339       outputContainer->Add(fhPtLeadingIsolated[p][i]);
 340
 341       fhPtLeadingSumPt[p][i]  = new TH2F(Form("h%sPtLeadingSumPt%s", particle[p].Data(), name[i].Data()),
 342                                          Form("%s: Leading of all particles%s",partTitl[p].Data(),title[i].Data()),
 343                                          nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
 344       fhPtLeadingSumPt[p][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 345       fhPtLeadingSumPt[p][i]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
 346       outputContainer->Add(fhPtLeadingSumPt[p][i]);
 347
 348
 349       // Leading or not loop
 350       for(Int_t j = 0; j < fgkNLead; j++)
 351       {
 352         if(p==0)
 353         {
 354           fhPtAcceptedGammaJet[j][i]  = new TH2F(Form("hPtAcceptedGammaJet%s%s",           leading[j].Data(), name[i].Data()),
 355                                                  Form("#gamma-jet: %s of all particles%s", leading[j].Data(), title[i].Data()),
 356                                                  nptbins,ptmin,ptmax,3,0,3);
 357           fhPtAcceptedGammaJet[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 358           fhPtAcceptedGammaJet[j][i]->SetYTitle("Parton type");
 359           fhPtAcceptedGammaJet[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
 360           fhPtAcceptedGammaJet[j][i]->GetYaxis()->SetBinLabel(2,"g");
 361           fhPtAcceptedGammaJet[j][i]->GetYaxis()->SetBinLabel(3,"q");
 362           outputContainer->Add(fhPtAcceptedGammaJet[j][i]);
 363         }
 364         // Near side parton
 365
 366         fhPtPartonTypeNear[p][j][i]  = new TH2F(Form("h%sPtPartonTypeNear%s%s",   particle[p].Data(), leading[j].Data(), name[i].Data()),
 367                                                 Form("%s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 368                                                 nptbins,ptmin,ptmax,3,0,3);
 369         fhPtPartonTypeNear[p][j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 370         fhPtPartonTypeNear[p][j][i]->SetYTitle("Parton type");
 371         fhPtPartonTypeNear[p][j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
 372         fhPtPartonTypeNear[p][j][i]->GetYaxis()->SetBinLabel(2,"g");
 373         fhPtPartonTypeNear[p][j][i]->GetYaxis()->SetBinLabel(3,"q");
 374         outputContainer->Add(fhPtPartonTypeNear[p][j][i]);
 375
 376         fhPtPartonTypeNearIsolated[p][j][i]  = new TH2F(Form("h%sPtPartonTypeNear%sIsolated%s",     particle[p].Data(), leading[j].Data(), name[i].Data()),
 377                                                         Form("%s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 378                                                         nptbins,ptmin,ptmax,3,0,3);
 379         fhPtPartonTypeNearIsolated[p][j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 380         fhPtPartonTypeNearIsolated[p][j][i]->SetYTitle("Parton type");
 381         fhPtPartonTypeNearIsolated[p][j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
 382         fhPtPartonTypeNearIsolated[p][j][i]->GetYaxis()->SetBinLabel(2,"g");
 383         fhPtPartonTypeNearIsolated[p][j][i]->GetYaxis()->SetBinLabel(3,"q");
 384         outputContainer->Add(fhPtPartonTypeNearIsolated[p][j][i]);
 385
 386
 387         // Away side parton
 388
 389         fhPtPartonTypeAway[p][j][i]  = new TH2F(Form("h%sPtPartonTypeAway%s%s",   particle[p].Data(), leading[j].Data(), name[i].Data()),
 390                                                 Form("%s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 391                                                 nptbins,ptmin,ptmax,3,0,3);
 392         fhPtPartonTypeAway[p][j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 393         fhPtPartonTypeAway[p][j][i]->SetYTitle("Parton type");
 394         fhPtPartonTypeAway[p][j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
 395         fhPtPartonTypeAway[p][j][i]->GetYaxis()->SetBinLabel(2,"g");
 396         fhPtPartonTypeAway[p][j][i]->GetYaxis()->SetBinLabel(3,"q");
 397         outputContainer->Add(fhPtPartonTypeAway[p][j][i]);
 398
 399         fhPtPartonTypeAwayIsolated[p][j][i]  = new TH2F(Form("h%sPtPartonTypeAway%sIsolated%s",     particle[p].Data(), leading[j].Data(), name[i].Data()),
 400                                                         Form("%s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 401                                                         nptbins,ptmin,ptmax,3,0,3);
 402         fhPtPartonTypeAwayIsolated[p][j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 403         fhPtPartonTypeAwayIsolated[p][j][i]->SetYTitle("Parton type");
 404         fhPtPartonTypeAwayIsolated[p][j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
 405         fhPtPartonTypeAwayIsolated[p][j][i]->GetYaxis()->SetBinLabel(2,"g");
 406         fhPtPartonTypeAwayIsolated[p][j][i]->GetYaxis()->SetBinLabel(3,"q");
 407         outputContainer->Add(fhPtPartonTypeAwayIsolated[p][j][i]);
 408
 409         // zHard
 410
 411         fhZHard[p][j][i]  = new TH2F(Form("h%sZHard%s%s",                               particle[p].Data(), leading[j].Data(), name[i].Data()),
 412                                      Form("#it{z}_{Hard} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 413                                      nptbins,ptmin,ptmax,200,0,2);
 414         fhZHard[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 415         fhZHard[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 416         outputContainer->Add(fhZHard[p][j][i]);
 417
 418         fhZHardIsolated[p][j][i]  = new TH2F(Form("h%sZHard%sIsolated%s",                                 particle[p].Data(), leading[j].Data(), name[i].Data()),
 419                                              Form("#it{z}_{Hard} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 420                                              nptbins,ptmin,ptmax,200,0,2);
 421         fhZHardIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 422         fhZHardIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 423         outputContainer->Add(fhZHardIsolated[p][j][i]);
 424
 425         // zHard
 426
 427         fhZParton[p][j][i]  = new TH2F(Form("h%sZParton%s%s",                               particle[p].Data(), leading[j].Data(), name[i].Data()),
 428                                        Form("#it{z}_{Parton} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 429                                        nptbins,ptmin,ptmax,200,0,2);
 430         fhZParton[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 431         fhZParton[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 432         outputContainer->Add(fhZParton[p][j][i]);
 433
 434         fhZPartonIsolated[p][j][i]  = new TH2F(Form("h%sZParton%sIsolated%s",                                 particle[p].Data(), leading[j].Data(), name[i].Data()),
 435                                                Form("#it{z}_{Parton} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 436                                                nptbins,ptmin,ptmax,200,0,2);
 437         fhZPartonIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 438         fhZPartonIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 439         outputContainer->Add(fhZPartonIsolated[p][j][i]);
 440
 441
 442         // zJet
 443
 444         fhZJet[p][j][i]  = new TH2F(Form("h%sZJet%s%s",                               particle[p].Data(), leading[j].Data(), name[i].Data()),
 445                                     Form("#it{z}_{Jet} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 446                                     nptbins,ptmin,ptmax,200,0,2);
 447         fhZJet[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 448         fhZJet[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 449         outputContainer->Add(fhZJet[p][j][i]);
 450
 451         fhZJetIsolated[p][j][i]  = new TH2F(Form("h%sZJet%sIsolated%s",                                 particle[p].Data(), leading[j].Data(), name[i].Data()),
 452                                             Form("#it{z}_{Jet} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 453                                             nptbins,ptmin,ptmax,200,0,2);
 454         fhZJetIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 455         fhZJetIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 456         outputContainer->Add(fhZJetIsolated[p][j][i]);
 457
 458
 459         // XE
 460
 461         fhXE[p][j][i]  = new TH2F(Form("h%sXE%s%s",                                 particle[p].Data(), leading[j].Data(), name[i].Data()),
 462                                   Form("#it{z}_{Jet} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 463                                   nptbins,ptmin,ptmax,200,0,2);
 464         fhXE[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 465         fhXE[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 466         outputContainer->Add(fhXE[p][j][i]);
 467
 468         fhXEIsolated[p][j][i]  = new TH2F(Form("h%sXE%sIsolated%s",                                   particle[p].Data(), leading[j].Data(), name[i].Data()),
 469                                           Form("#it{z}_{Jet} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 470                                           nptbins,ptmin,ptmax,200,0,2);
 471         fhXEIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 472         fhXEIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 473         outputContainer->Add(fhXEIsolated[p][j][i]);
 474
 475
 476         // XE from UE
 477
 478         fhXEUE[p][j][i]  = new TH2F(Form("h%sXEUE%s%s",                               particle[p].Data(), leading[j].Data(), name[i].Data()),
 479                                     Form("#it{z}_{Jet} of %s: %s of all particles%s", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 480                                     nptbins,ptmin,ptmax,200,0,2);
 481         fhXEUE[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 482         fhXEUE[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 483         outputContainer->Add(fhXEUE[p][j][i]);
 484
 485         fhXEUEIsolated[p][j][i]  = new TH2F(Form("h%sXEUE%sIsolated%s",                                 particle[p].Data(), leading[j].Data(), name[i].Data()),
 486                                             Form("#it{z}_{Jet} of %s: %s of all particles%s, isolated", partTitl[p].Data(), leading[j].Data(), title[i].Data()),
 487                                             nptbins,ptmin,ptmax,200,0,2);
 488         fhXEUEIsolated[p][j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
 489         fhXEUEIsolated[p][j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
 490         outputContainer->Add(fhXEUEIsolated[p][j][i]);
 491       }
 492     }
 493   }
 494
 495   return outputContainer;
 496
 497 }
 498
 499 //____________________________________________
 500 void  AliAnaGeneratorKine::GetPartonsAndJets()
 501 {
 502   // Fill data members with partons,jets and generated pt hard
 503
 504   AliDebug(1,"Start");
 505
 506 //  if( nPrimary > 2 ) fParton2 = fStack->Particle(2) ;
 507 //  if( nPrimary > 3 ) fParton3 = fStack->Particle(3) ;
 508
 509   Float_t p6phi = -1 ;
 510   Float_t p6eta = -10;
 511   Float_t p6pt  =  0 ;
 512   if( fNPrimaries > 6 )
 513   {
 514     p6pt  = fParton6.Pt();
 515     p6eta = fParton6.Eta();
 516     p6phi = fParton6.Phi();
 517     if(p6phi < 0) p6phi +=TMath::TwoPi();
 518   }
 519
 520   Float_t p7phi = -1 ;
 521   Float_t p7eta = -10;
 522   Float_t p7pt  =  0 ;
 523   if( fNPrimaries > 7 )
 524   {
 525     p7pt  = fParton7.Pt();
 526     p7phi = fParton7.Eta();
 527     p7phi = fParton7.Phi();
 528     if(p7phi < 0) p7phi +=TMath::TwoPi();
 529   }
 530
 531   //printf("parton6: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",p6pt,p6eta,p6phi, fParton6PDG);
 532   //printf("parton7: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",p7pt,p7eta,p7phi, fParton7PDG);
 533
 534   // Get the jets, only for pythia
 535   if(!strcmp(GetReader()->GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
 536   {
 537     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetReader()->GetGenEventHeader();
 538
 539     fPtHard = pygeh->GetPtHard();
 540
 541     //printf("pt Hard %2.2f\n",fPtHard);
 542
 543     const Int_t nTriggerJets =  pygeh->NTriggerJets();
 544
 545     Float_t tmpjet[]={0,0,0,0};
 546
 547     // select the closest jet to parton
 548     Float_t jet7R = 100;
 549     Float_t jet6R = 100;
 550
 551     for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
 552     {
 553       pygeh->TriggerJet(ijet, tmpjet);
 554
 555       fLVTmp.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
 556       Float_t jphi = fLVTmp.Phi();
 557       if(jphi < 0) jphi +=TMath::TwoPi();
 558
 559       Double_t radius6 = GetIsolationCut()->Radius(p6eta, p6phi, fLVTmp.Eta() , jphi) ;
 560       Double_t radius7 = GetIsolationCut()->Radius(p7eta, p7phi, fLVTmp.Eta() , jphi) ;
 561
 562       //printf("jet %d: pt %2.2f, eta %2.2f, phi %2.2f, r6 %2.2f, r7 %2.2f\n",ijet,jet.Pt(),jet.Eta(),jphi,radius6, radius7);
 563
 564       if (radius6 < jet6R)
 565       {
 566         jet6R = radius6;
 567         fJet6 = fLVTmp;
 568
 569       }
 570
 571       if (radius7 < jet7R)
 572       {
 573         jet7R = radius7;
 574         fJet7 = fLVTmp;
 575       }
 576
 577     } // jet loop
 578
 579     //printf("jet6: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet6.Pt(),fJet6.Eta(),fJet6.Phi());
 580     //printf("jet7: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet7.Pt(),fJet7.Eta(),fJet6.Phi());
 581
 582   } // pythia header
 583
 584   fhPtHard   ->Fill(fPtHard);
 585   fhPtJet    ->Fill(fJet6.Pt());
 586   fhPtJet    ->Fill(fJet7.Pt());
 587   fhPtParton ->Fill(p6pt);
 588   fhPtParton ->Fill(p7pt);
 589
 590   if( fPtHard > 0 )
 591   {
 592     fhPtPartonPtHard->Fill(fPtHard, p6pt/fPtHard);
 593     fhPtPartonPtHard->Fill(fPtHard, p7pt/fPtHard);
 594     fhPtJetPtHard   ->Fill(fPtHard, fJet6.Pt()/fPtHard);
 595     fhPtJetPtHard   ->Fill(fPtHard, fJet7.Pt()/fPtHard);
 596   }
 597
 598   if( p6pt > 0 ) fhPtJetPtParton ->Fill(fPtHard, fJet6.Pt()/p6pt);
 599   if( p7pt > 0 ) fhPtJetPtParton ->Fill(fPtHard, fJet7.Pt()/p7pt);
 600
 601   AliDebug(1,"End");
 602
 603 }
 604
 605 //_____________________________________________________
 606 void AliAnaGeneratorKine::GetXE(Int_t   indexTrig,
 607                                 Int_t   partType,
 608                                 Bool_t  leading [fgkNIso],
 609                                 Bool_t  isolated[fgkNIso],
 610                                 Int_t   iparton)
 611 {
 612
 613   // Calculate the real XE and the UE XE
 614
 615   AliDebug(1,"Start");
 616
 617   Float_t ptTrig  = fTrigger.Pt();
 618   Float_t phiTrig = fTrigger.Phi();
 619   if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
 620
 621   Int_t  pdg    = 0;
 622   Int_t  status = 0;
 623   Int_t  ipartonAway = 0;
 624   Int_t  charge = 0;
 625   //Loop on primaries, start from position 8, no partons
 626   for(Int_t ipr = 8; ipr < fNPrimaries; ipr ++ )
 627   {
 628
 629     if(ipr==indexTrig) continue;
 630
 631     // Get ESD particle kinematics
 632     if(fStack)
 633     {
 634       TParticle * particle = fStack->Particle(ipr) ;
 635
 636       pdg    = particle->GetPdgCode();
 637       status = particle->GetStatusCode();
 638
 639       // Compare trigger with final state particles
 640       if( status != 1) continue ; // do it here to avoid crashes
 641
 642       particle->Momentum(fLVTmp);
 643
 644       charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
 645     }
 646     else // AOD particle kinematics
 647     {
 648       AliAODMCParticle * particle = (AliAODMCParticle*) fAODMCparticles->At(ipr) ;
 649
 650       pdg    = particle->GetPdgCode();
 651       status = particle->GetStatus();
 652
 653       // Compare trigger with final state particles
 654       if( status != 1) continue ; // do it here to avoid crashes
 655
 656       fLVTmp.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),particle->E());
 657
 658       charge = particle->Charge();
 659     }
 660
 661     // construct xe only with charged
 662     if( charge == 0 ) continue;
 663
 664     Float_t pt  = fLVTmp.Pt();
 665     Float_t phi = fLVTmp.Phi();
 666     if(phi < 0 ) phi += TMath::TwoPi();
 667
 668     if( pt < fMinChargedPt)    continue ;
 669
 670     Bool_t inTPC = GetFiducialCut()->IsInFiducialCut(fLVTmp.Eta(),fLVTmp.Phi(),kCTS) ;
 671
 672     if(!inTPC) continue;
 673
 674     // ---------------------------------------------------
 675     // Get the index of the mother, get from what parton
 676     // ESD
 677     if(fStack)
 678     {
 679       ipartonAway =  fStack->Particle(ipr)->GetFirstMother();
 680       if(ipartonAway < 0)
 681       {
 682         AliDebug(1,"End, no mother index");
 683         return;
 684       }
 685
 686       TParticle * mother = fStack->Particle(ipartonAway);
 687       while (ipartonAway > 7)
 688       {
 689         ipartonAway   = mother->GetFirstMother();
 690         if(ipartonAway < 0) break;
 691         mother = fStack->Particle(ipartonAway);
 692       }
 693     }
 694     else // AOD
 695     {
 696       ipartonAway =  ((AliAODMCParticle*) fAODMCparticles->At(ipr))->GetMother();
 697       if(ipartonAway < 0)
 698       {
 699         AliDebug(1,"End, no mother index");
 700         return;
 701       }
 702
 703       AliAODMCParticle * mother = (AliAODMCParticle*) fAODMCparticles->At(ipartonAway);
 704       while (ipartonAway > 7)
 705       {
 706         ipartonAway   = mother->GetMother();
 707         if(ipartonAway < 0) break;
 708         mother = (AliAODMCParticle*) fAODMCparticles->At(ipartonAway);
 709       }
 710     }
 711
 712     //-----------------------------------------
 713     // Get XE of particles belonging to the jet
 714     // on the opposite side of the trigger
 715
 716     Float_t xe = -pt/ptTrig*TMath::Cos(phi-phiTrig);
 717
 718     if((ipartonAway==6 || ipartonAway==7) && iparton!=ipartonAway)
 719     {
 720       for( Int_t i = 0; i < fgkNIso; i++ )
 721       {
 722         fhXE[partType][leading[i]][i]  ->Fill(ptTrig,xe);
 723
 724         if(isolated[i])
 725         {
 726           fhXEIsolated[partType][leading[i]][i]  ->Fill(ptTrig,xe);
 727         }
 728       } // conditions loop
 729     } // Away side
 730
 731     //----------------------------------------------------------
 732     // Get the XE from particles not attached to any of the jets
 733     if(ipartonAway!=6 && ipartonAway!=7)
 734     {
 735       for( Int_t i = 0; i < fgkNIso; i++ )
 736       {
 737         fhXEUE[partType][leading[i]][i]  ->Fill(ptTrig,xe);
 738
 739         if(isolated[i])
 740         {
 741           fhXEUEIsolated[partType][leading[i]][i]  ->Fill(ptTrig,xe);
 742         }
 743       } // conditions loop
 744     } // Away side
 745
 746   } // primary loop
 747
 748   AliDebug(1,"End");
 749
 750 }
 751
 752 //________________________________________
 753 void AliAnaGeneratorKine::InitParameters()
 754 {
 755
 756   //Initialize the parameters of the analysis.
 757   AddToHistogramsName("AnaGenKine_");
 758
 759   fTriggerDetector = kEMCAL;
 760
 761   fMinChargedPt    = 0.2;
 762   fMinNeutralPt    = 0.3;
 763
 764 }
 765
 766 //_____________________________________________________________________
 767 void  AliAnaGeneratorKine::IsLeadingAndIsolated(Int_t indexTrig,
 768                                                 Int_t partType,
 769                                                 Bool_t leading[fgkNIso],
 770                                                 Bool_t isolated[fgkNIso])
 771 {
 772   // Check if the trigger is the leading particle and if it is isolated
 773   // In case of neutral particles check all neutral or neutral in EMCAL acceptance
 774
 775   AliDebug(1,"Start");
 776
 777   Float_t ptMaxCharged       = 0; // all charged
 778   Float_t ptMaxNeutral       = 0; // all neutral
 779   Float_t ptMaxNeutEMCAL     = 0; // for neutral, select them in EMCAL acceptance
 780   Float_t ptMaxNeutPhot      = 0; // for neutral, take only photons
 781   Float_t ptMaxNeutEMCALPhot = 0; // for neutral, take only photons in EMCAL acceptance
 782
 783   leading[0] = 0;
 784   leading[1] = 0;
 785   leading[2] = 0;
 786   leading[3] = 0;
 787
 788   isolated[0] = 0;
 789   isolated[1] = 0;
 790   isolated[2] = 0;
 791   isolated[3] = 0;
 792
 793   Float_t ptTrig  = fTrigger.Pt();
 794   Float_t etaTrig = fTrigger.Eta();
 795   Float_t phiTrig = fTrigger.Phi();
 796   if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
 797
 798   // Minimum track or cluster energy
 799
 800   //Isolation cuts
 801   Float_t ptThresIC    = GetIsolationCut()->GetPtThreshold();
 802   Float_t sumThresIC   = GetIsolationCut()->GetPtThreshold();
 803   Float_t rThresIC     = GetIsolationCut()->GetConeSize();
 804   Float_t isoMethod    = GetIsolationCut()->GetICMethod();
 805
 806   //Counters
 807   Int_t   nICTrack     = 0;
 808   Int_t   nICNeutral   = 0;
 809   Int_t   nICNeutEMCAL = 0;
 810   Int_t   nICNeutPhot  = 0;
 811   Int_t   nICNeutEMCALPhot = 0;
 812
 813   // Sum of pT
 814   Float_t sumNePt         = 0;
 815   Float_t sumNePtPhot     = 0;
 816   Float_t sumNePtEMC      = 0;
 817   Float_t sumNePtEMCPhot  = 0;
 818   Float_t sumChPt         = 0;
 819
 820   //Loop on primaries, start from position 8, no partons
 821
 822   Int_t imother = -1;
 823   Int_t pdg     =  0;
 824   Int_t status  =  0;
 825   Int_t charge  =  0;
 826   for(Int_t ipr = 8; ipr < fNPrimaries; ipr ++ )
 827   {
 828     if(ipr == indexTrig) continue;
 829
 830     if(fStack) // ESD
 831     {
 832       TParticle * particle = fStack->Particle(ipr) ;
 833
 834       imother = particle->GetFirstMother();
 835       pdg     = particle->GetPdgCode();
 836       status  = particle->GetStatusCode();
 837
 838       // Compare trigger with final state particles
 839       if( status != 1) continue ; // do it here to avoid crashes
 840
 841       charge  = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
 842       particle->Momentum(fLVTmp);
 843     }
 844     else // AOD
 845     {
 846       AliAODMCParticle * particle = (AliAODMCParticle*) fAODMCparticles->At(ipr) ;
 847
 848       imother = particle->GetMother();
 849       pdg     = particle->GetPdgCode();
 850       status  = particle->GetStatus();
 851
 852       // Compare trigger with final state particles
 853       if( status != 1) continue ; // do it here to avoid crashes
 854
 855       charge  = particle->Charge();
 856       fLVTmp.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),particle->E());
 857     }
 858
 859     // Do not consider the photon decays from pi0 and eta
 860     //printf("Leading ipr %d - mother %d - iTrig\n",ipr, imother,indexTrig);
 861     if( imother == indexTrig)  continue ;
 862
 863     Float_t pt  = fLVTmp.Pt();
 864     Float_t eta = fLVTmp.Eta();
 865     Float_t phi = fLVTmp.Phi();
 866     if(phi < 0 ) phi += TMath::TwoPi();
 867
 868     // Select all particles in at least the TPC acceptance
 869     Bool_t inTPC = GetFiducialCut()->IsInFiducialCut(eta,phi,kCTS) ;
 870     if(!inTPC) continue;
 871
 872     //Isolation
 873     Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
 874
 875     if(charge==0)
 876     {
 877       if(pt < fMinNeutralPt)  continue ;
 878
 879       if( ptMaxNeutral < pt ) ptMaxNeutral = pt;
 880
 881       if( radius < rThresIC )
 882       {
 883         if( pt > ptThresIC ) nICNeutral++ ;
 884         sumNePt+= pt;
 885       }
 886
 887       Bool_t phPDG = kFALSE;
 888       if(pdg==22 || pdg==111) phPDG = kTRUE;
 889
 890 //      if(pt > ptTrig) printf(" --- pdg %d, phPDG %d pT %2.2f, pTtrig %2.2f, eta %2.2f, phi %2.2f\n",pdg,phPDG,pt,ptTrig,eta, phi*TMath::RadToDeg());
 891       if(phPDG)
 892       {
 893         if( ptMaxNeutPhot < pt) ptMaxNeutPhot = pt;
 894
 895         if( radius < rThresIC )
 896         {
 897           if(pt > ptThresIC) nICNeutPhot++ ;
 898           sumNePtPhot += pt;
 899         }
 900       }
 901
 902       //Calorimeter acceptance
 903       Bool_t inCalo = GetFiducialCut()->IsInFiducialCut(eta,phi,GetCalorimeter()) ;
 904       if(!inCalo) continue;
 905
 906       if( ptMaxNeutEMCAL < pt ) ptMaxNeutEMCAL = pt;
 907       if( radius < rThresIC )
 908       {
 909         if( pt > ptThresIC ) nICNeutEMCAL++ ;
 910         sumNePtEMC += pt;
 911       }
 912
 913       if(phPDG)
 914       {
 915         if( ptMaxNeutEMCALPhot < pt ) ptMaxNeutEMCALPhot = pt;
 916         if(  radius < rThresIC )
 917         {
 918           if (pt > ptThresIC) nICNeutEMCALPhot++ ;
 919           sumNePtEMCPhot += pt;
 920         }
 921       }
 922     }
 923     else
 924     {
 925       if( pt < fMinChargedPt)  continue ;
 926
 927       if( ptMaxCharged < pt )   ptMaxCharged   = pt;
 928
 929       if( radius < rThresIC )
 930       {
 931 //        printf("UE track? pTtrig %2.2f, pt %2.2f, etaTrig %2.2f,  eta %2.2f, phiTrig %2.2f,  phi %2.2f, radius %2.2f\n",
 932 //               ptTrig, pt,etaTrig, eta, phiTrig*TMath::RadToDeg(), phi*TMath::RadToDeg(),radius);
 933         if( pt > ptThresIC ) nICTrack++ ;
 934         sumChPt += pt;
 935       }
 936     }
 937
 938   } // particle loop
 939
 940
 941   // Leding decision
 942   if(ptTrig > ptMaxCharged)
 943   {
 944 //    printf("pt charged %2.2f, pt neutral %2.2f, pt neutral emcal %2.2f, pt photon %2.2f, pt photon emcal %2.2f\n",
 945 //           ptMaxCharged, ptMaxNeutral, ptMaxNeutEMCAL,ptMaxNeutPhot, ptMaxNeutEMCALPhot);
 946     if(ptTrig > ptMaxNeutral      ) leading[0] = kTRUE ;
 947     if(ptTrig > ptMaxNeutEMCAL    ) leading[1] = kTRUE ;
 948     if(ptTrig > ptMaxNeutPhot     ) leading[2] = kTRUE ;
 949     if(ptTrig > ptMaxNeutEMCALPhot) leading[3] = kTRUE ;
 950   }
 951
 952 //  printf("N in cone over threshold: tracks  %d, neutral %d, neutral emcal %d, photon %d, photon emcal %d\n",
 953 //         nICTrack, nICNeutral ,nICNeutEMCAL,nICNeutPhot, nICNeutEMCALPhot);
 954
 955   //------------------
 956   // Isolation decision
 957   if( isoMethod == AliIsolationCut::kPtThresIC )
 958   {
 959     if( nICTrack == 0 )
 960     {
 961       if(nICNeutral       == 0 ) isolated[0] = kTRUE ;
 962       if(nICNeutEMCAL     == 0 ) isolated[1] = kTRUE ;
 963       if(nICNeutPhot      == 0 ) isolated[2] = kTRUE ;
 964       if(nICNeutEMCALPhot == 0 ) isolated[3] = kTRUE ;
 965     }
 966   }
 967   else if( isoMethod == AliIsolationCut::kSumPtIC )
 968   {
 969     if(sumChPt + sumNePt        < sumThresIC ) isolated[0] = kTRUE ;
 970     if(sumChPt + sumNePtEMC     < sumThresIC ) isolated[1] = kTRUE ;
 971     if(sumChPt + sumNePtPhot    < sumThresIC ) isolated[2] = kTRUE ;
 972     if(sumChPt + sumNePtEMCPhot < sumThresIC ) isolated[3] = kTRUE ;
 973   }
 974
 975
 976   //----------------------------------------------------
 977   // Fill histograms if conditions apply for all 4 cases
 978   for( Int_t i = 0; i < fgkNIso; i++ )
 979   {
 980     if(leading[i])
 981     {
 982       fhPtLeading[partType][i]->Fill(ptTrig);
 983
 984       if     (i == 0) fhPtLeadingSumPt[partType][i]->Fill(ptTrig, sumChPt + sumNePt);
 985       else if(i == 1) fhPtLeadingSumPt[partType][i]->Fill(ptTrig, sumChPt + sumNePtEMC);
 986       else if(i == 2) fhPtLeadingSumPt[partType][i]->Fill(ptTrig, sumChPt + sumNePtPhot);
 987       else if(i == 3) fhPtLeadingSumPt[partType][i]->Fill(ptTrig, sumChPt + sumNePtEMCPhot);
 988
 989       if(isolated[i]) fhPtLeadingIsolated[partType][i]->Fill(ptTrig);
 990     }
 991   } // conditions loop
 992
 993   AliDebug(1,"End");
 994
 995 }
 996
 997 //_____________________________________________________
 998 void  AliAnaGeneratorKine::MakeAnalysisFillHistograms()
 999 {
1000   // Particle-Parton/Jet/Hadron Correlation Analysis, main method
1001
1002   AliDebug(1,"Start");
1003
1004   fParton6.SetPxPyPzE(0,0,0,0);
1005   fParton7.SetPxPyPzE(0,0,0,0);
1006   fParton6PDG = 0;
1007   fParton7PDG = 0;
1008
1009   //
1010   // Get the ESD MC particles container
1011   fStack = 0;
1012   if( GetReader()->ReadStack() )
1013   {
1014     fStack = GetMCStack();
1015     if( !fStack )
1016     {
1017       AliFatal("Stack not available, is the MC handler called? STOP");
1018       return;
1019     }
1020
1021     fNPrimaries = fStack->GetNprimary(); // GetNtrack();
1022
1023     if(fNPrimaries > 6)
1024     {
1025       (fStack->Particle(6))->Momentum(fParton6) ;
1026       fParton6PDG =  (fStack->Particle(6))->GetPdgCode();
1027     }
1028
1029     if(fNPrimaries > 7)
1030     {
1031       (fStack->Particle(7))->Momentum(fParton7) ;
1032       fParton7PDG =  (fStack->Particle(7))->GetPdgCode();
1033     }
1034   }
1035
1036   // Get the AOD MC particles container
1037   fAODMCparticles = 0;
1038   if( GetReader()->ReadAODMCParticles() )
1039   {
1040     fAODMCparticles = GetReader()->GetAODMCParticles();
1041     if( !fAODMCparticles )
1042     {
1043       AliFatal("Standard MCParticles not available!");
1044       return;
1045     }
1046
1047     fNPrimaries = fAODMCparticles->GetEntriesFast();
1048     AliAODMCParticle * primAOD = 0;
1049     if(fNPrimaries > 6)
1050     {
1051       primAOD = (AliAODMCParticle *) fAODMCparticles->At(6);
1052       fParton6.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
1053
1054       fParton6PDG =  primAOD->GetPdgCode();
1055     }
1056
1057     if(fNPrimaries > 7)
1058     {
1059       primAOD = (AliAODMCParticle *) fAODMCparticles->At(7);
1060       fParton7.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
1061
1062       fParton7PDG =  primAOD->GetPdgCode();
1063     }
1064   }
1065   //
1066
1067   GetPartonsAndJets();
1068
1069   // Main particle loop
1070   Int_t   pdgTrig    = 0;
1071   Int_t   statusTrig = 0;
1072   Int_t   imother    = 0;
1073   Float_t ptTrig     = 0;
1074   Int_t   momStatus  = 0;
1075   Int_t   momPdg     = 0;
1076   Int_t   pdg0       = 0;
1077   Int_t   pdg1       = 0;
1078   Int_t   id0        = 0;
1079   Int_t   id1        = 0;
1080   Int_t   nDaughters = 0;
1081
1082   for(Int_t ipr = 0; ipr < fNPrimaries; ipr ++ )
1083   {
1084     if(fStack) // ESD
1085     {
1086       TParticle * particle = fStack->Particle(ipr) ;
1087
1088       pdgTrig    = particle->GetPdgCode();
1089       statusTrig = particle->GetStatusCode();
1090       imother    = particle->GetFirstMother();
1091       ptTrig     = particle->Pt();
1092       nDaughters = particle->GetNDaughters();
1093       id0        = particle->GetDaughter(0);
1094       id1        = particle->GetDaughter(1);
1095       // Recover the kinematics:
1096       particle->Momentum(fTrigger);
1097     }
1098     else // AOD
1099     {
1100       AliAODMCParticle* particle = (AliAODMCParticle*) fAODMCparticles->At(ipr) ;
1101
1102       pdgTrig    = particle->GetPdgCode();
1103       statusTrig = particle->GetStatus();
1104       imother    = particle->GetMother();
1105       nDaughters = particle->GetNDaughters();
1106       id0        = particle->GetDaughter(0);
1107       id1        = particle->GetDaughter(1);
1108       // Recover the kinematics:
1109       fTrigger.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),particle->E());
1110     }
1111
1112     // Select final state photons or pi0s or eta's
1113
1114     if( pdgTrig == 22 && statusTrig != 1 ) continue ;
1115
1116     if( pdgTrig != 111 && pdgTrig != 22 && pdgTrig !=221 ) continue ;
1117
1118     // Acceptance and kinematical cuts
1119
1120     Float_t ptTrig  = fTrigger.Pt();
1121
1122     if( ptTrig < GetMinPt() ) continue ;
1123
1124     Bool_t in = GetFiducialCutForTrigger()->IsInFiducialCut(fTrigger.Eta(),fTrigger.Phi(),fTriggerDetector) ;
1125     if(! in )  continue ;
1126
1127     // Identify the particle to fill appropriate histogram
1128     Int_t partType = -1;
1129
1130     if     (pdgTrig==22 )
1131     {
1132       if(imother > 0 )
1133       {
1134         if(fStack) // ESD
1135         {
1136           momStatus = (fStack->Particle(imother))->GetStatusCode();
1137           momPdg    = (fStack->Particle(imother))->GetPdgCode();
1138         }
1139         else // AOD
1140         {
1141           momStatus = ((AliAODMCParticle*) fAODMCparticles->At(imother))->GetStatus();
1142           momPdg    = ((AliAODMCParticle*) fAODMCparticles->At(imother))->GetPdgCode();
1143         }
1144
1145         if     (imother < 8 && statusTrig == 1)  partType = kmcPrimPhoton ;
1146         else if(momPdg == 111 ) partType = kmcPrimPi0Decay   ;
1147         else if(momPdg == 221 ) partType = kmcPrimEtaDecay   ;
1148         else if(momStatus > 0 ) partType = kmcPrimOtherDecay ;
1149       }
1150     }
1151     else if( (pdgTrig==111 || pdgTrig==221) && nDaughters == 2 )
1152     {
1153       if(fStack) // ESD
1154       {
1155         pdg0 = fStack->Particle(id0)->GetPdgCode();
1156         pdg1 = fStack->Particle(id1)->GetPdgCode();
1157       }
1158       else // AOD
1159       {
1160         pdg0 = ((AliAODMCParticle*) fAODMCparticles->At(id0))->GetPdgCode();
1161         pdg1 = ((AliAODMCParticle*) fAODMCparticles->At(id1))->GetPdgCode();
1162       }
1163
1164       if( pdg0 == 22 && pdg1== 22 )
1165       {
1166         if     ( pdgTrig==111 ) partType = kmcPrimPi0;
1167         else if( pdgTrig==221 ) partType = kmcPrimEta;
1168       }
1169     }
1170
1171     if(partType < 0 ) continue ;
1172
1173     AliDebug(1,Form("Select trigger particle %d: pdg %d, type %d, status %d, mother index %d, pT %2.2f, eta %2.2f, phi %2.2f",
1174                     ipr, pdgTrig, partType, statusTrig, imother, ptTrig, fTrigger.Eta(), fTrigger.Phi()*TMath::RadToDeg()));
1175
1176     //    if(pdgTrig==111)
1177     //    {
1178     //      printf("\t pi0 daughters %d, %d\n", particle->GetDaughter(0), particle->GetDaughter(1));
1179     //    }
1180
1181     // Fill histograms do analysis
1182
1183     fhPt[partType]->Fill(ptTrig);
1184
1185     // Check if it is leading
1186     Bool_t leading [fgkNIso] ;
1187     Bool_t isolated[fgkNIso] ;
1188
1189     IsLeadingAndIsolated(ipr, partType, leading, isolated);
1190
1191     Int_t iparton = -1;
1192     Int_t ok = CorrelateWithPartonOrJet(ipr, partType, leading, isolated, iparton);
1193     if(!ok) continue;
1194
1195     GetXE(ipr,partType,leading,isolated,iparton) ;
1196
1197   }
1198
1199   AliDebug(1,"End fill histograms");
1200
1201 }
1202
1203 //_________________________________________________________
1204 void AliAnaGeneratorKine::SetTriggerDetector(TString & det)
1205 {
1206   // Set the detrimeter for the analysis
1207
1208   fTriggerDetectorString = det;
1209
1210   if     (det=="EMCAL") fTriggerDetector = kEMCAL;
1211   else if(det=="PHOS" ) fTriggerDetector = kPHOS;
1212   else if(det=="CTS")   fTriggerDetector = kCTS;
1213   else if(det=="DCAL")  fTriggerDetector = kDCAL;
1214   else if(det.Contains("DCAL") && det.Contains("PHOS")) fTriggerDetector = kDCALPHOS;
1215   else AliFatal(Form("Detector < %s > not known!", det.Data()));
1216
1217 }
1218
1219 //_____________________________________________________
1220 void AliAnaGeneratorKine::SetTriggerDetector(Int_t det)
1221 {
1222   // Set the detrimeter for the analysis
1223
1224   fTriggerDetector = det;
1225
1226   if     (det==kEMCAL)    fTriggerDetectorString = "EMCAL";
1227   else if(det==kPHOS )    fTriggerDetectorString = "PHOS";
1228   else if(det==kCTS)      fTriggerDetectorString = "CTS";
1229   else if(det==kDCAL)     fTriggerDetectorString = "DCAL";
1230   else if(det==kDCALPHOS) fTriggerDetectorString = "DCAL_PHOS";
1231   else AliFatal(Form("Detector < %d > not known!", det));
1232
1233 }
1234