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