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   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 517
 518   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 519   if(!aodRecJets){
 520     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 521     return;
 522   }
 523
 524   // ==== General variables needed
 525
 526
 527   // We use statice array, not to fragment the memory
 528   AliAODJet genJets[kMaxJets];
 529   Int_t nGenJets = 0;
 530   AliAODJet recJets[kMaxJets];
 531   Int_t nRecJets = 0;
 532   ///////////////////////////
 533   Int_t nTracks = 0;
 534   //////////////////////////
 535
 536   Double_t eventW = 1;
 537   Double_t ptHard = 0;
 538   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
 539
 540   if(fUseExternalWeightOnly){
 541     eventW = fExternalWeight;
 542   }
 543
 544
 545   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 546   if((fAnalysisType&kAnaMC)==kAnaMC){
 547     // this is the part we only use when we have MC information
 548     AliMCEvent* mcEvent = MCEvent();
 549     //    AliStack *pStack = 0;
 550     if(!mcEvent){
 551       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 552       return;
 553     }
 554     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 555     if(!pythiaGenHeader){
 556       return;
 557     }
 558
 559     nTrials = pythiaGenHeader->Trials();
 560     ptHard  = pythiaGenHeader->GetPtHard();
 561     int iProcessType = pythiaGenHeader->ProcessType();
 562     // 11 f+f -> f+f
 563     // 12 f+barf -> f+barf
 564     // 13 f+barf -> g+g
 565     // 28 f+g -> f+g
 566     // 53 g+g -> f+barf
 567     // 68 g+g -> g+g
 568     /*
 569     if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 570     //    if(iProcessType != 13 && iProcessType != 68){ // allow only glue
 571     if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark
 572     //    if(iProcessType != 28){ // allow only -> f+g
 573       PostData(1, fHistList);
 574       return;
 575     }
 576     */
 577     if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 578
 579     if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 580
 581     if(!fUseExternalWeightOnly){
 582         // case were we combine more than one p_T hard bin...
 583     }
 584
 585     // fetch the pythia generated jets only to be used here
 586     Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 587     AliAODJet pythiaGenJets[kMaxJets];
 588     Int_t iCount = 0;
 589     for(int ip = 0;ip < nPythiaGenJets;++ip){
 590       if(iCount>=kMaxJets)continue;
 591       Float_t p[4];
 592       pythiaGenHeader->TriggerJet(ip,p);
 593       pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 594
 595       if(fLimitGenJetEta){
 596         if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 597            pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 598       }
 599
 600
 601       if(fBranchGen.Length()==0){
 602         // if we have MC particles and we do not read from the aod branch
 603         // use the pythia jets
 604         genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 605       }
 606       iCount++;
 607     }
 608     if(fBranchGen.Length()==0)nGenJets = iCount;
 609
 610   }// (fAnalysisType&kMC)==kMC)
 611
 612   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 613   fh1PtHard->Fill(ptHard,eventW);
 614   fh1PtHardNoW->Fill(ptHard,1);
 615   fh1PtHardTrials->Fill(ptHard,nTrials);
 616
 617   // If we set a second branch for the input jets fetch this
 618   if(fBranchGen.Length()>0){
 619     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 620     if(aodGenJets){
 621       Int_t iCount = 0;
 622       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 623         if(iCount>=kMaxJets)continue;
 624         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 625         if(!tmp)continue;
 626         if(fLimitGenJetEta){
 627           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 628              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 629         }
 630         genJets[iCount] = *tmp;
 631         iCount++;
 632       }
 633       nGenJets = iCount;
 634     }
 635     else{
 636       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 637     }
 638   }
 639
 640   fh1NGenJets->Fill(nGenJets);
 641   // We do not want to exceed the maximum jet number
 642   nGenJets = TMath::Min(nGenJets,kMaxJets);
 643
 644   // Fetch the reconstructed jets...
 645
 646
 647   nRecJets = aodRecJets->GetEntries();
 648   fh1NRecJets->Fill(nRecJets);
 649   nRecJets = TMath::Min(nRecJets,kMaxJets);
 650   //////////////////////////////////////////
 651   nTracks  = fAOD->GetNumberOfTracks();
 652   ///////////////////////////////////////////
 653
 654   for(int ir = 0;ir < nRecJets;++ir){
 655     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 656     if(!tmp)continue;
 657     recJets[ir] = *tmp;
 658   }
 659
 660
 661   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 662   // Relate the jets
 663   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 664   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 665
 666   for(int i = 0;i<kMaxJets;++i){
 667     iGenIndex[i] = iRecIndex[i] = -1;
 668   }
 669
 670
 671   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 672                  iGenIndex,iRecIndex,fDebug);
 673   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 674
 675   if(fDebug){
 676     for(int i = 0;i<kMaxJets;++i){
 677       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 678       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 679     }
 680   }
 681
 682   // loop over reconstructed jets
 683   for(int ir = 0;ir < nRecJets;++ir){
 684     Double_t ptRec = recJets[ir].Pt();
 685     Double_t phiRec = recJets[ir].Phi();
 686     Double_t phiRecRan = TMath::Pi()*gRandom->Rndm(); // better take real jet axis from previous events (TPC acceptance in phi)
 687     if(phiRec<0)phiRec+=TMath::Pi()*2.;
 688     Double_t etaRec = recJets[ir].Eta();
 689     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 690     fh1E[ir]->Fill(recJets[ir].E(),eventW);
 691     fh1PtRecIn[ir]->Fill(ptRec,eventW);
 692     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 693     for(int irr = ir+1;irr<nRecJets;irr++){
 694       fh2PtRecDeltaR[ir]->Fill(recJets[irr].Pt(),recJets[ir].DeltaR(&recJets[irr]));
 695     }
 696     // Fill Correlation
 697     Int_t ig = iGenIndex[ir];
 698     if(ig>=0 && ig<nGenJets){
 699       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 700       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 701       fh1PtRecOut[ir]->Fill(ptRec,eventW);
 702       Double_t ptGen  = genJets[ig].Pt();
 703       Double_t phiGen = genJets[ig].Phi();
 704       if(phiGen<0)phiGen+=TMath::Pi()*2.;
 705       Double_t etaGen = genJets[ig].Eta();
 706
 707       //
 708       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
 709       //
 710
 711       if(TMath::Abs(etaRec)<fRecEtaWindow){
 712
 713       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
 714       fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
 715       fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
 716       fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
 717       fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
 718       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
 719       fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
 720       /////////////////////////////////////////////////////
 721
 722       //      Double_t eRec = recJets[ir].E();
 723       //      Double_t eGen = genJets[ig].E();
 724       // CKB use p_T not Energy
 725       // TODO recname variabeles and histos
 726       Double_t eRec = recJets[ir].E();
 727       Double_t eGen = genJets[ig].E();
 728
 729       fh2Efficiency->Fill(eGen, eRec/eGen);
 730
 731       if (eGen>=0. && eGen<=250.){
 732         Double_t eLeading = -1;
 733         Double_t ptleading = -1;
 734         Int_t nPart=0;
 735         // loop over tracks
 736         for (Int_t it = 0; it< nTracks; it++){
 737           //      if (fAOD->GetTrack(it)->E() > eGen) continue; // CKB. Not allowed! cannot cut on gen properties in real events!
 738           // find leading particle
 739           //  if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
 740           // TODO implement esd filter flag to be the same as in the jet finder
 741           // allow also for MC particles...
 742           Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it));
 743           if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){
 744             eLeading  = fAOD->GetTrack(it)->E();
 745             ptleading = fAOD->GetTrack(it)->Pt();
 746           }
 747           //          if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7) // CKB cannot cut on gen properties
 748           if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7)
 749             nPart++;
 750
 751
 752           // correlate jet axis of leading jet with particles
 753           if(ir==0){
 754             Float_t phi =  fAOD->GetTrack(it)->Phi();
 755             Float_t dPhi = phi - phiRec;
 756             if(dPhi>TMath::Pi()/1.5)dPhi = dPhi - 2.*TMath::Pi();
 757             if(dPhi<(-0.5*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 758             Float_t dPhiRan = phi - phiRecRan;
 759             if(dPhiRan>TMath::Pi()/1.5)dPhiRan = dPhiRan - 2.*TMath::Pi();
 760             if(dPhiRan<(-0.5*TMath::Pi()))dPhiRan = dPhiRan + 2.*TMath::Pi();
 761             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorr"))->Fill(ptRec,dPhi);
 762             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrRan"))->Fill(ptRec,dPhiRan);
 763             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPt"))->Fill(ptRec,dPhi,fAOD->GetTrack(it)->Pt());
 764             ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPtRan"))->Fill(ptRec,dPhiRan,fAOD->GetTrack(it)->Pt());
 765
 766           }
 767         }
 768         if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 769
 770         // fill Response Map (4D histogram) and Energy vs z distributions
 771         Double_t var[4] = {eGen, eRec, ptleading/eGen, ptleading/eRec};
 772         fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
 773         fh2EGenZGen->Fill(var[0],var[2]);
 774         fh1JetMultiplicity->Fill(nPart);
 775         fh3EGenERecN->Fill(eGen, eRec, nPart);
 776         for(int ic = 0;ic <kMaxCorrelation;ic++){
 777           if (nPart<=fgkJetNpartCut[ic]){ // is this corrected for CKB
 778             fhnCorrelation[ic]->Fill(var);
 779             break;
 780           }
 781         }
 782       }
 783
 784     }// if etarec in window
 785
 786     }
 787     ////////////////////////////////////////////////////
 788     else{
 789       fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
 790     }
 791   }// loop over reconstructed jets
 792
 793
 794   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 795   for(int ig = 0;ig < nGenJets;++ig){
 796     Double_t ptGen = genJets[ig].Pt();
 797     // Fill Correlation
 798     Double_t phiGen = genJets[ig].Phi();
 799     if(phiGen<0)phiGen+=TMath::Pi()*2.;
 800     Double_t etaGen = genJets[ig].Eta();
 801     fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 802     fh1PtGenIn[ig]->Fill(ptGen,eventW);
 803     for(int igg = ig+1;igg<nGenJets;igg++){
 804       fh2PtGenDeltaR[ig]->Fill(genJets[igg].Pt(),genJets[ig].DeltaR(&genJets[igg]));
 805     }
 806     Int_t ir = iRecIndex[ig];
 807     if(ir>=0&&ir<nRecJets){
 808       fh1PtGenOut[ig]->Fill(ptGen,eventW);
 809     }
 810     else{
 811       fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
 812     }
 813   }// loop over reconstructed jets
 814
 815   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 816   PostData(1, fHistList);
 817 }
 818
 819 void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
 820 {
 821 // Terminate analysis
 822 //
 823     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
 824 }