PWG4/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 "AliLog.h"
  31 #include "AliStack.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   AliDebug(1,"Init parton histograms");
  90
  91   TList * outputContainer = new TList() ;
  92   outputContainer->SetName("ParticlePartonHistos") ;
  93
  94   fhDeltaPhiNearParton  = new TH2F
  95     ("DeltaPhiNearParton","#phi_{particle} - #phi_{parton} vs p_{T particle}",
  96      200,0,120,200,0,6.4);
  97   fhDeltaPhiNearParton->SetYTitle("#Delta #phi");
  98   fhDeltaPhiNearParton->SetXTitle("p_{T particle} (GeV/c)");
  99   outputContainer->Add(fhDeltaPhiNearParton);
 100
 101   fhDeltaEtaNearParton  = new TH2F
 102     ("DeltaEtaNearParton","#eta_{particle} - #eta_{parton} vs p_{T particle}",
 103      200,0,120,200,-2,2);
 104   fhDeltaEtaNearParton->SetYTitle("#Delta #eta");
 105   fhDeltaEtaNearParton->SetXTitle("p_{T particle} (GeV/c)");
 106   outputContainer->Add(fhDeltaEtaNearParton);
 107
 108   fhDeltaPtNearParton  = new TH2F
 109     ("DeltaPtNearParton","#p_{T particle} - #p_{T parton} vs p_{T particle}",
 110      200,0,120,100,-10,10);
 111   fhDeltaPtNearParton->SetYTitle("#Delta #p_{T}");
 112   fhDeltaPtNearParton->SetXTitle("p_{T particle} (GeV/c)");
 113   outputContainer->Add(fhDeltaPtNearParton);
 114
 115   fhPtRatNearParton  = new TH2F
 116     ("PtRatNearParton","#p_{T parton} / #p_{T particle} vs p_{T particle}",
 117      200,0,120,200,0,5);
 118   fhPtRatNearParton->SetYTitle("ratio");
 119   fhPtRatNearParton->SetXTitle("p_{T particle} (GeV/c)");
 120   outputContainer->Add(fhPtRatNearParton);
 121
 122   fhDeltaPhiAwayParton  = new TH2F
 123     ("DeltaPhiAwayParton","#phi_{particle} - #phi_{parton} vs p_{T particle}",
 124      200,0,120,200,0,6.4);
 125   fhDeltaPhiAwayParton->SetYTitle("#Delta #phi");
 126   fhDeltaPhiAwayParton->SetXTitle("p_{T particle} (GeV/c)");
 127   outputContainer->Add(fhDeltaPhiAwayParton);
 128
 129   fhDeltaEtaAwayParton  = new TH2F
 130     ("DeltaEtaAwayParton","#eta_{particle} - #eta_{parton} vs p_{T particle}",
 131      200,0,120,200,-2,2);
 132   fhDeltaEtaAwayParton->SetYTitle("#Delta #eta");
 133   fhDeltaEtaAwayParton->SetXTitle("p_{T particle} (GeV/c)");
 134   outputContainer->Add(fhDeltaEtaAwayParton);
 135
 136   fhDeltaPtAwayParton  = new TH2F
 137     ("DeltaPtAwayParton","#p_{T particle} - #p_{T parton} vs p_{T particle}",
 138      200,0,120,100,-10,10);
 139   fhDeltaPtAwayParton->SetYTitle("#Delta #p_{T}");
 140   fhDeltaPtAwayParton->SetXTitle("p_{T particle} (GeV/c)");
 141   outputContainer->Add(fhDeltaPtAwayParton);
 142
 143   fhPtRatAwayParton  = new TH2F
 144     ("PtRatAwayParton","#p_{T parton} / #p_{T particle} vs p_{T particle}",
 145      200,0,120,200,0,5);
 146   fhPtRatAwayParton->SetYTitle("ratio");
 147   fhPtRatAwayParton->SetXTitle("p_{T particle} (GeV/c)");
 148   outputContainer->Add(fhPtRatAwayParton);
 149
 150   return outputContainer;
 151 }
 152
 153 //____________________________________________________________________________
 154 void AliAnaParticlePartonCorrelation::InitParameters()
 155 {
 156
 157   //Initialize the parameters of the analysis.
 158   ;
 159
 160 }
 161
 162 //__________________________________________________________________
 163 void AliAnaParticlePartonCorrelation::Print(const Option_t * opt) const
 164 {
 165
 166   //Print some relevant parameters set for the analysis
 167   if(! opt)
 168     return;
 169
 170 }
 171
 172 //__________________________________________________________________
 173 void  AliAnaParticlePartonCorrelation::MakeAnalysisFillAOD()
 174 {
 175   //Particle-Parton Correlation Analysis, create AODs
 176   //Add partons to the reference list of the trigger particle
 177   //Partons are considered those in the first eight possitions in the stack
 178   //being 0, and 1 the 2 protons, and 6 and 7 the outgoing final partons.
 179
 180   if(GetDebug() > 1){
 181     printf("Begin parton correlation analysis, fill AODs \n");
 182     printf("In particle branch aod entries %d\n", GetAODBranch()->GetEntriesFast());
 183   }
 184
 185   //Loop on stored AOD particles
 186   Int_t naod = GetAODBranch()->GetEntriesFast();
 187   for(Int_t iaod = 0; iaod < naod ; iaod++){
 188     AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
 189
 190     AliStack * stack =  GetMCStack() ;
 191     if(!stack) AliFatal("No Stack available, STOP");
 192
 193     if(stack->GetNtrack() < 8) {
 194       printf("*** small number of particles, not a PYTHIA simulation? ***:  n tracks %d \n", stack->GetNprimary());
 195       continue ;
 196     }
 197
 198     //Fill AOD reference only with partons
 199     TParticle * parton = new TParticle ;
 200
 201     for(Int_t ipr = 0;ipr < 8; ipr ++ ){
 202       parton = stack->Particle(ipr) ;
 203       particle->AddTrack(parton);
 204       //parton->Print();
 205     }
 206
 207   }//Aod branch loop
 208
 209  if(GetDebug() > 1) printf("End parton correlation analysis, fill AODs \n");
 210 }
 211
 212 //__________________________________________________________________
 213 void  AliAnaParticlePartonCorrelation::MakeAnalysisFillHistograms()
 214 {
 215   //Particle-Parton Correlation Analysis, fill histograms
 216   if(GetDebug() > 1){
 217     printf("Begin parton correlation analysis, fill histograms \n");
 218     printf("In particle branch aod entries %d\n", GetAODBranch()->GetEntriesFast());
 219   }
 220
 221   AliStack * stack =  GetMCStack() ;
 222   if(!stack) AliFatal("No Stack available, STOP");
 223
 224   //Loop on stored AOD particles
 225   Int_t naod = GetAODBranch()->GetEntriesFast();
 226   TParticle *  mom =new TParticle ;
 227
 228   for(Int_t iaod = 0; iaod < naod ; iaod++){
 229     AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
 230
 231     Float_t ptTrigg  = particle->Pt();
 232     Float_t phiTrigg = particle->Phi();
 233     Float_t etaTrigg = particle->Eta();
 234     Int_t imom = particle->GetLabel();
 235     Int_t iparent  = 2000;
 236     Int_t iawayparent = -1;
 237
 238     if(!(particle->GetRefTracks()) || (particle->GetRefTracks())->GetEntriesFast() < 7) AliFatal("Reference list with partons not filled, STOP analysis");
 239
 240     //Check and get indeces of mother and parton
 241     if(imom < 8 ) iparent = imom ;   //mother is already a parton
 242     else if (imom <  stack->GetNtrack()) {
 243       mom =  stack->Particle(imom);
 244       iparent=mom->GetFirstMother();
 245       cout<<" iparent "<<iparent<<endl;
 246       while(iparent > 7 ){
 247         mom = stack->Particle(iparent);
 248         imom = iparent ; //Mother label is of the inmediate parton daughter
 249         iparent = mom->GetFirstMother();
 250         cout<<" while iparent "<<iparent<<endl;
 251       }
 252     }
 253
 254     if(GetDebug() > 1) printf("N reference partons %d; labels:  mother %d, parent %d \n", (particle->GetRefTracks())->GetEntriesFast(), imom, iparent);
 255
 256
 257     if(iparent < 0 || iparent > 8) {
 258       if(GetDebug() > 0 ) printf("Failed to find appropriate parton, index %d", iparent);
 259       continue ;
 260     }
 261
 262     //Near parton is the parton that fragmented and created the mother
 263     TParticle * nearParton = (TParticle*) (particle->GetRefTracks())->At(iparent);
 264     Float_t  ptNearParton    = nearParton->Pt();
 265     Float_t  phiNearParton   = nearParton->Phi() ;
 266     Float_t  etaNearParton   = nearParton->Eta() ;
 267
 268     fhDeltaEtaNearParton->Fill(ptTrigg,etaTrigg-etaNearParton);
 269     fhDeltaPhiNearParton->Fill(ptTrigg,phiTrigg-phiNearParton);
 270     fhDeltaPtNearParton->Fill(ptTrigg,ptTrigg-ptNearParton);
 271     fhPtRatNearParton->Fill(ptTrigg,ptNearParton/ptTrigg);
 272
 273     if(iparent == 7) iawayparent =6;
 274     else if(iparent == 6) iawayparent =7;
 275     else{
 276       printf("Parent parton is not final state, skip \n");
 277       continue;
 278     }
 279
 280     //Away parton is the other final parton.
 281     TParticle * awayParton = (TParticle*) (particle->GetRefTracks())->At(iawayparent);
 282     Float_t  ptAwayParton    = awayParton->Pt();
 283     Float_t  phiAwayParton   = awayParton->Phi() ;
 284     Float_t  etaAwayParton   = awayParton->Eta() ;
 285     fhDeltaEtaAwayParton->Fill(ptTrigg,etaTrigg-etaAwayParton);
 286     fhDeltaPhiAwayParton->Fill(ptTrigg,phiTrigg-phiAwayParton);
 287     fhDeltaPtAwayParton->Fill(ptTrigg,ptTrigg-ptAwayParton);
 288     fhPtRatAwayParton->Fill(ptTrigg,ptAwayParton/ptTrigg);
 289
 290
 291   }
 292
 293   if(GetDebug() > 1) printf("End parton correlation analysis, fill histograms \n");
 294
 295 }