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