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   fJetEtaMin(-.5),
  41   fJetEtaMax(.5),
  42   fJetPtMin(20.),
  43   fJetPtFractionMin(0.5),
  44   fNMatchJets(4),
  45   //fJetDeltaEta(.4),
  46   //fJetDeltaPhi(.4),
  47   fkNbranches(2),
  48   fkEvtClasses(5),
  49   fOutputList(0x0),
  50   fHistEvtSelection(0x0),
  51   fHistEvtClass(0x0),
  52   fHistCentrality(0x0),
  53   fHistPtJet(new TH1F*[fkNbranches]),
  54   fHistEtaPhiJet(new TH2F*[fkNbranches]),
  55   fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
  56   fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
  57   fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
  58   fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
  59   fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
  60   fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
  61   fHistPtFraction(new TH2F*[fkEvtClasses]),
  62   fHistPtPtExtra(0x0),
  63   fHistPtResponse(new TH2F*[fkEvtClasses]),
  64   fHistPtSmearing(new TH2F*[fkEvtClasses]),
  65   fHistDeltaR(new TH2F*[fkEvtClasses]),
  66   fHistArea(new TH2F*[fkEvtClasses]),
  67   fHistDeltaArea(new TH2F*[fkEvtClasses])
  68 {
  69   // default Constructor
  70
  71   fJetBranchName[0] = "";
  72   fJetBranchName[1] = "";
  73
  74   fListJets[0] = new TList;
  75   fListJets[1] = new TList;
  76 }
  77
  78 AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
  79   AliAnalysisTaskSE(name),
  80   fESD(0x0),
  81   fAOD(0x0),
  82   fOfflineTrgMask(AliVEvent::kAny),
  83   fMinContribVtx(1),
  84   fVtxZMin(-8.),
  85   fVtxZMax(8.),
  86   fEvtClassMin(1),
  87   fEvtClassMax(4),
  88   fCentMin(0.),
  89   fCentMax(100.),
  90   fJetEtaMin(-.5),
  91   fJetEtaMax(.5),
  92   fJetPtMin(20.),
  93   fJetPtFractionMin(0.5),
  94   fNMatchJets(4),
  95   //fJetDeltaEta(.4),
  96   //fJetDeltaPhi(.4),
  97   fkNbranches(2),
  98   fkEvtClasses(5),
  99   fOutputList(0x0),
 100   fHistEvtSelection(0x0),
 101   fHistEvtClass(0x0),
 102   fHistCentrality(0x0),
 103   fHistPtJet(new TH1F*[fkNbranches]),
 104   fHistEtaPhiJet(new TH2F*[fkNbranches]),
 105   fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
 106   fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
 107   fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
 108   fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
 109   fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
 110   fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
 111   fHistPtFraction(new TH2F*[fkEvtClasses]),
 112   fHistPtPtExtra(0x0),
 113   fHistPtResponse(new TH2F*[fkEvtClasses]),
 114   fHistPtSmearing(new TH2F*[fkEvtClasses]),
 115   fHistDeltaR(new TH2F*[fkEvtClasses]),
 116   fHistArea(new TH2F*[fkEvtClasses]),
 117   fHistDeltaArea(new TH2F*[fkEvtClasses])
 118 {
 119   // Constructor
 120
 121   fJetBranchName[0] = "";
 122   fJetBranchName[1] = "";
 123
 124   fListJets[0] = new TList;
 125   fListJets[1] = new TList;
 126
 127   DefineOutput(1, TList::Class());
 128 }
 129
 130 AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
 131 {
 132   delete fListJets[0];
 133   delete fListJets[1];
 134 }
 135
 136 void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
 137 {
 138   fJetBranchName[0] = branch1;
 139   fJetBranchName[1] = branch2;
 140 }
 141
 142 void AliAnalysisTaskJetResponse::Init()
 143 {
 144
 145    // check for jet branches
 146    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 147       AliError("Jet branch name not set.");
 148    }
 149
 150 }
 151
 152 void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
 153 {
 154   // Create histograms
 155   // Called once
 156   OpenFile(1);
 157   if(!fOutputList) fOutputList = new TList;
 158   fOutputList->SetOwner(kTRUE);
 159
 160   Bool_t oldStatus = TH1::AddDirectoryStatus();
 161   TH1::AddDirectory(kFALSE);
 162
 163
 164   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 5, -0.5, 4.5);
 165   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 166   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 167   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 168   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 169   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 170
 171   fHistEvtClass   = new TH1I("fHistEvtClass", "event classes", fkEvtClasses, -0.5, fkEvtClasses-0.5);
 172   fHistCentrality = new TH1F("fHistCentrality", "event centrality", 100, 0., 100.);
 173
 174   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
 175     fHistPtJet[iBranch] = new TH1F(Form("fHistPtJet%i", iBranch),
 176                                           Form("p_{T} of jets from branch %i;p_{T} (GeV/c);counts", iBranch),
 177                                           250, 0., 250.);
 178     fHistPtJet[iBranch]->SetMarkerStyle(kFullCircle);
 179
 180     fHistEtaPhiJet[iBranch]    = new TH2F(Form("fHistEtaPhiJet%i", iBranch),
 181                                                          Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
 182                                                          48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
 183         fHistEtaPhiJetCut[iBranch] = new TH2F(Form("fHistEtaPhiJetCut%i", iBranch),
 184                                                          Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
 185                                                          48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
 186
 187   }
 188
 189
 190   fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
 191                             250, 0., 250., 250, 0., 250.);
 192
 193   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
 194     fHistDeltaEtaDeltaPhiJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJet%i",iEvtClass),
 195                                                                       "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
 196                                                                       101, -1.01, 1.01, 101, -1.01, 1.01);
 197     fHistDeltaEtaDeltaPhiJetCut[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJetCut%i",iEvtClass),
 198                                                                       "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
 199                                                                       101, -1.01, 1.01, 101, -1.01, 1.01);
 200     fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass] = new TH2F("fHistDeltaEtaDeltaPhiJetNOMatching",
 201                                                                                  "#Delta#eta - #Delta#phi of jets which do not match;#Delta#eta;#Delta#phi",
 202                                                                                  100, -2., 2., 100, TMath::Pi(), TMath::Pi());
 203     fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
 204                                                               250, 0., 250., 250, 0., 250.);
 205     fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;#Deltap_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
 206                                                               200, -50., 150., 250, 0., 250.);
 207     fHistDeltaR[iEvtClass]     = new TH2F(Form("hist_DeltaR%i",iEvtClass), "#DeltaR of matched jets;#Deltap_{T};#DeltaR", 200, -50.,150., 60, 0.,.6);
 208     fHistArea[iEvtClass]       = new TH2F(Form("hist_Area%i",iEvtClass), "jet area;#Deltap_{T};jet area", 200, -50.,150., 100, 0.,1.0);
 209     fHistDeltaArea[iEvtClass]  = new TH2F(Form("hist_DeltaArea%i",iEvtClass), "#DeltaArea of matched jets;#Deltap_{T};#Deltaarea", 200, -50., 150., 60, 0., .3);
 210
 211     fHistDeltaEtaEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaJet%i", iEvtClass),
 212                                                 Form("#eta - #Delta#eta of matched jets from event class %i;#eta;#Delta#eta", iEvtClass),
 213                                                 60, -.6, .6, 100, -.5, .5);
 214     fHistDeltaPtEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaJet%i", iEvtClass),
 215                                                 Form("#eta - #Deltap_{T} of matched jets from event class %i;#eta;#Deltap_{T}", iEvtClass),
 216                                                 60, -.6, .6, 200, -50., 150);
 217     fHistPtFraction[iEvtClass] = new TH2F(Form("fHistPtFraction%i", iEvtClass),
 218                                               Form("fraction from embedded jet in reconstructed jet;fraction (event class %i);p_{T}^{emb}", iEvtClass),
 219                                                                                   100, 0., 1., 250, 0., 250.);
 220   }
 221
 222
 223   fOutputList->Add(fHistEvtSelection);
 224   fOutputList->Add(fHistEvtClass);
 225   fOutputList->Add(fHistCentrality);
 226
 227   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
 228     fOutputList->Add(fHistPtJet[iBranch]);
 229     fOutputList->Add(fHistEtaPhiJet[iBranch]);
 230         fOutputList->Add(fHistEtaPhiJetCut[iBranch]);
 231   }
 232
 233   fOutputList->Add(fHistPtPtExtra);
 234
 235   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
 236     fOutputList->Add(fHistDeltaEtaDeltaPhiJet[iEvtClass]);
 237     fOutputList->Add(fHistDeltaEtaDeltaPhiJetCut[iEvtClass]);
 238     fOutputList->Add(fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass]);
 239         fOutputList->Add(fHistDeltaEtaEtaJet[iEvtClass]);
 240     fOutputList->Add(fHistDeltaPtEtaJet[iEvtClass]);
 241         fOutputList->Add(fHistPtFraction[iEvtClass]);
 242
 243     fOutputList->Add(fHistPtResponse[iEvtClass]);
 244     fOutputList->Add(fHistPtSmearing[iEvtClass]);
 245     fOutputList->Add(fHistDeltaR[iEvtClass]);
 246     fOutputList->Add(fHistArea[iEvtClass]);
 247     fOutputList->Add(fHistDeltaArea[iEvtClass]);
 248   }
 249
 250   // =========== Switch on Sumw2 for all histos ===========
 251   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 252     TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 253     if (h1){
 254       h1->Sumw2();
 255       continue;
 256     }
 257   }
 258   TH1::AddDirectory(oldStatus);
 259
 260   PostData(1, fOutputList);
 261 }
 262
 263 void AliAnalysisTaskJetResponse::UserExec(Option_t *)
 264 {
 265   // load events, apply event cuts, then compare leading jets
 266
 267   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 268       AliError("Jet branch name not set.");
 269       return;
 270    }
 271
 272   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 273   if (!fESD) {
 274     AliError("ESD not available");
 275     return;
 276   }
 277   fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
 278   if (!fAOD) {
 279     AliError("AOD not available");
 280     return;
 281   }
 282
 283   // -- event selection --
 284   fHistEvtSelection->Fill(1); // number of events before event selection
 285
 286   // physics selection
 287   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 288     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 289   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 290     if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 291     fHistEvtSelection->Fill(2);
 292     PostData(1, fOutputList);
 293     return;
 294   }
 295
 296   // vertex selection
 297   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 298   Int_t nTracksPrim = primVtx->GetNContributors();
 299   if ((nTracksPrim < fMinContribVtx) ||
 300       (primVtx->GetZ() < fVtxZMin) ||
 301       (primVtx->GetZ() > fVtxZMax) ){
 302     if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 303     fHistEvtSelection->Fill(3);
 304     PostData(1, fOutputList);
 305     return;
 306   }
 307
 308   // event class selection (from jet helper task)
 309   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 310   if(fDebug) Printf("Event class %d", eventClass);
 311   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 312     fHistEvtSelection->Fill(4);
 313     PostData(1, fOutputList);
 314     return;
 315   }
 316
 317   // centrality selection
 318   AliCentrality *cent = fESD->GetCentrality();
 319   Float_t centValue = cent->GetCentralityPercentile("TRK");
 320   if(fDebug) printf("centrality: %f\n", centValue);
 321   if (centValue < fCentMin || centValue > fCentMax){
 322     fHistEvtSelection->Fill(4);
 323     PostData(1, fOutputList);
 324     return;
 325   }
 326
 327   fHistEvtSelection->Fill(0); // accepted events
 328   fHistEvtClass->Fill(eventClass);
 329   fHistCentrality->Fill(centValue);
 330   // -- end event selection --
 331
 332   // fetch jets
 333   TClonesArray *aodJets[2];
 334   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
 335   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
 336
 337   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 338     fListJets[iJetType]->Clear();
 339     if (!aodJets[iJetType]) continue;
 340
 341     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 342       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 343       if (jet) fListJets[iJetType]->Add(jet);
 344     }
 345     fListJets[iJetType]->Sort();
 346   }
 347
 348   // jet matching
 349   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
 350   static TArrayF aPtFraction(fListJets[0]->GetEntries());
 351   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
 352   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
 353
 354   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
 355   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
 356                                             fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
 357                                             aMatchIndex, aPtFraction, fDebug);
 358
 359   // loop over matched jets
 360   Int_t ir = -1; // index of matched reconstruced jet
 361   Float_t fraction = 0.;
 362   AliAODJet *jet[2]  = { 0x0, 0x0 };
 363   Float_t jetEta[2]  = { 0., 0. };
 364   Float_t jetPhi[2]  = { 0., 0. };
 365   Float_t jetPt[2]   = { 0., 0. };
 366   Float_t jetArea[2] = { 0., 0. };
 367
 368   for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
 369      ir = aMatchIndex[ig];
 370          if(ir<0) continue;
 371          fraction = aPtFraction[ig];
 372
 373          // fetch jets
 374          jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
 375          jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
 376          if(!jet[0] || !jet[1]) continue;
 377
 378          for(Int_t i=0; i<fkNbranches; ++i){
 379             jetEta[i]  = jet[i]->Eta();
 380                 jetPhi[i]  = jet[i]->Phi();
 381                 jetPt[i]   = jet[i]->Pt();
 382                 jetArea[i] = jet[i]->EffectiveAreaCharged();
 383          }
 384          if(eventClass > -1 && eventClass < fkEvtClasses){
 385         fHistPtFraction[eventClass] -> Fill(fraction, jetPt[0]);
 386          }
 387
 388          if(fraction<fJetPtFractionMin) continue;
 389
 390          // calculate parameters of associated jets
 391          Float_t deltaPt    = jetPt[1]-jetPt[0];
 392          Float_t deltaEta   = jetEta[1]-jetEta[0];
 393          Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
 394          Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
 395          Float_t deltaArea  = jetArea[1]-jetArea[0];
 396
 397          // fill histograms before acceptance cut
 398          for(Int_t i=0; i<fkNbranches; ++i){
 399              fHistEtaPhiJet[i] -> Fill(jetEta[i], jetPhi[i]);
 400          }
 401          if(eventClass > -1 && eventClass < fkEvtClasses){
 402              fHistDeltaEtaDeltaPhiJet[eventClass] -> Fill(deltaEta, deltaPhi);
 403          }
 404
 405          // jet acceptance + minimum pT check
 406          if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
 407             jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
 408
 409                 if(fDebug){
 410                 Printf("Jet not in eta acceptance.");
 411                         Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
 412                         Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
 413                 }
 414                 continue;
 415      }
 416          if(jetPt[1] < fJetPtMin){
 417             if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
 418             continue;
 419          }
 420
 421
 422
 423          // fill histograms
 424          for(Int_t i=0; i<fkNbranches; ++i){
 425              fHistPtJet[i]        -> Fill(jetPt[i]);
 426              fHistEtaPhiJetCut[i] -> Fill(jetEta[i], jetPhi[i]);
 427          }
 428
 429          fHistPtPtExtra->Fill(jetPt[0], jetPt[1]);
 430
 431          if(eventClass > -1 && eventClass < fkEvtClasses){
 432                  fHistDeltaEtaDeltaPhiJetCut[eventClass] -> Fill(deltaEta, deltaPhi);
 433
 434                  fHistPtResponse[eventClass]             -> Fill(jetPt[1], jetPt[0]);
 435                  fHistPtSmearing[eventClass]             -> Fill(deltaPt,  jetPt[0]);
 436
 437                  fHistDeltaPtEtaJet[eventClass]          -> Fill(jetEta[0], deltaPt);
 438                  fHistDeltaEtaEtaJet[eventClass]         -> Fill(jetEta[0], deltaEta);
 439
 440                  fHistDeltaR[eventClass]                 -> Fill(deltaPt, deltaR);
 441                  fHistArea[eventClass]                   -> Fill(deltaPt, jetArea[1]);
 442                  fHistDeltaArea[eventClass]              -> Fill(deltaPt, deltaArea);
 443
 444          }
 445   }
 446
 447   PostData(1, fOutputList);
 448 }
 449
 450 void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
 451 {
 452   // Draw result to the screen
 453   // Called once at the end of the query
 454
 455   if (!GetOutputData(1))
 456     return;
 457 }
 458