PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx

   1 // **************************************
   2 // Task used for the correction of determiantion of reconstructed jet spectra
   3 // Compares input (gen) and output (rec) jets
   4 // *******************************************
   5
   6
   7 /**************************************************************************
   8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   9  *                                                                        *
  10  * Author: The ALICE Off-line Project.                                    *
  11  * Contributors are mentioned in the code where appropriate.              *
  12  *                                                                        *
  13  * Permission to use, copy, modify and distribute this software and its   *
  14  * documentation strictly for non-commercial purposes is hereby granted   *
  15  * without fee, provided that the above copyright notice appears in all   *
  16  * copies and that both the copyright notice and this permission notice   *
  17  * appear in the supporting documentation. The authors make no claims     *
  18  * about the suitability of this software for any purpose. It is          *
  19  * provided "as is" without express or implied warranty.                  *
  20  **************************************************************************/
  21
  22
  23 #include <TROOT.h>
  24 #include <TRandom.h>
  25 #include <TSystem.h>
  26 #include <TInterpreter.h>
  27 #include <TChain.h>
  28 #include <TFile.h>
  29 #include <TKey.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TH3F.h>
  33 #include <TProfile.h>
  34 #include <TList.h>
  35 #include <TLorentzVector.h>
  36 #include <TClonesArray.h>
  37 #include  "TDatabasePDG.h"
  38
  39 #include "AliAnalysisTaskJetSpectrum.h"
  40 #include "AliAnalysisManager.h"
  41 #include "AliJetFinder.h"
  42 #include "AliJetHeader.h"
  43 #include "AliJetReader.h"
  44 #include "AliJetReaderHeader.h"
  45 #include "AliUA1JetHeaderV1.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODHandler.h"
  49 #include "AliAODTrack.h"
  50 #include "AliAODJet.h"
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCEvent.h"
  53 #include "AliStack.h"
  54 #include "AliGenPythiaEventHeader.h"
  55 #include "AliJetKineReaderHeader.h"
  56 #include "AliGenCocktailEventHeader.h"
  57 #include "AliInputEventHandler.h"
  58
  59
  60 #include "AliAnalysisHelperJetTasks.h"
  61
  62 ClassImp(AliAnalysisTaskJetSpectrum)
  63
  64 const Float_t AliAnalysisTaskJetSpectrum::fgkJetNpartCut[AliAnalysisTaskJetSpectrum::kMaxCorrelation] = {5,10,1E+09};
  65
  66 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
  67   fJetHeaderRec(0x0),
  68   fJetHeaderGen(0x0),
  69   fAOD(0x0),
  70   fBranchRec("jets"),
  71   fBranchGen(""),
  72   fUseAODInput(kFALSE),
  73   fUseExternalWeightOnly(kFALSE),
  74   fLimitGenJetEta(kFALSE),
  75   fAnalysisType(0),
  76   fExternalWeight(1),
  77   fRecEtaWindow(0.5),
  78   fh1Xsec(0x0),
  79   fh1Trials(0x0),
  80   fh1PtHard(0x0),
  81   fh1PtHardNoW(0x0),
  82   fh1PtHardTrials(0x0),
  83   fh1NGenJets(0x0),
  84   fh1NRecJets(0x0),
  85   fHistList(0x0) ,
  86   ////////////////
  87   fh1JetMultiplicity(0x0) ,
  88   fh2ERecZRec(0x0) ,
  89   fh2EGenZGen(0x0) ,
  90   fh2Efficiency(0x0) ,
  91   fh3EGenERecN(0x0)
  92   ////////////////
  93 {
  94   // Default constructor
  95     for(int ij  = 0;ij<kMaxJets;++ij){
  96       fh1E[ij] =  fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
  97       fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] =   fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] =  fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] = 0;
  98       fh3PtRecGenHard[ij] =  fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
  99     }
 100     for(int ic = 0;ic < kMaxCorrelation;ic++){
 101       fhnCorrelation[ic] = 0;
 102     }
 103
 104 }
 105
 106 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
 107   AliAnalysisTaskSE(name),
 108   fJetHeaderRec(0x0),
 109   fJetHeaderGen(0x0),
 110   fAOD(0x0),
 111   fBranchRec("jets"),
 112   fBranchGen(""),
 113   fUseAODInput(kFALSE),
 114   fUseExternalWeightOnly(kFALSE),
 115   fLimitGenJetEta(kFALSE),
 116   fAnalysisType(0),
 117   fExternalWeight(1),
 118   fRecEtaWindow(0.5),
 119   fh1Xsec(0x0),
 120   fh1Trials(0x0),
 121   fh1PtHard(0x0),
 122   fh1PtHardNoW(0x0),
 123   fh1PtHardTrials(0x0),
 124   fh1NGenJets(0x0),
 125   fh1NRecJets(0x0),
 126   fHistList(0x0) ,
 127   ////////////////
 128   fh1JetMultiplicity(0x0) ,
 129   fh2ERecZRec(0x0) ,
 130   fh2EGenZGen(0x0) ,
 131   fh2Efficiency(0x0) ,
 132   fh3EGenERecN(0x0)
 133   ////////////////
 134 {
 135   // Default constructor
 136   for(int ij  = 0;ij<kMaxJets;++ij){
 137     fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
 138     fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] =   fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] =0;
 139
 140     fh3PtRecGenHard[ij] =  fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
 141   }
 142
 143     for(int ic = 0;ic < kMaxCorrelation;ic++){
 144       fhnCorrelation[ic] = 0;
 145     }
 146
 147   DefineOutput(1, TList::Class());
 148 }
 149
 150
 151
 152 Bool_t AliAnalysisTaskJetSpectrum::Notify()
 153 {
 154   //
 155   // Implemented Notify() to read the cross sections
 156   // and number of trials from pyxsec.root
 157   //
 158   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 159   Double_t xsection = 0;
 160   UInt_t   ntrials  = 0;
 161   Float_t   ftrials  = 0;
 162   if(tree){
 163     TFile *curfile = tree->GetCurrentFile();
 164     if (!curfile) {
 165       Error("Notify","No current file");
 166       return kFALSE;
 167     }
 168     if(!fh1Xsec||!fh1Trials){
 169       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 170       return kFALSE;
 171     }
 172
 173     TString fileName(curfile->GetName());
 174     if(fileName.Contains("AliESDs.root")){
 175         fileName.ReplaceAll("AliESDs.root", "");
 176     }
 177     else if(fileName.Contains("AliAOD.root")){
 178         fileName.ReplaceAll("AliAOD.root", "");
 179     }
 180     else if(fileName.Contains("AliAODs.root")){
 181         fileName.ReplaceAll("AliAODs.root", "");
 182     }
 183     else if(fileName.Contains("galice.root")){
 184         // for running with galice and kinematics alone...
 185         fileName.ReplaceAll("galice.root", "");
 186     }
 187     TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
 188     if(!fxsec){
 189       if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec.root"));
 190       // next trial fetch the histgram file
 191       fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
 192       if(!fxsec){
 193         // not a severe condition
 194         if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
 195         return kTRUE;
 196       }
 197       else{
 198         // find the tlist we want to be independtent of the name so use the Tkey
 199         TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 200         if(!key){
 201           if(fDebug>0)Printf("%s:%d key not found in the file",(char*)__FILE__,__LINE__);
 202           return kTRUE;
 203         }
 204         TList *list = dynamic_cast<TList*>(key->ReadObj());
 205         if(!list){
 206           if(fDebug>0)Printf("%s:%d key is not a tlist",(char*)__FILE__,__LINE__);
 207           return kTRUE;
 208         }
 209         xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 210         ftrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 211       }
 212     }
 213     else{
 214       TTree *xtree = (TTree*)fxsec->Get("Xsection");
 215       if(!xtree){
 216         Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
 217       }
 218       xtree->SetBranchAddress("xsection",&xsection);
 219       xtree->SetBranchAddress("ntrials",&ntrials);
 220       ftrials = ntrials;
 221       xtree->GetEntry(0);
 222     }
 223     fh1Xsec->Fill("<#sigma>",xsection);
 224     fh1Trials->Fill("#sum{ntrials}",ftrials);
 225   }
 226
 227   return kTRUE;
 228 }
 229
 230 void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
 231 {
 232
 233   //
 234   // Create the output container
 235   //
 236
 237
 238   // Connect the AOD
 239
 240   if(fUseAODInput){
 241     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 242     if(!fAOD){
 243       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 244       return;
 245     }
 246     // fethc the header
 247     fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 248     if(!fJetHeaderRec){
 249       Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
 250     }
 251   }
 252   else{
 253     //  assume that the AOD is in the general output...
 254     fAOD  = AODEvent();
 255     if(!fAOD){
 256       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 257       return;
 258     }
 259     fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
 260     if(!fJetHeaderRec){
 261       Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
 262     }
 263     else{
 264       if(fDebug>10)fJetHeaderRec->Dump();
 265     }
 266   }
 267
 268
 269
 270   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
 271
 272   OpenFile(1);
 273   if(!fHistList)fHistList = new TList();
 274
 275   Bool_t oldStatus = TH1::AddDirectoryStatus();
 276   TH1::AddDirectory(kFALSE);
 277
 278   //
 279   //  Histogram
 280
 281   const Int_t nBinPt = 100;
 282   Double_t binLimitsPt[nBinPt+1];
 283   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 284     if(iPt == 0){
 285       binLimitsPt[iPt] = 0.0;
 286     }
 287     else {// 1.0
 288       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
 289     }
 290   }
 291
 292   const Int_t nBinEta = 26;
 293   Double_t binLimitsEta[nBinEta+1] = {
 294     -1.6,-1.4,-1.2,-1.0,
 295     -0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
 296     0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
 297     1.0, 1.2, 1.4, 1.6
 298   };
 299
 300
 301   const Int_t nBinPhi = 30;
 302   Double_t binLimitsPhi[nBinPhi+1];
 303   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 304     if(iPhi==0){
 305       binLimitsPhi[iPhi] = 0;
 306     }
 307     else{
 308       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 309     }
 310   }
 311
 312   const Int_t nBinFrag = 25;
 313
 314
 315   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 316   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 317
 318   fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
 319   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 320
 321   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 322
 323   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 324
 325   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 326
 327   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 328
 329   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 330
 331
 332   for(int ij  = 0;ij<kMaxJets;++ij){
 333     fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
 334     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 335     fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
 336     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 337     fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
 338
 339
 340
 341     fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
 342                              nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 343
 344     fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};#phi_{gen}",
 345                              nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
 346
 347     fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};#eta_{gen}",
 348                              nBinEta,binLimitsEta,nBinEta,binLimitsEta);
 349
 350
 351     fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-#phi_{rec}",
 352                                     nBinPt,binLimitsPt,100,-1.0,1.0);
 353     fh2PtGenDeltaEta[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-#eta_{rec}",
 354                                     nBinPt,binLimitsPt,100,-1.0,1.0);
 355
 356
 357     fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,rec,j};#Delta R",
 358                                   nBinPt,binLimitsPt,60,0,6.0);
 359     fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,gen,j};#Delta R",
 360                                   nBinPt,binLimitsPt,60,0,6.0);
 361
 362
 363
 364     fh3PtRecGenHard[ij] = new TH3F(Form("fh3PtRecGenHard_j%d",ij), "Pt hard vs. pt gen vs. pt rec;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 365
 366
 367
 368     fh3PtRecGenHardNoW[ij] = new TH3F(Form("fh3PtRecGenHardNoW_j%d",ij), "Pt hard vs. pt gen vs. pt rec no weight;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 369
 370
 371     fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
 372                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 373
 374     fh2FragLn[ij] = new TH2F(Form("fh2FragLn_j%d",ij),"Jet Fragmentation Ln;#xi=ln(E_{jet}/E_{i});E_{jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
 375                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 376
 377     fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 378                                   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 379
 380
 381
 382     fh3RecEtaPhiPtNoGen[ij] = new TH3F(Form("fh3RecEtaPhiPtNoGen_j%d",ij),"No generated for found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
 383                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 384
 385
 386     fh3GenEtaPhiPtNoFound[ij] = new TH3F(Form("fh3GenEtaPhiPtNoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
 387                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 388
 389
 390
 391     fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
 392                                  nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 393
 394   }
 395
 396
 397   // tmp histos do not add to the header
 398   TH2F *hCorrPt = new TH2F("fh2PtRecPhiCorrPt","#Delta#phi correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
 399   fHistList->Add(hCorrPt);
 400   TH2F *hCorrRanPt = new TH2F("fh2PtRecPhiCorrPtRan","#Delta#phi Random correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
 401   fHistList->Add(hCorrRanPt);
 402
 403   TH2F *hCorr = new TH2F("fh2PtRecPhiCorr","#Delta#phi correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
 404   fHistList->Add(hCorr);
 405   TH2F *hCorrRan = new TH2F("fh2PtRecPhiCorrRan","#Delta#phi Random correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
 406   fHistList->Add(hCorrRan);
 407
 408
 409   /////////////////////////////////////////////////////////////////
 410   fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
 411
 412   fh2ERecZRec   = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
 413   fh2EGenZGen   = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
 414   fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);
 415
 416   fh3EGenERecN  = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100, 0., 250., 100, 0., 250., 51, 0., 50.);
 417
 418   // Response map
 419   //arrays for bin limits
 420   const Int_t nbin[4] = {100, 100, 100, 100};
 421   Double_t vLowEdge[4] = {0.,0.,0.,0.};
 422   Double_t vUpEdge[4] = {250., 250., 1., 1.};
 423
 424   for(int ic = 0;ic < kMaxCorrelation;ic++){
 425     fhnCorrelation[ic]  = new THnSparseF(Form("fhnCorrelation_%d",ic),  "Response Map", 4, nbin, vLowEdge, vUpEdge);
 426     if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
 427     else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
 428   }
 429   const Int_t saveLevel = 3; // large save level more histos
 430   if(saveLevel>0){
 431     fHistList->Add(fh1Xsec);
 432     fHistList->Add(fh1Trials);
 433     fHistList->Add(fh1PtHard);
 434     fHistList->Add(fh1PtHardNoW);
 435     fHistList->Add(fh1PtHardTrials);
 436     fHistList->Add(fh1NGenJets);
 437     fHistList->Add(fh1NRecJets);
 438     ////////////////////////
 439     fHistList->Add(fh1JetMultiplicity);
 440     fHistList->Add(fh2ERecZRec);
 441     fHistList->Add(fh2EGenZGen);
 442     fHistList->Add(fh2Efficiency);
 443     fHistList->Add(fh3EGenERecN);
 444
 445     for(int ic = 0;ic < kMaxCorrelation;++ic){
 446       fHistList->Add(fhnCorrelation[ic]);
 447     }
 448     ////////////////////////
 449     for(int ij  = 0;ij<kMaxJets;++ij){
 450       fHistList->Add(fh1E[ij]);
 451       fHistList->Add(fh1PtRecIn[ij]);
 452       fHistList->Add(fh1PtRecOut[ij]);
 453       fHistList->Add(fh1PtGenIn[ij]);
 454       fHistList->Add(fh1PtGenOut[ij]);
 455       fHistList->Add(fh2PtFGen[ij]);
 456       fHistList->Add(fh2PhiFGen[ij]);
 457       fHistList->Add(fh2EtaFGen[ij]);
 458       fHistList->Add(fh2PtGenDeltaEta[ij]);
 459       fHistList->Add(fh2PtGenDeltaPhi[ij]);
 460       fHistList->Add(fh2PtRecDeltaR[ij]);
 461       fHistList->Add(fh2PtGenDeltaR[ij]);
 462       fHistList->Add(fh3RecEtaPhiPt[ij]);
 463       fHistList->Add(fh3GenEtaPhiPt[ij]);
 464       if(saveLevel>2){
 465         fHistList->Add(fh3RecEtaPhiPtNoGen[ij]);
 466         fHistList->Add(fh3GenEtaPhiPtNoFound[ij]);
 467       }
 468     }
 469   }
 470
 471   // =========== Switch on Sumw2 for all histos ===========
 472   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 473     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 474     if (h1){
 475       // Printf("%s ",h1->GetName());
 476       h1->Sumw2();
 477       continue;
 478     }
 479     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 480     if(hn)hn->Sumw2();
 481   }
 482
 483   TH1::AddDirectory(oldStatus);
 484
 485 }
 486
 487 void AliAnalysisTaskJetSpectrum::Init()
 488 {
 489   //
 490   // Initialization
 491   //
 492
 493   Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
 494   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
 495
 496 }
 497
 498 void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
 499 {
 500   //
 501   // Execute analysis for current event
 502   //
 503
 504
 505
 506   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 507
 508
 509   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 510
 511   if(!aodH){
 512     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 513     return;
 514   }
 515
 516
 517   //  aodH->SelectEvent(kTRUE);
 518
 519   // ========= These pointers need to be valid in any case =======
 520
 521
 522   /*
 523   AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
 524   if(!jhRec){
 525     Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
 526     return;
 527   }
 528   */
 529   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 530   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 531   if(!aodRecJets){
 532     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 533     return;
 534   }
 535
 536   // ==== General variables needed
 537
 538
 539   // We use statice array, not to fragment the memory
 540   AliAODJet genJets[kMaxJets];
 541   Int_t nGenJets = 0;
 542   AliAODJet recJets[kMaxJets];
 543   Int_t nRecJets = 0;
 544   ///////////////////////////
 545   Int_t nTracks = 0;
 546   //////////////////////////
 547
 548   Double_t eventW = 1;
 549   Double_t ptHard = 0;
 550   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
 551
 552   if(fUseExternalWeightOnly){
 553     eventW = fExternalWeight;
 554   }
 555
 556
 557   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 558   if((fAnalysisType&kAnaMC)==kAnaMC){
 559     // this is the part we only use when we have MC information
 560     AliMCEvent* mcEvent = MCEvent();
 561     //    AliStack *pStack = 0;
 562     if(!mcEvent){
 563       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 564       return;
 565     }
 566     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 567     if(!pythiaGenHeader){
 568       return;
 569     }
 570
 571     nTrials = pythiaGenHeader->Trials();
 572     ptHard  = pythiaGenHeader->GetPtHard();
 573     int iProcessType = pythiaGenHeader->ProcessType();
 574     // 11 f+f -> f+f
 575     // 12 f+barf -> f+barf
 576     // 13 f+barf -> g+g
 577     // 28 f+g -> f+g
 578     // 53 g+g -> f+barf
 579     // 68 g+g -> g+g
 580     /*
 581     if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 582     //    if(iProcessType != 13 && iProcessType != 68){ // allow only glue
 583     if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark
 584     //    if(iProcessType != 28){ // allow only -> f+g
 585       PostData(1, fHistList);
 586       return;
 587     }
 588     */
 589     if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 590
 591     if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 592
 593     if(!fUseExternalWeightOnly){
 594         // case were we combine more than one p_T hard bin...
 595     }
 596
 597     // fetch the pythia generated jets only to be used here
 598     Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 599     AliAODJet pythiaGenJets[kMaxJets];
 600     Int_t iCount = 0;
 601     for(int ip = 0;ip < nPythiaGenJets;++ip){
 602       if(iCount>=kMaxJets)continue;
 603       Float_t p[4];
 604       pythiaGenHeader->TriggerJet(ip,p);
 605       pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 606
 607       if(fLimitGenJetEta){
 608         if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 609            pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 610       }
 611
 612
 613       if(fBranchGen.Length()==0){
 614         // if we have MC particles and we do not read from the aod branch
 615         // use the pythia jets
 616         genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 617       }
 618       iCount++;
 619     }
 620     if(fBranchGen.Length()==0)nGenJets = iCount;
 621
 622   }// (fAnalysisType&kMC)==kMC)
 623
 624   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 625   fh1PtHard->Fill(ptHard,eventW);
 626   fh1PtHardNoW->Fill(ptHard,1);
 627   fh1PtHardTrials->Fill(ptHard,nTrials);
 628
 629   // If we set a second branch for the input jets fetch this
 630   if(fBranchGen.Length()>0){
 631     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 632     if(aodGenJets){
 633       Int_t iCount = 0;
 634       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 635         if(iCount>=kMaxJets)continue;
 636         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 637         if(!tmp)continue;
 638         if(fLimitGenJetEta){
 639           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 640              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 641         }
 642         genJets[iCount] = *tmp;
 643         iCount++;
 644       }
 645       nGenJets = iCount;
 646     }
 647     else{
 648       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 649     }
 650   }
 651
 652   fh1NGenJets->Fill(nGenJets);
 653   // We do not want to exceed the maximum jet number
 654   nGenJets = TMath::Min(nGenJets,kMaxJets);
 655
 656   // Fetch the reconstructed jets...
 657
 658
 659   nRecJets = aodRecJets->GetEntries();
 660   fh1NRecJets->Fill(nRecJets);
 661   nRecJets = TMath::Min(nRecJets,kMaxJets);
 662   //////////////////////////////////////////
 663   nTracks  = fAOD->GetNumberOfTracks();
 664   ///////////////////////////////////////////
 665
 666   for(int ir = 0;ir < nRecJets;++ir){
 667     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 668     if(!tmp)continue;
 669     recJets[ir] = *tmp;
 670   }
 671
 672
 673   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 674   // Relate the jets
 675   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 676   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 677
 678   for(int i = 0;i<kMaxJets;++i){
 679     iGenIndex[i] = iRecIndex[i] = -1;
 680   }
 681
 682
 683   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 684                  iGenIndex,iRecIndex,fDebug);
 685   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 686
 687   if(fDebug){
 688     for(int i = 0;i<kMaxJets;++i){
 689       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 690       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 691     }
 692   }
 693
 694   // loop over reconstructed jets
 695   for(int ir = 0;ir < nRecJets;++ir){
 696     Double_t ptRec = recJets[ir].Pt();
 697     Double_t phiRec = recJets[ir].Phi();
 698     Double_t phiRecRan = TMath::Pi()*gRandom->Rndm(); // better take real jet axis from previous events (TPC acceptance in phi)
 699     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 700     Double_t etaRec = recJets[ir].Eta();
 701     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 702     fh1E[ir]->Fill(recJets[ir].E(),eventW);
 703     fh1PtRecIn[ir]->Fill(ptRec,eventW);
 704     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 705     for(int irr = ir+1;irr<nRecJets;irr++){
 706       fh2PtRecDeltaR[ir]->Fill(recJets[irr].Pt(),recJets[ir].DeltaR(&recJets[irr]));
 707     }
 708     // Fill Correlation
 709     Int_t ig = iGenIndex[ir];
 710     if(ig>=0 && ig<nGenJets){
 711       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 712       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 713       fh1PtRecOut[ir]->Fill(ptRec,eventW);
 714       Double_t ptGen  = genJets[ig].Pt();
 715       Double_t phiGen = genJets[ig].Phi();
 716       if(phiGen<0)phiGen+=TMath::Pi()*2.;
 717       Double_t etaGen = genJets[ig].Eta();
 718
 719       //
 720       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
 721       //
 722
 723       if(TMath::Abs(etaRec)<fRecEtaWindow){
 724
 725       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
 726       fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
 727       fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
 728       fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
 729       fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
 730       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
 731       fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
 732       /////////////////////////////////////////////////////
 733
 734       //      Double_t eRec = recJets[ir].E();
 735       //      Double_t eGen = genJets[ig].E();
 736       // CKB use p_T not Energy
 737       // TODO recname variabeles and histos
 738       Double_t eRec = recJets[ir].E();
 739       Double_t eGen = genJets[ig].E();
 740
 741       fh2Efficiency->Fill(eGen, eRec/eGen);
 742
 743       if (eGen>=0. && eGen<=250.){
 744         Double_t eLeading = -1;
 745         Double_t ptleading = -1;
 746         Int_t nPart=0;
 747         // loop over tracks
 748         for (Int_t it = 0; it< nTracks; it++){
 749           //      if (fAOD->GetTrack(it)->E() > eGen) continue; // CKB. Not allowed! cannot cut on gen properties in real events!
 750           // find leading particle
 751           //  if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
 752           // TODO implement esd filter flag to be the same as in the jet finder
 753           // allow also for MC particles...
 754           Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it));
 755           if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){
 756             eLeading  = fAOD->GetTrack(it)->E();
 757             ptleading = fAOD->GetTrack(it)->Pt();
 758           }
 759           //          if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7) // CKB cannot cut on gen properties
 760           if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7)
 761             nPart++;
 762
 763
 764           // correlate jet axis of leading jet with particles
 765           if(ir==0){
 766             Float_t phi =  fAOD->GetTrack(it)->Phi();
 767             Float_t dPhi = phi - phiRec;
 768             if(dPhi>TMath::Pi()/1.5)dPhi = dPhi - 2.*TMath::Pi();
 769             if(dPhi<(-0.5*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 770             Float_t dPhiRan = phi - phiRecRan;
 771             if(dPhiRan>TMath::Pi()/1.5)dPhiRan = dPhiRan - 2.*TMath::Pi();
 772             if(dPhiRan<(-0.5*TMath::Pi()))dPhiRan = dPhiRan + 2.*TMath::Pi();
 773             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorr"))->Fill(ptRec,dPhi);
 774             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrRan"))->Fill(ptRec,dPhiRan);
 775             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPt"))->Fill(ptRec,dPhi,fAOD->GetTrack(it)->Pt());
 776             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPtRan"))->Fill(ptRec,dPhiRan,fAOD->GetTrack(it)->Pt());
 777
 778           }
 779         }
 780         if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 781
 782         // fill Response Map (4D histogram) and Energy vs z distributions
 783         Double_t var[4] = {eGen, eRec, ptleading/eGen, ptleading/eRec};
 784         fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
 785         fh2EGenZGen->Fill(var[0],var[2]);
 786         fh1JetMultiplicity->Fill(nPart);
 787         fh3EGenERecN->Fill(eGen, eRec, nPart);
 788         for(int ic = 0;ic <kMaxCorrelation;ic++){
 789           if (nPart<=fgkJetNpartCut[ic]){ // is this corrected for CKB
 790             fhnCorrelation[ic]->Fill(var);
 791             break;
 792           }
 793         }
 794       }
 795
 796     }// if etarec in window
 797
 798     }
 799     ////////////////////////////////////////////////////
 800     else{
 801       fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 802     }
 803   }// loop over reconstructed jets
 804
 805
 806   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 807   for(int ig = 0;ig < nGenJets;++ig){
 808     Double_t ptGen = genJets[ig].Pt();
 809     // Fill Correlation
 810     Double_t phiGen = genJets[ig].Phi();
 811     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 812     Double_t etaGen = genJets[ig].Eta();
 813     fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 814     fh1PtGenIn[ig]->Fill(ptGen,eventW);
 815     for(int igg = ig+1;igg<nGenJets;igg++){
 816       fh2PtGenDeltaR[ig]->Fill(genJets[igg].Pt(),genJets[ig].DeltaR(&genJets[igg]));
 817     }
 818     Int_t ir = iRecIndex[ig];
 819     if(ir>=0&&ir<nRecJets){
 820       fh1PtGenOut[ig]->Fill(ptGen,eventW);
 821     }
 822     else{
 823       fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 824     }
 825   }// loop over reconstructed jets
 826
 827   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 828   PostData(1, fHistList);
 829 }
 830
 831 void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
 832 {
 833 // Terminate analysis
 834 //
 835     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
 836 }