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