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