PWGJE/AliAnalysisTaskJetResponse.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 // task compares jets in two branches,
  18 // written for analysis of jet embedding in HI events
  19 //
  20 // newer class: AliAnalysisTaskJetResponseV2
  21 //
  22
  23 #include "TChain.h"
  24 #include "TTree.h"
  25 #include "TMath.h"
  26 #include "TH1F.h"
  27 #include "TH2F.h"
  28 #include "TH3F.h"
  29 #include "TCanvas.h"
  30
  31 #include "AliLog.h"
  32
  33 #include "AliAnalysisTask.h"
  34 #include "AliAnalysisManager.h"
  35
  36 #include "AliVEvent.h"
  37 #include "AliESDEvent.h"
  38 #include "AliESDInputHandler.h"
  39 #include "AliCentrality.h"
  40 #include "AliAnalysisHelperJetTasks.h"
  41 #include "AliInputEventHandler.h"
  42
  43 #include "AliAODEvent.h"
  44 #include "AliAODJet.h"
  45
  46 #include "AliAnalysisTaskJetResponse.h"
  47
  48 ClassImp(AliAnalysisTaskJetResponse)
  49
  50 AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse() :
  51   AliAnalysisTaskSE(),
  52   fESD(0x0),
  53   fAOD(0x0),
  54   fOfflineTrgMask(AliVEvent::kAny),
  55   fMinContribVtx(1),
  56   fVtxZMin(-8.),
  57   fVtxZMax(8.),
  58   fEvtClassMin(1),
  59   fEvtClassMax(4),
  60   fCentMin(0.),
  61   fCentMax(100.),
  62   fNInputTracksMin(0),
  63   fNInputTracksMax(-1),
  64   fReactionPlaneBin(-1),
  65   fJetEtaMin(-.5),
  66   fJetEtaMax(.5),
  67   fJetPtMin(20.),
  68   fJetPtFractionMin(0.5),
  69   fNMatchJets(4),
  70   //fJetDeltaEta(.4),
  71   //fJetDeltaPhi(.4),
  72   fEvtClassMode(0),
  73   fkNbranches(2),
  74   fkEvtClasses(12),
  75   fOutputList(0x0),
  76   fHistEvtSelection(0x0),
  77   fHistEvtClass(0x0),
  78   fHistCentrality(0x0),
  79   fHistNInputTracks(0x0),
  80   //fHistNInputTracks2(0x0),
  81   fHistCentVsTracks(0x0),
  82   fHistReactionPlane(0x0),
  83   fHistReactionPlaneWrtJets(0x0),
  84   fHistPtJet(new TH1F*[fkNbranches]),
  85   fHistEtaPhiJet(new TH2F*[fkNbranches]),
  86   fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
  87   fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
  88   fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
  89   //fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
  90   fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
  91   fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
  92   fHistPtFraction(new TH2F*[fkEvtClasses]),
  93   fHistPtResponse(new TH2F*[fkEvtClasses]),
  94   fHistPtSmearing(new TH2F*[fkEvtClasses]),
  95   fHistDeltaR(new TH2F*[fkEvtClasses]),
  96   fHistArea(new TH3F*[fkEvtClasses]),
  97   //fHistJetPtArea(new TH2F*[fkEvtClasses]),
  98   fHistDeltaArea(new TH2F*[fkEvtClasses]),
  99   fHistJetPtDeltaArea(new TH2F*[fkEvtClasses])
 100 {
 101   // default Constructor
 102
 103   fJetBranchName[0] = "";
 104   fJetBranchName[1] = "";
 105
 106   fListJets[0] = new TList;
 107   fListJets[1] = new TList;
 108 }
 109
 110 AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
 111   AliAnalysisTaskSE(name),
 112   fESD(0x0),
 113   fAOD(0x0),
 114   fOfflineTrgMask(AliVEvent::kAny),
 115   fMinContribVtx(1),
 116   fVtxZMin(-8.),
 117   fVtxZMax(8.),
 118   fEvtClassMin(1),
 119   fEvtClassMax(4),
 120   fCentMin(0.),
 121   fCentMax(100.),
 122   fNInputTracksMin(0),
 123   fNInputTracksMax(-1),
 124   fReactionPlaneBin(-1),
 125   fJetEtaMin(-.5),
 126   fJetEtaMax(.5),
 127   fJetPtMin(20.),
 128   fJetPtFractionMin(0.5),
 129   fNMatchJets(4),
 130   fEvtClassMode(0),
 131   fkNbranches(2),
 132   fkEvtClasses(12),
 133   fOutputList(0x0),
 134   fHistEvtSelection(0x0),
 135   fHistEvtClass(0x0),
 136   fHistCentrality(0x0),
 137   fHistNInputTracks(0x0),
 138   //fHistNInputTracks2(0x0),
 139   fHistCentVsTracks(0x0),
 140   fHistReactionPlane(0x0),
 141   fHistReactionPlaneWrtJets(0x0),
 142   fHistPtJet(new TH1F*[fkNbranches]),
 143   fHistEtaPhiJet(new TH2F*[fkNbranches]),
 144   fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
 145   fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
 146   fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
 147   //fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
 148   fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
 149   fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
 150   fHistPtFraction(new TH2F*[fkEvtClasses]),
 151   fHistPtResponse(new TH2F*[fkEvtClasses]),
 152   fHistPtSmearing(new TH2F*[fkEvtClasses]),
 153   fHistDeltaR(new TH2F*[fkEvtClasses]),
 154   fHistArea(new TH3F*[fkEvtClasses]),
 155   //fHistJetPtArea(new TH2F*[fkEvtClasses]),
 156   fHistDeltaArea(new TH2F*[fkEvtClasses]),
 157   fHistJetPtDeltaArea(new TH2F*[fkEvtClasses])
 158 {
 159   // Constructor
 160
 161   fJetBranchName[0] = "";
 162   fJetBranchName[1] = "";
 163
 164   fListJets[0] = new TList;
 165   fListJets[1] = new TList;
 166
 167   DefineOutput(1, TList::Class());
 168 }
 169
 170 AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
 171 {
 172   delete fListJets[0];
 173   delete fListJets[1];
 174 }
 175
 176 void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
 177 {
 178   fJetBranchName[0] = branch1;
 179   fJetBranchName[1] = branch2;
 180 }
 181
 182 void AliAnalysisTaskJetResponse::Init()
 183 {
 184
 185    // check for jet branches
 186    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 187       AliError("Jet branch name not set.");
 188    }
 189
 190 }
 191
 192 void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
 193 {
 194   // Create histograms
 195   // Called once
 196   OpenFile(1);
 197   if(!fOutputList) fOutputList = new TList;
 198   fOutputList->SetOwner(kTRUE);
 199
 200   Bool_t oldStatus = TH1::AddDirectoryStatus();
 201   TH1::AddDirectory(kFALSE);
 202
 203
 204   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 205   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 206   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 207   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 208   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 209   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 210   fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
 211
 212   fHistEvtClass   = new TH1I("fHistEvtClass", "event classes (centrality) ", fkEvtClasses, -0.5, fkEvtClasses-0.5);
 213   fHistCentrality = new TH1F("fHistCentrality", "event centrality", 100, 0., 100.);
 214   fHistNInputTracks  = new TH1F("fHistNInputTracks", "nb. of input tracks", 400, 0., 4000.);
 215   //fHistNInputTracks2 = new TH2F("fHistNInputTracks2", "nb. of input tracks;from B0;from Bx", 400, 0., 4000., 400, 0., 4000.);
 216   fHistCentVsTracks  = new TH2F("fHistCentVsTracks", "centrality vs. nb. of input tracks;centrality;input tracks", 100, 0., 100., 400, 0., 4000.);
 217   fHistReactionPlane = new TH1F("fHistReactionPlane", "reaction plane", 30, 0., TMath::Pi());
 218   fHistReactionPlaneWrtJets = new TH1F("fHistReactionPlaneWrtJets", "reaction plane wrt jets", 48, -TMath::Pi(), TMath::Pi());
 219
 220   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
 221     fHistPtJet[iBranch] = new TH1F(Form("fHistPtJet%i", iBranch),
 222                                           Form("p_{T} of jets from branch %i;p_{T} (GeV/c);counts", iBranch),
 223                                           250, 0., 250.);
 224     fHistPtJet[iBranch]->SetMarkerStyle(kFullCircle);
 225
 226     fHistEtaPhiJet[iBranch]    = new TH2F(Form("fHistEtaPhiJet%i", iBranch),
 227                                                          Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
 228                                                          28, -0.7, 0.7, 100, 0., 2*TMath::Pi());
 229         fHistEtaPhiJetCut[iBranch] = new TH2F(Form("fHistEtaPhiJetCut%i", iBranch),
 230                                                          Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
 231                                                          28, -0.7, 0.7, 100, 0., 2*TMath::Pi());
 232
 233   }
 234
 235   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
 236     fHistDeltaEtaDeltaPhiJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJet%i",iEvtClass),
 237                                                                       "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
 238                                                                       101, -1.01, 1.01, 101, -1.01, 1.01);
 239     fHistDeltaEtaDeltaPhiJetCut[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJetCut%i",iEvtClass),
 240                                                                       "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
 241                                                                       101, -1.01, 1.01, 101, -1.01, 1.01);
 242     //fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass] = new TH2F("fHistDeltaEtaDeltaPhiJetNOMatching",
 243         //                                                                       "#Delta#eta - #Delta#phi of jets which do not match;#Delta#eta;#Delta#phi",
 244         //                                                                       100, -2., 2., 100, TMath::Pi(), TMath::Pi());
 245     fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
 246                                                               250, 0., 250., 250, 0., 250.);
 247     fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;#Deltap_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
 248                                                               161, -80.5, 80.5, 250, 0., 250.);
 249     fHistDeltaR[iEvtClass]     = new TH2F(Form("hist_DeltaR%i",iEvtClass), "#DeltaR of matched jets;#Deltap_{T};#DeltaR", 161, -80.5,80.5, 85, 0.,3.4);
 250     fHistArea[iEvtClass]       = new TH3F(Form("hist_Area%i",iEvtClass), "jet area;probe p_{T};#Deltap_{T};jet area", 250, 0., 250., 161, -80.5,80.5, 100, 0.,1.0);
 251         //fHistJetPtArea[iEvtClass]       = new TH2F(Form("hist_JetPtArea%i",iEvtClass), "jet area vs. probe p_{T};jet p_{T};jet area", 250, 0.,250., 100, 0.,1.0);
 252     fHistDeltaArea[iEvtClass]  = new TH2F(Form("hist_DeltaArea%i",iEvtClass), "#DeltaArea of matched jets;#Deltap_{T};#Deltaarea", 161, -80.5, 80.5, 32, -.8, .8);
 253         fHistJetPtDeltaArea[iEvtClass]  = new TH2F(Form("hist_JetPtDeltaArea%i",iEvtClass), "#DeltaArea of matched jets vs. probe jet p_{T};#Deltap_{T};#Deltaarea", 250, 0., 250., 32, -.8, .8);
 254
 255     fHistDeltaEtaEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaJet%i", iEvtClass),
 256                                                 Form("#eta - #Delta#eta of matched jets from event class %i;#eta;#Delta#eta", iEvtClass),
 257                                                 28, -.7, .7, 101, -1.01, 1.01);
 258     fHistDeltaPtEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaJet%i", iEvtClass),
 259                                                 Form("#eta - #Deltap_{T} of matched jets from event class %i;#eta;#Deltap_{T}", iEvtClass),
 260                                                 28, -.7, .7, 161, -80.5, 80.5);
 261     fHistPtFraction[iEvtClass] = new TH2F(Form("fHistPtFraction%i", iEvtClass),
 262                                               Form("fraction from embedded jet in reconstructed jet;fraction (event class %i);p_{T}^{emb}", iEvtClass),
 263                                                                                   52, 0., 1.04, 250, 0., 250.);
 264   }
 265
 266
 267   fOutputList->Add(fHistEvtSelection);
 268   fOutputList->Add(fHistEvtClass);
 269   fOutputList->Add(fHistCentrality);
 270   fOutputList->Add(fHistNInputTracks);
 271   fOutputList->Add(fHistCentVsTracks);
 272   fOutputList->Add(fHistReactionPlane);
 273   fOutputList->Add(fHistReactionPlaneWrtJets);
 274
 275
 276   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
 277     fOutputList->Add(fHistPtJet[iBranch]);
 278     fOutputList->Add(fHistEtaPhiJet[iBranch]);
 279         fOutputList->Add(fHistEtaPhiJetCut[iBranch]);
 280   }
 281
 282   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
 283     fOutputList->Add(fHistDeltaEtaDeltaPhiJet[iEvtClass]);
 284     fOutputList->Add(fHistDeltaEtaDeltaPhiJetCut[iEvtClass]);
 285     //fOutputList->Add(fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass]);
 286         fOutputList->Add(fHistDeltaEtaEtaJet[iEvtClass]);
 287     fOutputList->Add(fHistDeltaPtEtaJet[iEvtClass]);
 288         fOutputList->Add(fHistPtFraction[iEvtClass]);
 289
 290     fOutputList->Add(fHistPtResponse[iEvtClass]);
 291     fOutputList->Add(fHistPtSmearing[iEvtClass]);
 292     fOutputList->Add(fHistDeltaR[iEvtClass]);
 293     fOutputList->Add(fHistArea[iEvtClass]);
 294         //fOutputList->Add(fHistJetPtArea[iEvtClass]);
 295     fOutputList->Add(fHistDeltaArea[iEvtClass]);
 296         fOutputList->Add(fHistJetPtDeltaArea[iEvtClass]);
 297   }
 298
 299   // =========== Switch on Sumw2 for all histos ===========
 300   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 301     TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 302     if (h1){
 303       h1->Sumw2();
 304       continue;
 305     }
 306   }
 307   TH1::AddDirectory(oldStatus);
 308
 309   PostData(1, fOutputList);
 310 }
 311
 312 void AliAnalysisTaskJetResponse::UserExec(Option_t *)
 313 {
 314   // load events, apply event cuts, then compare leading jets
 315   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 316       AliError("Jet branch name not set.");
 317       return;
 318    }
 319
 320   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 321   if (!fESD) {
 322     AliError("ESD not available");
 323     return;
 324   }
 325   fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
 326   if (!fAOD) {
 327     AliError("AOD not available");
 328     return;
 329   }
 330
 331   // -- event selection --
 332   fHistEvtSelection->Fill(1); // number of events before event selection
 333
 334   // physics selection
 335   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 336     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 337   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 338     if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 339     fHistEvtSelection->Fill(2);
 340     PostData(1, fOutputList);
 341     return;
 342   }
 343
 344   // vertex selection
 345   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 346   Int_t nTracksPrim = primVtx->GetNContributors();
 347   if ((nTracksPrim < fMinContribVtx) ||
 348       (primVtx->GetZ() < fVtxZMin) ||
 349       (primVtx->GetZ() > fVtxZMax) ){
 350     if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 351     fHistEvtSelection->Fill(3);
 352     PostData(1, fOutputList);
 353     return;
 354   }
 355
 356   // event class selection (from jet helper task)
 357   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 358   if(fDebug) Printf("Event class %d", eventClass);
 359   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 360     fHistEvtSelection->Fill(4);
 361     PostData(1, fOutputList);
 362     return;
 363   }
 364
 365   // centrality selection
 366   AliCentrality *cent = fESD->GetCentrality();
 367   Float_t centValue = cent->GetCentralityPercentile("V0M");
 368   if(fDebug) printf("centrality: %f\n", centValue);
 369   if (centValue < fCentMin || centValue > fCentMax){
 370     fHistEvtSelection->Fill(4);
 371     PostData(1, fOutputList);
 372     return;
 373   }
 374
 375
 376   // multiplicity due to input tracks
 377   Int_t nInputTracks = 0;
 378
 379   TString jbname(fJetBranchName[1]);
 380   for(Int_t i=1; i<=3; ++i){
 381     if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
 382   }
 383   // use only HI event
 384   if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
 385   if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
 386
 387   if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
 388   TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
 389
 390   if (tmpAODjets) {
 391     for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
 392       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
 393           if(!jet) continue;
 394           TRefArray *trackList = jet->GetRefTracks();
 395           Int_t nTracks = trackList->GetEntriesFast();
 396           nInputTracks += nTracks;
 397           if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
 398     }
 399   }
 400   if(fDebug) Printf("---> input tracks: %d", nInputTracks);
 401
 402   if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
 403        fHistEvtSelection->Fill(5);
 404        PostData(1, fOutputList);
 405        return;
 406   }
 407
 408   if(fEvtClassMode==1){ //multiplicity
 409      fHistEvtClass->SetTitle("event classes (multiplicity)");
 410      eventClass = fkEvtClasses-1 - (Int_t)(nInputTracks/250.);
 411      if(eventClass<0) eventClass = 0;
 412   }
 413
 414
 415
 416   fHistEvtSelection->Fill(0); // accepted events
 417   fHistEvtClass->Fill(eventClass);
 418   fHistCentrality->Fill(centValue);
 419   fHistNInputTracks->Fill(nInputTracks);
 420   fHistCentVsTracks->Fill(centValue,nInputTracks);
 421   Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE);
 422   // -- end event selection --
 423
 424
 425
 426   // fetch jets
 427   TClonesArray *aodJets[2];
 428   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
 429   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
 430
 431   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 432     fListJets[iJetType]->Clear();
 433     if (!aodJets[iJetType]) continue;
 434
 435     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 436       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 437       if (jet) fListJets[iJetType]->Add(jet);
 438     }
 439     fListJets[iJetType]->Sort();
 440   }
 441
 442   // jet matching
 443   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
 444   static TArrayF aPtFraction(fListJets[0]->GetEntries());
 445   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
 446   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
 447
 448   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
 449   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
 450                                             fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
 451                                             aMatchIndex, aPtFraction, fDebug);
 452
 453   // loop over matched jets
 454   Int_t ir = -1; // index of matched reconstruced jet
 455   Float_t fraction = 0.;
 456   AliAODJet *jet[2]  = { 0x0, 0x0 };
 457   Float_t jetEta[2]  = { 0., 0. };
 458   Float_t jetPhi[2]  = { 0., 0. };
 459   Float_t jetPt[2]   = { 0., 0. };
 460   Float_t jetArea[2] = { 0., 0. };
 461
 462   for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
 463      ir = aMatchIndex[ig];
 464          if(ir<0) continue;
 465          fraction = aPtFraction[ig];
 466
 467          // fetch jets
 468          jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
 469          jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
 470          if(!jet[0] || !jet[1]) continue;
 471
 472          for(Int_t i=0; i<fkNbranches; ++i){
 473             jetEta[i]  = jet[i]->Eta();
 474                 jetPhi[i]  = jet[i]->Phi();
 475                 jetPt[i]   = jet[i]->Pt();
 476                 jetArea[i] = jet[i]->EffectiveAreaCharged();
 477          }
 478
 479          // reaction plane;
 480          if(fReactionPlaneBin>=0){
 481             Int_t rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[1]), 3);
 482             if(rpBin!=fReactionPlaneBin) continue;
 483          }
 484          fHistReactionPlane->Fill(rp);
 485          fHistReactionPlaneWrtJets->Fill(TVector2::Phi_mpi_pi(rp-jetPhi[1]));
 486
 487
 488
 489          if(eventClass > -1 && eventClass < fkEvtClasses){
 490         fHistPtFraction[eventClass] -> Fill(fraction, jetPt[0]);
 491          }
 492
 493          if(fraction<fJetPtFractionMin) continue;
 494
 495          // calculate parameters of associated jets
 496          Float_t deltaPt    = jetPt[1]-jetPt[0];
 497          Float_t deltaEta   = jetEta[1]-jetEta[0];
 498          Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
 499          Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
 500          Float_t deltaArea  = jetArea[1]-jetArea[0];
 501
 502          // fill histograms before acceptance cut
 503          for(Int_t i=0; i<fkNbranches; ++i){
 504              fHistEtaPhiJet[i] -> Fill(jetEta[i], jetPhi[i]);
 505          }
 506          if(eventClass > -1 && eventClass < fkEvtClasses){
 507              fHistDeltaEtaDeltaPhiJet[eventClass] -> Fill(deltaEta, deltaPhi);
 508          }
 509
 510          // jet acceptance + minimum pT check
 511          if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
 512             jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
 513
 514                 if(fDebug){
 515                 Printf("Jet not in eta acceptance.");
 516                         Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
 517                         Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
 518                 }
 519                 continue;
 520      }
 521          if(jetPt[1] < fJetPtMin){
 522             if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
 523             continue;
 524          }
 525
 526
 527
 528          // fill histograms
 529          for(Int_t i=0; i<fkNbranches; ++i){
 530              fHistPtJet[i]        -> Fill(jetPt[i]);
 531              fHistEtaPhiJetCut[i] -> Fill(jetEta[i], jetPhi[i]);
 532          }
 533
 534          if(eventClass > -1 && eventClass < fkEvtClasses){
 535                  fHistDeltaEtaDeltaPhiJetCut[eventClass] -> Fill(deltaEta, deltaPhi);
 536
 537                  fHistPtResponse[eventClass]             -> Fill(jetPt[1], jetPt[0]);
 538                  fHistPtSmearing[eventClass]             -> Fill(deltaPt,  jetPt[0]);
 539
 540                  fHistDeltaPtEtaJet[eventClass]          -> Fill(jetEta[0], deltaPt);
 541                  fHistDeltaEtaEtaJet[eventClass]         -> Fill(jetEta[0], deltaEta);
 542
 543                  fHistDeltaR[eventClass]                 -> Fill(deltaPt, deltaR);
 544                  fHistArea[eventClass]                   -> Fill(jetPt[0], deltaPt, jetArea[1]);
 545                                  //fHistJetPtArea[eventClass]                   -> Fill(jetPt[0], jetArea[1]);
 546                  fHistDeltaArea[eventClass]              -> Fill(deltaPt, deltaArea);
 547                                  fHistJetPtDeltaArea[eventClass]              -> Fill(jetPt[0], deltaArea);
 548
 549          }
 550   }
 551
 552   PostData(1, fOutputList);
 553 }
 554
 555 void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
 556 {
 557   // Draw result to the screen
 558   // Called once at the end of the query
 559
 560   if (!GetOutputData(1))
 561     return;
 562 }
 563