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