Fix reading from AOD and changed warnings to AliWarning
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJFSystematics.cxx
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>
37 #include <TKey.h>
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"
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
65 ClassImp(AliAnalysisTaskJFSystematics)
66
67   const Int_t AliAnalysisTaskJFSystematics::fgkSysBins[AliAnalysisTaskJFSystematics::kSysTypes] = {0,AliAnalysisTaskJFSystematics::kMaxJets};
68 const char* AliAnalysisTaskJFSystematics::fgkSysName[AliAnalysisTaskJFSystematics::kSysTypes] = {"","j"};
69
70 AliAnalysisTaskJFSystematics::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
116 AliAnalysisTaskJFSystematics::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
168 Bool_t AliAnalysisTaskJFSystematics::Notify()
169 {
170 //
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;
177   Float_t   ftrials  = 0;
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")){
191         fileName.ReplaceAll("AliESDs.root", "");
192     }
193     else if(fileName.Contains("AliAOD.root")){
194         fileName.ReplaceAll("AliAOD.root", "");
195     }
196     else if(fileName.Contains("AliAODs.root")){
197         fileName.ReplaceAll("AliAODs.root", "");
198     }
199     else if(fileName.Contains("galice.root")){
200         // for running with galice and kinematics alone...                      
201         fileName.ReplaceAll("galice.root", "");
202     }
203     TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
204     if(!fxsec){
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       }
228     }
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);
238     }
239     fh1Xsec->Fill("<#sigma>",xsection);
240     fh1Trials->Fill("#sum{ntrials}",ftrials);
241   }
242   return kTRUE;
243 }
244
245 void AliAnalysisTaskJFSystematics::UserCreateOutputObjects()
246 {
247
248   //
249   // Create the output container
250   //
251
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);
327   fh2PtGenDeltaEta = new TH2F("fh2PtGenDeltaEta","delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-#eta_{rec}",
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
403 void 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
414 void AliAnalysisTaskJFSystematics::UserExec(Option_t */*option*/)
415 {
416   //
417   // Execute analysis for current event
418   //
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   
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
660 void AliAnalysisTaskJFSystematics::Terminate(Option_t */*option*/)
661 {
662 // Terminate analysis
663 //
664     if (fDebug > 1) printf("AnalysisTaskJFSystematics: Terminate() \n");
665     // Plot 
666
667
668 }