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