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