PWG4/JetTasks/AliAnalysisTaskJFSystematics.cxx

   1 // **************************************
   2 // Task used for the systematic study of jet finders
   3 //
   4 // Compares input (gen) and output (rec) jets
   5 // gen can also be another jet finder on the rec level, matching is done in eta phi
   6 //
   7 // *******************************************
   8
   9
  10 /**************************************************************************
  11  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  12  *                                                                        *
  13  * Author: The ALICE Off-line Project.                                    *
  14  * Contributors are mentioned in the code where appropriate.              *
  15  *                                                                        *
  16  * Permission to use, copy, modify and distribute this software and its   *
  17  * documentation strictly for non-commercial purposes is hereby granted   *
  18  * without fee, provided that the above copyright notice appears in all   *
  19  * copies and that both the copyright notice and this permission notice   *
  20  * appear in the supporting documentation. The authors make no claims     *
  21  * about the suitability of this software for any purpose. It is          *
  22  * provided "as is" without express or implied warranty.                  *
  23  **************************************************************************/
  24
  25
  26 #include <TROOT.h>
  27 #include <TRandom.h>
  28 #include <TSystem.h>
  29 #include <TInterpreter.h>
  30 #include <TChain.h>
  31 #include <TFile.h>
  32 #include <TH1F.h>
  33 #include <TH2F.h>
  34 #include <TH3F.h>
  35 #include <TProfile.h>
  36 #include <TList.h>
  37 #include <TLorentzVector.h>
  38 #include <TClonesArray.h>
  39 #include  "TDatabasePDG.h"
  40
  41 #include "AliAnalysisTaskJFSystematics.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliJetFinder.h"
  44 #include "AliJetHeader.h"
  45 #include "AliJetReader.h"
  46 #include "AliJetReaderHeader.h"
  47 #include "AliUA1JetHeaderV1.h"
  48 #include "AliESDEvent.h"
  49 #include "AliAODEvent.h"
  50 #include "AliAODHandler.h"
  51 #include "AliAODTrack.h"
  52 #include "AliAODJet.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliStack.h"
  56 #include "AliGenPythiaEventHeader.h"
  57 #include "AliJetKineReaderHeader.h"
  58 #include "AliGenCocktailEventHeader.h"
  59 #include "AliInputEventHandler.h"
  60
  61
  62 #include "AliAnalysisHelperJetTasks.h"
  63
  64 ClassImp(AliAnalysisTaskJFSystematics)
  65
  66   const Int_t AliAnalysisTaskJFSystematics::fgkSysBins[AliAnalysisTaskJFSystematics::kSysTypes] = {0,AliAnalysisTaskJFSystematics::kMaxJets};
  67 const char* AliAnalysisTaskJFSystematics::fgkSysName[AliAnalysisTaskJFSystematics::kSysTypes] = {"","j"};
  68
  69 AliAnalysisTaskJFSystematics::AliAnalysisTaskJFSystematics(): AliAnalysisTaskSE(),
  70   fJetHeaderRec(0x0),
  71   fJetHeaderGen(0x0),
  72   fAOD(0x0),
  73   fBranchRec("jets"),
  74   fBranchGen(""),
  75   fUseAODInput(kFALSE),
  76   fUseExternalWeightOnly(kFALSE),
  77   fLimitGenJetEta(kFALSE),
  78   fAnalysisType(0),
  79   fExternalWeight(1),
  80   fRecEtaWindow(0.5),
  81   fh1Xsec(0x0),
  82   fh1Trials(0x0),
  83   fh1PtHard(0x0),
  84   fh1PtHardNoW(0x0),
  85   fh1PtHardTrials(0x0),
  86   fh1NGenJets(0x0),
  87   fh1NRecJets(0x0),
  88   fh1PtRecIn(0x0),
  89   fh1PtRecOut(0x0),
  90   fh1PtGenIn(0x0),
  91   fh1PtGenOut(0x0),
  92   fh2PtFGen(0x0),
  93   fh2PhiFGen(0x0),
  94   fh2EtaFGen(0x0),
  95   fh2PtGenDeltaPhi(0x0),
  96   fh2PtGenDeltaEta(0x0),
  97   fh3RecInEtaPhiPt(0x0),
  98   fh3RecOutEtaPhiPt(0x0),
  99   fh3GenInEtaPhiPt(0x0),
 100   fh3GenOutEtaPhiPt(0x0),
 101   fhnCorrelation(0x0),
 102   fHistList(0x0)
 103 {
 104   // Default constructor
 105   /*
 106     for(int ij  = 0;ij<kMaxJets;++ij){
 107       fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
 108       fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] = 0;
 109       fh3RecInEtaPhiPt[ij] = fh3RecOutEtaPhiPt[ij] =fh3GenInEtaPhiPt[ij] =  fh3GenOutEtaPhiPt[ij] = 0;
 110       fhnCorrelation[ij] = 0;
 111     }
 112   */
 113 }
 114
 115 AliAnalysisTaskJFSystematics::AliAnalysisTaskJFSystematics(const char* name):
 116   AliAnalysisTaskSE(name),
 117   fJetHeaderRec(0x0),
 118   fJetHeaderGen(0x0),
 119   fAOD(0x0),
 120   fBranchRec("jets"),
 121   fBranchGen(""),
 122   fUseAODInput(kFALSE),
 123   fUseExternalWeightOnly(kFALSE),
 124   fLimitGenJetEta(kFALSE),
 125   fAnalysisType(0),
 126   fExternalWeight(1),
 127   fRecEtaWindow(0.5),
 128   fh1Xsec(0x0),
 129   fh1Trials(0x0),
 130   fh1PtHard(0x0),
 131   fh1PtHardNoW(0x0),
 132   fh1PtHardTrials(0x0),
 133   fh1NGenJets(0x0),
 134   fh1NRecJets(0x0),
 135   fh1PtRecIn(0x0),
 136   fh1PtRecOut(0x0),
 137   fh1PtGenIn(0x0),
 138   fh1PtGenOut(0x0),
 139   fh2PtFGen(0x0),
 140   fh2PhiFGen(0x0),
 141   fh2EtaFGen(0x0),
 142   fh2PtGenDeltaPhi(0x0),
 143   fh2PtGenDeltaEta(0x0),
 144   fh3RecInEtaPhiPt(0x0),
 145   fh3RecOutEtaPhiPt(0x0),
 146   fh3GenInEtaPhiPt(0x0),
 147   fh3GenOutEtaPhiPt(0x0),
 148   fhnCorrelation(0x0),
 149   fHistList(0x0)
 150 {
 151   // Default constructor
 152
 153   // Default constructor
 154   /*
 155   for(int ij  = 0;ij<kMaxJets;++ij){
 156     fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
 157     fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] = 0;
 158     fh3RecInEtaPhiPt[ij] = fh3RecOutEtaPhiPt[ij] =fh3GenInEtaPhiPt[ij] =  fh3GenOutEtaPhiPt[ij] = 0;
 159     fhnCorrelation[ij] = 0;
 160   }
 161   */
 162   DefineOutput(1, TList::Class());
 163 }
 164
 165
 166
 167 Bool_t AliAnalysisTaskJFSystematics::Notify()
 168 {
 169
 170   //
 171   // Implemented Notify() to read the cross sections
 172   // and number of trials from pyxsec.root
 173   //
 174   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 175   Double_t xsection = 0;
 176   UInt_t   ntrials  = 0;
 177   if(tree){
 178     TFile *curfile = tree->GetCurrentFile();
 179     if (!curfile) {
 180       Error("Notify","No current file");
 181       return kFALSE;
 182     }
 183     if(!fh1Xsec||!fh1Trials){
 184       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 185       return kFALSE;
 186     }
 187
 188     TString fileName(curfile->GetName());
 189     if(fileName.Contains("AliESDs.root")){
 190         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
 191     }
 192     else if(fileName.Contains("AliAOD.root")){
 193         fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
 194     }
 195     else if(fileName.Contains("galice.root")){
 196         // for running with galice and kinematics alone...
 197         fileName.ReplaceAll("galice.root", "pyxsec.root");
 198     }
 199
 200     TFile *fxsec = TFile::Open(fileName.Data());
 201     if(!fxsec){
 202       Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
 203       // no a severe condition, will happen for minbias data
 204       return kTRUE;
 205     }
 206     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 207     if(!xtree){
 208       Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
 209     }
 210     xtree->SetBranchAddress("xsection",&xsection);
 211     xtree->SetBranchAddress("ntrials",&ntrials);
 212     xtree->GetEntry(0);
 213     fh1Xsec->Fill("<#sigma>",xsection);
 214     fh1Trials->Fill("#sum{ntrials}",ntrials);
 215   }
 216   return kTRUE;
 217 }
 218
 219 void AliAnalysisTaskJFSystematics::UserCreateOutputObjects()
 220 {
 221
 222   //
 223   // Create the output container
 224   //
 225
 226
 227   // Connect the AOD
 228
 229   if(fUseAODInput){
 230     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 231     if(!fAOD){
 232       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 233       return;
 234     }
 235     // fetch the header
 236     fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 237     if(!fJetHeaderRec){
 238       Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
 239     }
 240   }
 241   else{
 242     //  assume that the AOD is in the general output...
 243     fAOD  = AODEvent();
 244     if(!fAOD){
 245       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 246       return;
 247     }
 248     fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 249     if(!fJetHeaderRec){
 250       Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
 251     }
 252     else{
 253       if(fDebug>10)fJetHeaderRec->Dump();
 254     }
 255   }
 256
 257
 258
 259   if (fDebug > 1) printf("AnalysisTaskJFSystematics::UserCreateOutputObjects() \n");
 260
 261   OpenFile(1);
 262   if(!fHistList)fHistList = new TList();
 263
 264   Bool_t oldStatus = TH1::AddDirectoryStatus();
 265   TH1::AddDirectory(kFALSE);
 266
 267   //
 268   //  Histograms
 269   //
 270
 271   const Int_t nBinPt = 100;
 272   Double_t binLimitsPt[nBinPt+1];
 273   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 274     if(iPt == 0){
 275       binLimitsPt[iPt] = 0.0;
 276     }
 277     else {// 1.0
 278       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
 279     }
 280   }
 281
 282   const Int_t nBinEta = 26;
 283   Double_t binLimitsEta[nBinEta+1] = {
 284     -1.6,-1.4,-1.2,-1.0,
 285     -0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
 286     0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
 287     1.0, 1.2, 1.4, 1.6
 288   };
 289
 290
 291   const Int_t nBinPhi = 18;
 292   Double_t binLimitsPhi[nBinPhi+1];
 293   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 294     if(iPhi==0){
 295       binLimitsPhi[iPhi] = 0;
 296     }
 297     else{
 298       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 299     }
 300   }
 301
 302
 303   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 304   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 305
 306   fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
 307   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 308
 309   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 310   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 311   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 312   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 313   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 314
 315   // book the single histograms these we clone for various systematics
 316   //
 317    fh1PtRecIn =  new TH1F("fh1PtRecIn","rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 318   fh1PtRecOut = new TH1F("fh1PtRecOut","rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
 319   fh1PtGenIn = new TH1F("fh1PtGenIn","found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 320   fh1PtGenOut = new TH1F("fh1PtGenOut","found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
 321
 322
 323
 324   fh2PtFGen = new TH2F("fh2PtFGen","Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
 325                              nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 326   fh2PhiFGen = new TH2F("fh2PhiFGen","#phi Found vs. gen;#phi_{rec};#phi_{gen}",
 327                         nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
 328   fh2EtaFGen = new TH2F("fh2EtaFGen","#eta Found vs. gen;#eta_{rec};#eta_{gen}",
 329                         nBinEta,binLimitsEta,nBinEta,binLimitsEta);
 330
 331   fh2PtGenDeltaPhi = new TH2F("fh2PtGenDeltaPhi","delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-#phi_{rec}",
 332                               nBinPt,binLimitsPt,100,-1.0,1.0);
 333   fh2PtGenDeltaEta = new TH2F("fh2PtGenDeltaEta_j%d","delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-#eta_{rec}",
 334                               nBinPt,binLimitsPt,100,-1.0,1.0);
 335
 336
 337   fh3RecInEtaPhiPt = new TH3F("fh3RecInEtaPhiPt","Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 338                               nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 339   fh3RecOutEtaPhiPt = new TH3F("fh3RecOutEtaPhiPt","generated found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 340                                    nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 341   fh3GenInEtaPhiPt = new TH3F("fh3GenInEtaPhiPt","generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
 342                                   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 343   fh3GenOutEtaPhiPt = new TH3F("fh3GenOutEtaPhiPt","reconstructed found for Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
 344                                    nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 345
 346   const int nbin[4] = {nBinPt,nBinPt,24,24};
 347   Double_t vLowEdge[4] = {0,0,-1.2,-1.2};
 348   Double_t vUpEdge[4] = {250,250,1.2,1.2};
 349
 350   fhnCorrelation  = new THnSparseF("fhnCorrelation","Response Map", 4, nbin, vLowEdge, vUpEdge);
 351
 352
 353
 354   fHistList->Add(fh1Xsec);
 355   fHistList->Add(fh1Trials);
 356   fHistList->Add(fh1PtHard);
 357   fHistList->Add(fh1PtHardNoW);
 358   fHistList->Add(fh1PtHardTrials);
 359   fHistList->Add(fh1NGenJets);
 360   fHistList->Add(fh1NRecJets);
 361   fHistList->Add(fh1PtRecIn);
 362   fHistList->Add(fh1PtRecOut);
 363   fHistList->Add(fh1PtGenIn);
 364   fHistList->Add(fh1PtGenOut);
 365   fHistList->Add(fh2PtFGen);
 366   fHistList->Add(fh2PhiFGen);
 367   fHistList->Add(fh2EtaFGen);
 368   fHistList->Add(fh2PtGenDeltaEta);
 369   fHistList->Add(fh2PtGenDeltaPhi);
 370   fHistList->Add(fh3RecOutEtaPhiPt);
 371   fHistList->Add(fh3GenOutEtaPhiPt);
 372   fHistList->Add(fh3RecInEtaPhiPt);
 373   fHistList->Add(fh3GenInEtaPhiPt);
 374   fHistList->Add(fhnCorrelation);
 375
 376
 377   if(fAnalysisType==kSysJetOrder){
 378     //
 379     for(int i = 0; i< fgkSysBins[kSysJetOrder];++i){
 380       TH1F *hTmp = (TH1F*)fh1PtRecIn->Clone(Form("%s_%s%d",fh1PtRecIn->GetName(),fgkSysName[kSysJetOrder],i));
 381       fHistList->Add(hTmp);
 382       hTmp = (TH1F*)fh1PtRecOut->Clone(Form("%s_%s%d",fh1PtRecOut->GetName(),fgkSysName[kSysJetOrder],i));
 383       fHistList->Add(hTmp);
 384       hTmp = (TH1F*)fh1PtGenIn->Clone(Form("%s_%s%d",fh1PtGenIn->GetName(),fgkSysName[kSysJetOrder],i));
 385       fHistList->Add(hTmp);
 386       hTmp = (TH1F*)fh1PtGenOut->Clone(Form("%s_%s%d",fh1PtGenOut->GetName(),fgkSysName[kSysJetOrder],i));
 387       fHistList->Add(hTmp);
 388       THnSparseF *hnTmp = (THnSparseF*)fhnCorrelation->Clone(Form("%s_%s%d",fhnCorrelation->GetName(),fgkSysName[kSysJetOrder],i));
 389       fHistList->Add(hnTmp);
 390     }
 391   }
 392
 393   // =========== Switch on Sumw2 for all histos ===========
 394   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 395     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 396     if (h1){
 397       // Printf("%s ",h1->GetName());
 398       h1->Sumw2();
 399       continue;
 400     }
 401     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 402     if(hn)hn->Sumw2();
 403   }
 404
 405   TH1::AddDirectory(oldStatus);
 406
 407   }
 408
 409 void AliAnalysisTaskJFSystematics::Init()
 410 {
 411   //
 412   // Initialization
 413   //
 414
 415   Printf(">>> AnalysisTaskJFSystematics::Init() debug level %d\n",fDebug);
 416   if (fDebug > 1) printf("AnalysisTaskJFSystematics::Init() \n");
 417
 418 }
 419
 420 void AliAnalysisTaskJFSystematics::UserExec(Option_t */*option*/)
 421 {
 422   //
 423   // Execute analysis for current event
 424   //
 425
 426
 427
 428   if (fDebug > 1)printf("AliAnalysisTaskJFSystematics::Analysing event # %5d\n", (Int_t) fEntry);
 429
 430   // ========= These pointers need to be valid in any case =======
 431
 432   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 433   if(!aodRecJets){
 434     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 435     return;
 436   }
 437
 438   // We use static arrays, not to fragment the memory
 439   //
 440   AliAODJet genJets[kMaxJets];
 441   Int_t nGenJets = 0;
 442   AliAODJet recJets[kMaxJets];
 443   Int_t nRecJets = 0;
 444
 445   Double_t eventW = 1;
 446   Double_t ptHard = 0;
 447
 448   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
 449   Int_t     iProcessType = 0;
 450   if(fUseExternalWeightOnly){
 451     eventW = fExternalWeight;
 452   }
 453
 454   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 455   // this is the part where when we need to have MC information
 456   // we can also work on Reconstructed only when just comparing two JF
 457   AliMCEvent* mcEvent = MCEvent();
 458   if(!mcEvent){
 459     Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 460   }
 461   else {
 462     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 463     if(pythiaGenHeader){
 464       nTrials = pythiaGenHeader->Trials();
 465       ptHard  = pythiaGenHeader->GetPtHard();
 466       iProcessType = pythiaGenHeader->ProcessType();
 467       // 11 f+f -> f+f
 468       // 12 f+barf -> f+barf
 469       // 13 f+barf -> g+g
 470       // 28 f+g -> f+g
 471       // 53 g+g -> f+barf
 472       // 68 g+g -> g+g
 473       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 474       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 475       // fetch the pythia generated jets only to be used here
 476       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 477       AliAODJet pythiaGenJets[kMaxJets];
 478       Int_t iCount = 0;
 479       for(int ip = 0;ip < nPythiaGenJets;++ip){
 480         if(iCount>=kMaxJets)continue;
 481         Float_t p[4];
 482         pythiaGenHeader->TriggerJet(ip,p);
 483         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 484         if(fLimitGenJetEta){
 485           if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 486              pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 487       }
 488         if(fBranchGen.Length()==0){
 489           // if we have MC particles and we do not read from the aod branch
 490           // use the pythia jets
 491           genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 492         }
 493       iCount++;
 494       }
 495       if(fBranchGen.Length()==0)nGenJets = iCount;
 496     }
 497   }// if we had the MCEvent
 498
 499   fh1PtHard->Fill(ptHard,eventW);
 500   fh1PtHardNoW->Fill(ptHard,1);
 501   fh1PtHardTrials->Fill(ptHard,nTrials);
 502
 503   // If we set a second branch for the input jets fetch this
 504   if(fBranchGen.Length()>0){
 505     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 506     if(aodGenJets){
 507       Int_t iCount = 0;
 508       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 509         if(iCount>=kMaxJets)continue;
 510         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 511         if(!tmp)continue;
 512         if(fLimitGenJetEta){
 513           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 514              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 515         }
 516         genJets[iCount] = *tmp;
 517         iCount++;
 518       }
 519       nGenJets = iCount;
 520     }
 521     else{
 522       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 523     }
 524   }
 525
 526   fh1NGenJets->Fill(nGenJets);
 527   // We do not want to exceed the maximum jet number
 528   nGenJets = TMath::Min(nGenJets,kMaxJets);
 529
 530
 531   //
 532   // Fetch the reconstructed jets...
 533   //
 534
 535   nRecJets = aodRecJets->GetEntries();
 536   fh1NRecJets->Fill(nRecJets);
 537   nRecJets = TMath::Min(nRecJets,kMaxJets);
 538
 539   for(int ir = 0;ir < nRecJets;++ir){
 540     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 541     if(!tmp)continue;
 542     recJets[ir] = *tmp;
 543   }
 544
 545
 546   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 547
 548   //
 549   // Relate the jets
 550   //
 551   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 552   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 553
 554   for(int i = 0;i<kMaxJets;++i){
 555     iGenIndex[i] = iRecIndex[i] = -1;
 556   }
 557
 558
 559   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 560                                             iGenIndex,iRecIndex,fDebug);
 561   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 562
 563   if(fDebug){
 564     for(int i = 0;i<kMaxJets;++i){
 565       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 566       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 567     }
 568   }
 569
 570   //
 571   // Now the premliminaries are over, lets do the jet analysis
 572   //
 573
 574
 575   Double_t value[4]; // for the thnsparse
 576   // loop over reconstructed jets
 577   for(int ir = 0;ir < nRecJets;++ir){
 578     Double_t ptRec = recJets[ir].Pt();
 579     Double_t phiRec = recJets[ir].Phi();
 580     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 581     Double_t etaRec = recJets[ir].Eta();
 582     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 583     fh1PtRecIn->Fill(ptRec,eventW);
 584     if(fAnalysisType==kSysJetOrder)((TH1F*)fHistList->FindObject(Form("fh1PtRecIn_%s%d",fgkSysName[kSysJetOrder],ir)))->Fill(ptRec,eventW);
 585     fh3RecInEtaPhiPt->Fill(etaRec,phiRec,ptRec,eventW);
 586     // Fill Correlation
 587     Int_t ig = iGenIndex[ir];
 588     if(ig>=0 && ig<nGenJets){
 589       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 590       fh1PtRecOut->Fill(ptRec,eventW);
 591       if(fAnalysisType==kSysJetOrder)((TH1F*)fHistList->FindObject(Form("fh1PtRecOut_%s%d",fgkSysName[kSysJetOrder],ir)))->Fill(ptRec,eventW);
 592       Double_t ptGen  = genJets[ig].Pt();
 593       Double_t phiGen = genJets[ig].Phi();
 594       if(phiGen<0)phiGen+=TMath::Pi()*2.;
 595       Double_t etaGen = genJets[ig].Eta();
 596
 597       fh3RecOutEtaPhiPt->Fill(etaRec,phiRec,ptRec,eventW);
 598
 599       value[0] = ptRec;
 600       value[1] = ptGen;
 601       value[2] = etaRec;
 602       value[3] = etaGen;
 603
 604       fhnCorrelation->Fill(value,eventW);
 605       if(fAnalysisType==kSysJetOrder)((THnSparseF*)fHistList->FindObject(Form("fhnCorrelation_%s%d",fgkSysName[kSysJetOrder],ir)))->Fill(value,eventW);
 606       //
 607       // we accept only jets which are detected within a smaller window, to
 608       // avoid ambigious pair association at the edges of the acceptance
 609       //
 610
 611       if(TMath::Abs(etaRec)<fRecEtaWindow){
 612         fh2PtFGen->Fill(ptRec,ptGen,eventW);
 613         fh2PhiFGen->Fill(phiRec,phiGen,eventW);
 614         fh2EtaFGen->Fill(etaRec,etaGen,eventW);
 615         fh2PtGenDeltaEta->Fill(ptGen,etaGen-etaRec,eventW);
 616         fh2PtGenDeltaPhi->Fill(ptGen,phiGen-phiRec,eventW);
 617       }// if etarec in window
 618     }//if ig valid
 619   }// loop over reconstructed jets
 620
 621
 622   // Now llop over generated jets
 623
 624   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 625   for(int ig = 0;ig < nGenJets;++ig){
 626     Double_t ptGen = genJets[ig].Pt();
 627     // Fill Correlation
 628     Double_t phiGen = genJets[ig].Phi();
 629     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 630     Double_t etaGen = genJets[ig].Eta();
 631     fh3GenInEtaPhiPt->Fill(etaGen,phiGen,ptGen,eventW);
 632     fh1PtGenIn->Fill(ptGen,eventW);
 633     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 634     if(fAnalysisType==kSysJetOrder)((TH1F*)fHistList->FindObject(Form("fh1PtGenIn_%s%d",fgkSysName[kSysJetOrder],ig)))->Fill(ptGen,eventW);
 635     Int_t ir = iRecIndex[ig];
 636     if(ir>=0&&ir<nRecJets){
 637       fh1PtGenOut->Fill(ptGen,eventW);
 638       fh3GenOutEtaPhiPt->Fill(etaGen,phiGen,ptGen,eventW);
 639       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 640       if(fAnalysisType==kSysJetOrder)((TH1F*)fHistList->FindObject(Form("fh1PtGenOut_%s%d",fgkSysName[kSysJetOrder],ig)))->Fill(ptGen,eventW);
 641     }
 642   }// loop over reconstructed jets
 643
 644   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 645   PostData(1, fHistList);
 646 }
 647
 648 void AliAnalysisTaskJFSystematics::Terminate(Option_t */*option*/)
 649 {
 650 // Terminate analysis
 651 //
 652     if (fDebug > 1) printf("AnalysisTaskJFSystematics: Terminate() \n");
 653     // Plot
 654
 655
 656 }