PWG4/PartCorrDep/AliAnaParticlePartonCorrelation.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 /* $Id:  $ */
  16
  17 //_________________________________________________________________________
  18 // Class for the analysis of particle-parton correlations
  19 // Particle (for example direct gamma) must be found in a previous analysis
  20 // -- Author: Gustavo Conesa (LNF-INFN)
  21 //////////////////////////////////////////////////////////////////////////////
  22
  23 // --- ROOT system ---
  24 //#include "Riostream.h"
  25 #include "TH2F.h"
  26 #include "TParticle.h"
  27
  28 //---- ANALYSIS system ----
  29 #include "AliAnaParticlePartonCorrelation.h"
  30 #include "AliStack.h"
  31 #include "AliAODPWG4ParticleCorrelation.h"
  32
  33   ClassImp(AliAnaParticlePartonCorrelation)
  34
  35
  36 //____________________________________________________________________________
  37   AliAnaParticlePartonCorrelation::AliAnaParticlePartonCorrelation() :
  38     AliAnaPartCorrBaseClass(),
  39     fhDeltaEtaNearParton(0), fhDeltaPhiNearParton(0),
  40     fhDeltaPtNearParton(0), fhPtRatNearParton(0),
  41     fhDeltaEtaAwayParton(0), fhDeltaPhiAwayParton(0),
  42     fhDeltaPtAwayParton(0), fhPtRatAwayParton(0)
  43 {
  44   //Default Ctor
  45
  46   //Initialize parameters
  47   InitParameters();
  48 }
  49
  50 //____________________________________________________________________________
  51 AliAnaParticlePartonCorrelation::AliAnaParticlePartonCorrelation(const AliAnaParticlePartonCorrelation & g) :
  52   AliAnaPartCorrBaseClass(g),
  53   fhDeltaEtaNearParton(g.fhDeltaEtaNearParton), fhDeltaPhiNearParton(g.fhDeltaPhiNearParton),
  54   fhDeltaPtNearParton(g.fhDeltaPtNearParton), fhPtRatNearParton(g.fhPtRatNearParton),
  55   fhDeltaEtaAwayParton(g.fhDeltaEtaAwayParton), fhDeltaPhiAwayParton(g.fhDeltaPhiAwayParton),
  56   fhDeltaPtAwayParton(g.fhDeltaPtAwayParton), fhPtRatAwayParton(g.fhPtRatAwayParton)
  57 {
  58   // cpy ctor
  59
  60 }
  61
  62 //_________________________________________________________________________
  63 AliAnaParticlePartonCorrelation & AliAnaParticlePartonCorrelation::operator = (const AliAnaParticlePartonCorrelation & source)
  64 {
  65   // assignment operator
  66
  67   if(this == &source)return *this;
  68   ((AliAnaPartCorrBaseClass *)this)->operator=(source);
  69
  70   fhDeltaEtaAwayParton = source.fhDeltaEtaAwayParton;
  71   fhDeltaPhiAwayParton = source.fhDeltaPhiAwayParton;
  72   fhDeltaPtAwayParton = source.fhDeltaPtAwayParton;
  73   fhPtRatAwayParton = source.fhPtRatAwayParton;
  74   fhDeltaEtaNearParton = source.fhDeltaEtaNearParton;
  75   fhDeltaPhiNearParton = source.fhDeltaPhiNearParton;
  76   fhDeltaPtNearParton = source.fhDeltaPtNearParton;
  77   fhPtRatNearParton = source.fhPtRatNearParton;
  78
  79   return *this;
  80
  81 }
  82
  83
  84 //________________________________________________________________________
  85 TList *  AliAnaParticlePartonCorrelation::GetCreateOutputObjects()
  86 {
  87   // Create histograms to be saved in output file
  88
  89   TList * outputContainer = new TList() ;
  90   outputContainer->SetName("ParticlePartonHistos") ;
  91
  92   fhDeltaPhiNearParton  = new TH2F
  93     ("DeltaPhiNearParton","#phi_{particle} - #phi_{parton} vs p_{T particle}",
  94      200,0,120,200,0,6.4);
  95   fhDeltaPhiNearParton->SetYTitle("#Delta #phi");
  96   fhDeltaPhiNearParton->SetXTitle("p_{T particle} (GeV/c)");
  97   outputContainer->Add(fhDeltaPhiNearParton);
  98
  99   fhDeltaEtaNearParton  = new TH2F
 100     ("DeltaEtaNearParton","#eta_{particle} - #eta_{parton} vs p_{T particle}",
 101      200,0,120,200,-2,2);
 102   fhDeltaEtaNearParton->SetYTitle("#Delta #eta");
 103   fhDeltaEtaNearParton->SetXTitle("p_{T particle} (GeV/c)");
 104   outputContainer->Add(fhDeltaEtaNearParton);
 105
 106   fhDeltaPtNearParton  = new TH2F
 107     ("DeltaPtNearParton","#p_{T particle} - #p_{T parton} vs p_{T particle}",
 108      200,0,120,100,-10,10);
 109   fhDeltaPtNearParton->SetYTitle("#Delta #p_{T}");
 110   fhDeltaPtNearParton->SetXTitle("p_{T particle} (GeV/c)");
 111   outputContainer->Add(fhDeltaPtNearParton);
 112
 113   fhPtRatNearParton  = new TH2F
 114     ("PtRatNearParton","#p_{T parton} / #p_{T particle} vs p_{T particle}",
 115      200,0,120,200,0,5);
 116   fhPtRatNearParton->SetYTitle("ratio");
 117   fhPtRatNearParton->SetXTitle("p_{T particle} (GeV/c)");
 118   outputContainer->Add(fhPtRatNearParton);
 119
 120   fhDeltaPhiAwayParton  = new TH2F
 121     ("DeltaPhiAwayParton","#phi_{particle} - #phi_{parton} vs p_{T particle}",
 122      200,0,120,200,0,6.4);
 123   fhDeltaPhiAwayParton->SetYTitle("#Delta #phi");
 124   fhDeltaPhiAwayParton->SetXTitle("p_{T particle} (GeV/c)");
 125   outputContainer->Add(fhDeltaPhiAwayParton);
 126
 127   fhDeltaEtaAwayParton  = new TH2F
 128     ("DeltaEtaAwayParton","#eta_{particle} - #eta_{parton} vs p_{T particle}",
 129      200,0,120,200,-2,2);
 130   fhDeltaEtaAwayParton->SetYTitle("#Delta #eta");
 131   fhDeltaEtaAwayParton->SetXTitle("p_{T particle} (GeV/c)");
 132   outputContainer->Add(fhDeltaEtaAwayParton);
 133
 134   fhDeltaPtAwayParton  = new TH2F
 135     ("DeltaPtAwayParton","#p_{T particle} - #p_{T parton} vs p_{T particle}",
 136      200,0,120,100,-10,10);
 137   fhDeltaPtAwayParton->SetYTitle("#Delta #p_{T}");
 138   fhDeltaPtAwayParton->SetXTitle("p_{T particle} (GeV/c)");
 139   outputContainer->Add(fhDeltaPtAwayParton);
 140
 141   fhPtRatAwayParton  = new TH2F
 142     ("PtRatAwayParton","#p_{T parton} / #p_{T particle} vs p_{T particle}",
 143      200,0,120,200,0,5);
 144   fhPtRatAwayParton->SetYTitle("ratio");
 145   fhPtRatAwayParton->SetXTitle("p_{T particle} (GeV/c)");
 146   outputContainer->Add(fhPtRatAwayParton);
 147
 148   return outputContainer;
 149 }
 150
 151 //____________________________________________________________________________
 152 void AliAnaParticlePartonCorrelation::InitParameters()
 153 {
 154
 155   //Initialize the parameters of the analysis.
 156   SetInputAODName("photons");
 157
 158 }
 159
 160 //__________________________________________________________________
 161 void AliAnaParticlePartonCorrelation::Print(const Option_t * opt) const
 162 {
 163
 164   //Print some relevant parameters set for the analysis
 165   if(! opt)
 166     return;
 167
 168 }
 169
 170 //__________________________________________________________________
 171 void  AliAnaParticlePartonCorrelation::MakeAnalysisFillAOD()
 172 {
 173   //Particle-Parton Correlation Analysis, create AODs
 174   //Add partons to the reference list of the trigger particle
 175   //Partons are considered those in the first eight possitions in the stack
 176   //being 0, and 1 the 2 protons, and 6 and 7 the outgoing final partons.
 177   if(!GetInputAODBranch()){
 178     printf("ParticlePartonCorrelation::FillAOD: No input particles in AOD with name branch < %s > \n",GetInputAODName().Data());
 179     abort();
 180   }
 181   if(GetDebug() > 1){
 182     printf("Begin parton correlation analysis, fill AODs \n");
 183     printf("In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
 184   }
 185
 186   //Loop on stored AOD particles
 187   Int_t naod = GetInputAODBranch()->GetEntriesFast();
 188   for(Int_t iaod = 0; iaod < naod ; iaod++){
 189     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 190
 191     AliStack * stack =  GetMCStack() ;
 192     if(!stack){
 193       printf("No Stack available, STOP");
 194       abort();
 195     }
 196     if(stack->GetNtrack() < 8) {
 197       printf("*** small number of particles, not a PYTHIA simulation? ***:  n tracks %d \n", stack->GetNprimary());
 198       continue ;
 199     }
 200
 201     //Fill AOD reference only with partons
 202     TParticle * parton = new TParticle ;
 203     Bool_t first = kTRUE;
 204     for(Int_t ipr = 0;ipr < 8; ipr ++ ){
 205       parton = stack->Particle(ipr) ;
 206
 207       if(first) {
 208         new (particle->GetRefTracks()) TRefArray(TProcessID::GetProcessWithUID(parton));
 209         first = kFALSE;
 210       }
 211
 212       particle->AddTrack(parton);
 213     }
 214
 215   }//Aod branch loop
 216
 217   if(GetDebug() > 1) printf("End parton correlation analysis, fill AODs \n");
 218 }
 219
 220 //__________________________________________________________________
 221 void  AliAnaParticlePartonCorrelation::MakeAnalysisFillHistograms()
 222 {
 223   //Particle-Parton Correlation Analysis, fill histograms
 224   if(!GetInputAODBranch()){
 225     printf("ParticlePartonCorrelation::FillHistos: No input particles in AOD with name branch < %s > \n",GetInputAODName().Data());
 226     abort();
 227   }
 228   if(GetDebug() > 1){
 229     printf("Begin parton correlation analysis, fill histograms \n");
 230     printf("In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
 231   }
 232
 233   AliStack * stack =  GetMCStack() ;
 234   if(!stack) {
 235     printf("ParticlePartonCorrelation::FillHistos - No Stack available, STOP");
 236     abort();
 237   }
 238
 239   //Loop on stored AOD particles
 240   Int_t naod = GetInputAODBranch()->GetEntriesFast();
 241   TParticle *  mom =new TParticle ;
 242
 243   for(Int_t iaod = 0; iaod < naod ; iaod++){
 244     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 245
 246     Float_t ptTrigg  = particle->Pt();
 247     Float_t phiTrigg = particle->Phi();
 248     Float_t etaTrigg = particle->Eta();
 249     Int_t imom = particle->GetLabel();
 250     Int_t iparent  = 2000;
 251     Int_t iawayparent = -1;
 252
 253     if(!(particle->GetRefTracks()) || (particle->GetRefTracks())->GetEntriesFast() < 7) {
 254       printf("ParticlePartonCorrelation::FillHistos - Reference list with partons not filled, STOP analysis");
 255       abort();
 256     }
 257     //Check and get indeces of mother and parton
 258     if(imom < 8 ) iparent = imom ;   //mother is already a parton
 259     else if (imom <  stack->GetNtrack()) {
 260       mom =  stack->Particle(imom);
 261       iparent=mom->GetFirstMother();
 262       //cout<<" iparent "<<iparent<<endl;
 263       while(iparent > 7 ){
 264         mom = stack->Particle(iparent);
 265         imom = iparent ; //Mother label is of the inmediate parton daughter
 266         iparent = mom->GetFirstMother();
 267         //cout<<" while iparent "<<iparent<<endl;
 268       }
 269     }
 270
 271     if(GetDebug() > 1) printf("N reference partons %d; labels:  mother %d, parent %d \n", (particle->GetRefTracks())->GetEntriesFast(), imom, iparent);
 272
 273
 274     if(iparent < 0 || iparent > 8) {
 275       if(GetDebug() > 0 ) printf("Failed to find appropriate parton, index %d", iparent);
 276       continue ;
 277     }
 278
 279     //Near parton is the parton that fragmented and created the mother
 280     TParticle * nearParton = (TParticle*) (particle->GetRefTracks())->At(iparent);
 281     Float_t  ptNearParton    = nearParton->Pt();
 282     Float_t  phiNearParton   = nearParton->Phi() ;
 283     Float_t  etaNearParton   = nearParton->Eta() ;
 284
 285     fhDeltaEtaNearParton->Fill(ptTrigg,etaTrigg-etaNearParton);
 286     fhDeltaPhiNearParton->Fill(ptTrigg,phiTrigg-phiNearParton);
 287     fhDeltaPtNearParton->Fill(ptTrigg,ptTrigg-ptNearParton);
 288     fhPtRatNearParton->Fill(ptTrigg,ptNearParton/ptTrigg);
 289
 290     if(iparent == 7) iawayparent =6;
 291     else if(iparent == 6) iawayparent =7;
 292     else{
 293       printf("Parent parton is not final state, skip \n");
 294       continue;
 295     }
 296
 297     //Away parton is the other final parton.
 298     TParticle * awayParton = (TParticle*) (particle->GetRefTracks())->At(iawayparent);
 299     Float_t  ptAwayParton    = awayParton->Pt();
 300     Float_t  phiAwayParton   = awayParton->Phi() ;
 301     Float_t  etaAwayParton   = awayParton->Eta() ;
 302     fhDeltaEtaAwayParton->Fill(ptTrigg,etaTrigg-etaAwayParton);
 303     fhDeltaPhiAwayParton->Fill(ptTrigg,phiTrigg-phiAwayParton);
 304     fhDeltaPtAwayParton->Fill(ptTrigg,ptTrigg-ptAwayParton);
 305     fhPtRatAwayParton->Fill(ptTrigg,ptAwayParton/ptTrigg);
 306
 307
 308   }
 309
 310   if(GetDebug() > 1) printf("End parton correlation analysis, fill histograms \n");
 311
 312 }