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