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