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