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