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