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