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