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