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,100,-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,100,-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,100,-50,50);
 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,100,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,100,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, AliAODEvent *event) const
 203 {
 204   //Returns the index of the jet that is opposite to the particle
 205
 206   Int_t njets = event->GetNJets() ;
 207
 208   AliAODJet * jet = 0 ;
 209   Int_t index = -1;
 210   for(Int_t ijet = 0; ijet < njets ; ijet++){
 211     jet = event->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
 234   //Get the event, check if there are AODs available, if not, skip this analysis
 235   AliAODEvent * event = NULL;
 236   if(GetReader()->GetOutputEvent())
 237   {
 238     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
 239   }
 240   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
 241   {
 242     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
 243   }
 244   else
 245   {
 246     if(GetDebug() > 3) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - There are no jets available for this analysis\n");
 247     return;
 248   }
 249
 250   if(!GetInputAODBranch() || !event){
 251     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s > \n",GetInputAODName().Data());
 252     abort();
 253   }
 254
 255   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
 256     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());
 257     abort();
 258   }
 259
 260   Int_t ntrig =  GetInputAODBranch()->GetEntriesFast() ;
 261   if(GetDebug() > 3){
 262     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Begin jet finder  correlation analysis, fill AODs \n");
 263     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", ntrig);
 264     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In jet      branch aod entries %d\n", event->GetNJets());
 265   }
 266
 267   //Loop on stored AOD particles, trigger
 268   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
 269     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 270
 271     //Correlate with jets
 272     Int_t ijet = SelectJet(particle,event);
 273     if(ijet > -1){
 274       if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",ijet);
 275       AliAODJet *jet = event->GetJet(ijet);
 276       particle->SetRefJet(jet);
 277     }
 278   } // input aod loop
 279
 280   if(GetDebug() > 1 ) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
 281 }
 282
 283 //__________________________________________________________________
 284 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
 285 {
 286   //Particle-Jet Correlation Analysis, fill histograms
 287
 288   //Get the event, check if there are AODs available, if not, skip this analysis
 289   AliAODEvent * event = NULL;
 290   if(GetReader()->GetOutputEvent())
 291   {
 292     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
 293   }
 294   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
 295   {
 296     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
 297   }
 298   else {
 299     if(GetDebug() > 3) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - There are no jets available for this analysis\n");
 300     return;
 301   }
 302
 303   if(!GetInputAODBranch() || !event){
 304     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s > \n",GetInputAODName().Data());
 305     abort();
 306   }
 307
 308   Int_t ntrig   =  GetInputAODBranch()->GetEntriesFast() ;
 309   if(GetDebug() > 1){
 310     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Begin jet finder  correlation analysis, fill histograms \n");
 311     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", ntrig);
 312     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In jet output branch aod entries %d\n", event->GetNJets());
 313   }
 314
 315   //Loop on stored AOD particles, trigger
 316   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
 317     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 318
 319     if(OnlyIsolated() && !particlecorr->IsIsolated()) continue;
 320
 321     //Recover the jet correlated, found previously.
 322     AliAODJet   * jet = particlecorr->GetJet();
 323     //If correlation not made before, do it now.
 324     if(fMakeCorrelationInHistoMaker){
 325       //Correlate with jets
 326       Int_t ijet = SelectJet(particlecorr,event);
 327       if(ijet > -1){
 328         if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Jet with index %d selected \n",ijet);
 329         jet = event->GetJet(ijet);
 330         particlecorr->SetRefJet(jet);
 331       }
 332     }
 333
 334     if (!jet) continue ;
 335
 336     //Fill Histograms
 337
 338     Double_t ptTrig = particlecorr->Pt();
 339     Double_t ptJet = jet->Pt();
 340     Double_t phiTrig = particlecorr->Phi();
 341     Double_t phiJet = jet->Phi();
 342     Double_t etaTrig = particlecorr->Eta();
 343     Double_t etaJet = jet->Eta();
 344     //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",
 345     //  ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
 346     fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
 347     fhDeltaPhi->Fill(ptTrig, phiJet-phiTrig);
 348     fhDeltaEta->Fill(ptTrig, etaJet-etaTrig);
 349     fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
 350     fhPt      ->Fill(ptTrig, ptJet);
 351
 352     //Fragmentation function
 353     Float_t      rad = 0, pt = 0, eta = 0, phi = 0;
 354     Int_t        npartcone = 0;
 355     TVector3 p3;
 356
 357     Int_t ntracks =  0;
 358     if(!fUseJetRefTracks)
 359       ntracks =GetAODCTS()->GetEntriesFast();
 360     else //If you want to use jet tracks from JETAN
 361       ntracks =  (jet->GetRefTracks())->GetEntriesFast();
 362
 363     AliAODTrack* track = 0x0 ;
 364     for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
 365       if(!fUseJetRefTracks)
 366         track = (AliAODTrack *) (GetAODCTS()->At(ipr)) ;
 367       else //If you want to use jet tracks from JETAN
 368         track = (AliAODTrack *) ((jet->GetRefTracks())->At(ipr));
 369
 370       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 371       pt    = p3.Pt();
 372       eta  = p3.Eta();
 373       phi  = p3.Phi() ;
 374       if(phi < 0) phi+=TMath::TwoPi();
 375
 376       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 377       rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
 378       if(rad < fConeSize  && pt > fPtThresholdInCone){
 379         //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
 380         npartcone++;
 381         fhFFz ->Fill(ptTrig, pt/ptTrig);
 382         fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
 383         fhFFpt->Fill(ptTrig, pt);
 384       }
 385     }//Tracks
 386     fhNTracksInCone->Fill(ptTrig, npartcone);
 387   }//AOD trigger particle loop
 388   if(GetDebug() > 1) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
 389 }
 390
 391
 392 //__________________________________________________________________
 393 void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
 394 {
 395
 396   //Print some relevant parameters set for the analysis
 397   if(! opt)
 398     return;
 399
 400   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 401   AliAnaPartCorrBaseClass::Print(" ");
 402
 403   printf("Phi trigger-jet        <     %3.2f\n", fDeltaPhiMaxCut) ;
 404   printf("Phi trigger-jet        >     %3.2f\n", fDeltaPhiMinCut) ;
 405   printf("pT Ratio trigger/jet   <     %3.2f\n", fRatioMaxCut) ;
 406   printf("pT Ratio trigger/jet   >     %3.2f\n", fRatioMinCut) ;
 407   printf("fConeSize              =     %3.2f\n", fConeSize) ;
 408   printf("fPtThresholdInCone     =     %3.2f\n", fPtThresholdInCone) ;
 409   printf("fUseJetRefTracks         =     %d\n",    fUseJetRefTracks) ;
 410   printf("fMakeCorrelationInHistoMaker     =     %d\n",    fMakeCorrelationInHistoMaker) ;
 411   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
 412
 413 }
 414