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