PWG4/AliAnaParticleJetFinderCorrelation.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: AliAnaParticleJetFinderCorrelation.cxx 22232 2007-11-17 16:39:49Z gustavo $ */
  16
  17 //_________________________________________________________________________
  18 // Class for the analysis of particle (direct gamma) -jet (jet found with finder) correlations
  19 //*-- Author: Gustavo Conesa (LNF-INFN)
  20 //////////////////////////////////////////////////////////////////////////////
  21
  22
  23 // --- ROOT system ---
  24 #include "TH2F.h"
  25 #include "Riostream.h"
  26
  27 //---- AliRoot system ----
  28 #include "AliCaloTrackReader.h"
  29 #include "AliAODJet.h"
  30 #include "AliAnaParticleJetFinderCorrelation.h"
  31 #include "AliLog.h"
  32
  33   ClassImp(AliAnaParticleJetFinderCorrelation)
  34
  35
  36 //____________________________________________________________________________
  37   AliAnaParticleJetFinderCorrelation::AliAnaParticleJetFinderCorrelation() :
  38         AliAnaPartCorrBaseClass(),
  39     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fRatioMaxCut(0.),  fRatioMinCut(0.),
  40         fConeSize(0), fPtThresholdInCone(0),fUseJetRefTracks(0), fMakeCorrelationInHistoMaker(0),
  41         fhDeltaEta(0), fhDeltaPhi(0), fhDeltaPt(0), fhPtRatio(0), fhPt(0),
  42         fhFFz(0),fhFFxi(0),fhFFpt(0),fhNTracksInCone(0)
  43 {
  44         //Default Ctor
  45
  46         //Initialize parameters
  47         InitParameters();
  48 }
  49
  50 //____________________________________________________________________________
  51 AliAnaParticleJetFinderCorrelation::AliAnaParticleJetFinderCorrelation(const AliAnaParticleJetFinderCorrelation & pjf) :
  52         AliAnaPartCorrBaseClass(pjf),
  53         fDeltaPhiMaxCut(pjf.fDeltaPhiMaxCut), fDeltaPhiMinCut(pjf.fDeltaPhiMinCut),
  54         fRatioMaxCut(pjf.fRatioMaxCut), fRatioMinCut(pjf.fRatioMinCut),
  55         fConeSize(pjf.fConeSize), fPtThresholdInCone(pjf.fPtThresholdInCone),
  56         fUseJetRefTracks(pjf.fUseJetRefTracks), fMakeCorrelationInHistoMaker(pjf.fMakeCorrelationInHistoMaker),
  57         fhDeltaEta(pjf.fhDeltaEta), fhDeltaPhi(pjf.fhDeltaPhi),
  58         fhDeltaPt(pjf.fhDeltaPt), fhPtRatio(pjf.fhPtRatio), fhPt(pjf.fhPt),
  59         fhFFz(pjf.fhFFz),fhFFxi(pjf.fhFFxi),fhFFpt(pjf.fhFFpt),
  60         fhNTracksInCone(pjf.fhNTracksInCone)
  61 {
  62         // cpy ctor
  63
  64 }
  65
  66 //_________________________________________________________________________
  67 AliAnaParticleJetFinderCorrelation & AliAnaParticleJetFinderCorrelation::operator = (const AliAnaParticleJetFinderCorrelation & pjf)
  68 {
  69         // assignment operator
  70
  71         if(this == &pjf)return *this;
  72         ((AliAnaPartCorrBaseClass *)this)->operator=(pjf);
  73
  74         fDeltaPhiMaxCut    = pjf.fDeltaPhiMaxCut ;
  75         fDeltaPhiMinCut    = pjf.fDeltaPhiMinCut ;
  76         fRatioMaxCut       = pjf.fRatioMaxCut ;
  77         fRatioMinCut       = pjf.fRatioMinCut ;
  78         fConeSize          = pjf.fConeSize ;
  79         fPtThresholdInCone = pjf.fPtThresholdInCone ;
  80         fUseJetRefTracks   = pjf.fUseJetRefTracks ;
  81         fMakeCorrelationInHistoMaker = pjf.fMakeCorrelationInHistoMaker ;
  82
  83         //Histograms
  84         fhDeltaEta = pjf.fhDeltaEta;
  85         fhDeltaPhi = pjf.fhDeltaPhi;
  86         fhDeltaPt  = pjf.fhDeltaPt;
  87         fhPtRatio  = pjf.fhPtRatio;
  88         fhPt       = pjf.fhPt;
  89
  90         fhFFz   = pjf.fhFFz;
  91         fhFFxi  = pjf.fhFFxi;
  92         fhFFpt  = pjf.fhFFpt;
  93         fhNTracksInCone = pjf.fhNTracksInCone;
  94
  95   return *this;
  96
  97 }
  98
  99 //____________________________________________________________________________
 100 //AliAnaParticleJetFinderCorrelation::~AliAnaParticleJetFinderCorrelation()
 101 //{
 102 //   // Remove all pointers except analysis output pointers.
 103 //
 104 //}
 105
 106
 107 //________________________________________________________________________
 108 TList *  AliAnaParticleJetFinderCorrelation::GetCreateOutputObjects()
 109 {
 110   // Create histograms to be saved in output file and
 111   // store them in fOutputContainer
 112
 113   if(GetDebug()>1) printf("Init histograms \n");
 114
 115   TList * outputContainer = new TList() ;
 116   outputContainer->SetName("ParticleJetFinderHistos") ;
 117
 118   fhDeltaPhi  = new TH2F("DeltaPhi","#phi_{jet} - #phi_{trigger} vs p_{T trigger}",120,0,120,200,-4,4);
 119   fhDeltaPhi->SetYTitle("#Delta #phi");
 120   fhDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 121   outputContainer->Add(fhDeltaPhi);
 122
 123   fhDeltaEta  = new TH2F("DeltaEta","#eta_{jet} - #eta_{trigger} vs p_{T trigger}",120,0,120,200,-5,5);
 124   fhDeltaEta->SetYTitle("#Delta #eta");
 125   fhDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 126   outputContainer->Add(fhDeltaEta);
 127
 128   fhDeltaPt  = new TH2F("DeltaPt","#p_{T trigger} - #p_{T jet} vs p_{T trigger}",120,0,120,200,-100,100);
 129   fhDeltaPt->SetYTitle("#Delta #p_{T}");
 130   fhDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
 131   outputContainer->Add(fhDeltaPt);
 132
 133   fhPtRatio  = new TH2F("PtRatio","#p_{T jet} / #p_{T trigger} vs p_{T trigger}",120,0,120,200,0,2.);
 134   fhPtRatio->SetYTitle("ratio");
 135   fhPtRatio->SetXTitle("p_{T trigger} (GeV/c)");
 136   outputContainer->Add(fhPtRatio);
 137
 138   fhPt  = new TH2F("Pt","#p_{T jet} vs p_{T trigger}",120,0,120,120,0,120);
 139   fhPt->SetYTitle("#p_{T jet}(GeV/c)");
 140   fhPt->SetXTitle("p_{T trigger} (GeV/c)");
 141   outputContainer->Add(fhPt);
 142
 143   fhFFz  = new TH2F("FFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2);
 144   fhFFz->SetYTitle("z");
 145   fhFFz->SetXTitle("p_{T trigger}");
 146   outputContainer->Add(fhFFz) ;
 147
 148   fhFFxi  = new TH2F("FFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.);
 149   fhFFxi->SetYTitle("#xi");
 150   fhFFxi->SetXTitle("p_{T trigger}");
 151   outputContainer->Add(fhFFxi) ;
 152
 153   fhFFpt  = new TH2F("FFpt","#xi = p_{T i charged}) vs p_{T trigger}", 120,0.,120.,200,0.,50.);
 154   fhFFpt->SetYTitle("p_{T charged hadron}");
 155   fhFFpt->SetXTitle("p_{T trigger}");
 156   outputContainer->Add(fhFFpt) ;
 157
 158   fhNTracksInCone  = new TH2F("NTracksInCone","#xi = p_{T i charged}) vs p_{T trigger}", 120,0.,120.,200,0.,50.);
 159   fhNTracksInCone->SetYTitle("p_{T charged hadron}");
 160   fhNTracksInCone->SetXTitle("p_{T trigger}");
 161   outputContainer->Add(fhNTracksInCone) ;
 162
 163   return outputContainer;
 164 }
 165
 166  //____________________________________________________________________________
 167 void AliAnaParticleJetFinderCorrelation::InitParameters()
 168 {
 169
 170   //Initialize the parameters of the analysis.
 171
 172   fDeltaPhiMinCut    = 1.5 ;
 173   fDeltaPhiMaxCut    = 4.5 ;
 174   fRatioMaxCut       = 1.2 ;
 175   fRatioMinCut       = 0.3 ;
 176   fConeSize          = 0.3 ;
 177   fPtThresholdInCone = 0.5 ;
 178   fUseJetRefTracks   = kFALSE ;
 179   fMakeCorrelationInHistoMaker = kFALSE ;
 180   }
 181
 182 //__________________________________________________________________
 183 Int_t  AliAnaParticleJetFinderCorrelation::SelectJet(AliAODParticleCorrelation * particle)
 184 {
 185         //Returns the index of the jet that is opposite to the particle
 186
 187         Int_t njets = (GetReader()->GetAOD())->GetNJets() ;
 188         cout<<">> Loop on Jets"<<endl;
 189         AliAODJet * jet = new AliAODJet ;
 190         Int_t index = -1;
 191         for(Int_t ijet = 0; ijet < njets ; ijet++){
 192                 jet = (GetReader()->GetAOD())->GetJet(ijet) ;
 193                 Float_t dphi  = TMath::Abs(particle->Phi()-jet->Phi());
 194                 Float_t ratio = jet->Pt()/particle->Pt();
 195                 if(GetDebug() > 3)
 196                         printf("Jet %d, Ratio pT %2.3f, Delta phi %2.3f\n",ijet,ratio,dphi);
 197                 Float_t dphiprev= 10000;
 198                 if((dphi > fDeltaPhiMinCut) && (dphi<fDeltaPhiMaxCut) &&
 199                         (ratio > fRatioMinCut) && (ratio < fRatioMaxCut)  &&
 200                         (TMath::Abs(dphi-3.14159) < TMath::Abs(dphiprev-3.14159))){//In case there is more than one jet select the most opposite.
 201                         dphiprev = dphi;
 202                         index = ijet ;
 203                 }//Selection cuts
 204         }//AOD Jet loop
 205
 206         return index ;
 207
 208 }
 209
 210 //__________________________________________________________________
 211 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD()
 212 {
 213         //Particle-Jet Correlation Analysis, fill AODs
 214
 215         Int_t ntrig =  GetAODBranch()->GetEntriesFast() ;
 216     if(GetDebug() > -3){
 217     printf("Begin jet finder  correlation analysis, fill AODs \n");
 218     printf("In particle branch aod entries %d\n", ntrig);
 219     printf("In jet      branch aod entries %d\n", (GetReader()->GetAOD())->GetNJets());
 220         }
 221
 222         //Loop on stored AOD particles, trigger
 223         for(Int_t iaod = 0; iaod < ntrig ; iaod++){
 224                 AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
 225                 //cout<<">> Test reference "<<particle->GetRefJet()<<endl;
 226                         //Correlate with jets
 227                 Int_t ijet = SelectJet(particle);
 228                 if(ijet > -1){
 229                         if(GetDebug() > 2) printf ("Jet with index %d selected \n",ijet);
 230                         AliAODJet *jet = (GetReader()->GetAOD())->GetJet(ijet);
 231                         particle->SetRefJet(jet);
 232                 }
 233         }
 234
 235         if(GetDebug() >1)printf("End of jet leading cone analysis, fill AODs \n");
 236 }
 237
 238 //__________________________________________________________________
 239 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
 240 {
 241         //Particle-Jet Correlation Analysis, fill histograms
 242         Int_t ntrig   =  GetAODBranch()->GetEntriesFast() ;
 243         if(GetDebug() > 1){
 244                 printf("Begin jet finder  correlation analysis, fill histograms \n");
 245                 printf("In particle branch aod entries %d\n", ntrig);
 246                 printf("In jet      branch aod entries %d\n", (GetReader()->GetAOD())->GetNJets());
 247         }
 248
 249         //Loop on stored AOD particles, trigger
 250         for(Int_t iaod = 0; iaod < ntrig ; iaod++){
 251                 AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
 252                 //cout<<">> Test reference "<<particle->GetRefJet()<<endl;
 253
 254                 //Recover the jet correlated, found previously.
 255                 AliAODJet       * jet = particle->GetJet();
 256                 //If correlation not made before, do it now.
 257                 if(fMakeCorrelationInHistoMaker){
 258                         //Correlate with jets
 259                         Int_t ijet = SelectJet(particle);
 260                         if(ijet > -1){
 261                                 if(GetDebug() > 2) printf ("Jet with index %d selected \n",ijet);
 262                                 jet = (GetReader()->GetAOD())->GetJet(ijet);
 263                                 particle->SetRefJet(jet);
 264                         }
 265                 }
 266
 267                 if (!jet) continue ;
 268
 269                 //Fill Histograms
 270
 271                 Double_t ptTrig = particle->Pt();
 272                 Double_t ptJet = jet->Pt();
 273                 Double_t phiTrig = particle->Phi();
 274                 Double_t phiJet = jet->Phi();
 275                 Double_t etaTrig = particle->Eta();
 276                 Double_t etaJet = jet->Eta();
 277                 //printf("pT trigger %2.3f, pT jet %2.3f, Delta phi %2.3f, Delta eta %2.3f, Delta pT %2.3f, ratio %2.3f \n",
 278                 //      ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
 279                 fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
 280                 fhDeltaPhi->Fill(ptTrig, phiJet-phiTrig);
 281                 fhDeltaEta->Fill(ptTrig, etaJet-etaTrig);
 282                 fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
 283                 fhPt      ->Fill(ptTrig, ptJet);
 284
 285                 //Fragmentation function
 286                 Float_t  rad = 0, pt = 0, eta = 0, phi = 0;
 287                 Int_t    npartcone = 0;
 288                 TVector3 p3;
 289                 AliAODTrack* track = new AliAODTrack ;
 290
 291                 Int_t ntracks =  0;
 292                 if(!fUseJetRefTracks)
 293                         ntracks =GetAODCTS()->GetEntriesFast();
 294                 else //If you want to use jet tracks from JETAN
 295                         ntracks =  (jet->GetRefTracks())->GetEntriesFast();
 296
 297                 for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
 298                         if(!fUseJetRefTracks)
 299                                 track = dynamic_cast<AliAODTrack *>(GetAODCTS()->At(ipr)) ;
 300                         else //If you want to use jet tracks from JETAN
 301                                 track = dynamic_cast<AliAODTrack *>((jet->GetRefTracks())->At(ipr));
 302
 303                         p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 304                         pt    = p3.Pt();
 305                         eta  = p3.Eta();
 306                         phi  = p3.Phi() ;
 307                         if(phi < 0) phi+=TMath::TwoPi();
 308
 309                         //Check if there is any particle inside cone with pt larger than  fPtThreshold
 310                         rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
 311                         if(rad < fConeSize  && pt > fPtThresholdInCone){
 312                                 //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
 313                                 npartcone++;
 314                                 fhFFz ->Fill(ptTrig, pt/ptTrig);
 315                                 fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
 316                                 fhFFpt->Fill(ptTrig, pt);
 317                         }
 318                 }//Tracks
 319                 fhNTracksInCone->Fill(ptTrig, npartcone);
 320         }//AOD trigger particle loop
 321         if(GetDebug() >1) printf("End of jet leading cone analysis, fill histograms \n");
 322 }
 323
 324
 325 //__________________________________________________________________
 326 void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
 327 {
 328
 329         //Print some relevant parameters set for the analysis
 330         if(! opt)
 331                 return;
 332
 333         printf("Phi trigger-jet        <     %3.2f\n", fDeltaPhiMaxCut) ;
 334         printf("Phi trigger-jet        >     %3.2f\n", fDeltaPhiMinCut) ;
 335         printf("pT Ratio trigger/jet   <     %3.2f\n", fRatioMaxCut) ;
 336         printf("pT Ratio trigger/jet   >     %3.2f\n", fRatioMinCut) ;
 337         printf("fConeSize              =     %3.2f\n", fConeSize) ;
 338         printf("fPtThresholdInCone     =     %3.2f\n", fPtThresholdInCone) ;
 339         printf("fUseJetRefTracks           =     %d\n",    fUseJetRefTracks) ;
 340         printf("fMakeCorrelationInHistoMaker       =     %d\n",    fMakeCorrelationInHistoMaker) ;
 341 }
 342