PWG4/JetTasks/AliAnalysisTaskJetResponse.cxx

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