New book-keeping class
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetCorrections.cxx
1 // **************************************
2 // Task used for estimating a charged to neutral correction
3 // sona.pochybova@cern.ch
4 // *******************************************
5
6
7 /**************************************************************************
8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9  *                                                                        *
10  * Author: The ALICE Off-line Project.                                    *
11  * Contributors are mentioned in the code where appropriate.              *
12  *                                                                        *
13  * Permission to use, copy, modify and distribute this software and its   *
14  * documentation strictly for non-commercial purposes is hereby granted   *
15  * without fee, provided that the above copyright notice appears in all   *
16  * copies and that both the copyright notice and this permission notice   *
17  * appear in the supporting documentation. The authors make no claims     *
18  * about the suitability of this software for any purpose. It is          *
19  * provided "as is" without express or implied warranty.                  *
20  **************************************************************************/
21
22
23 #include <TFile.h>
24 #include <TH1F.h>
25 #include <TH2F.h>
26 #include <TList.h>
27 #include <TTree.h>
28 #include <TLorentzVector.h>
29 #include <TClonesArray.h>
30 #include <TRefArray.h>
31 #include <TVector3.h>
32 #include <TProfile.h>
33
34 #include "AliAnalysisTaskJetCorrections.h"
35 #include "AliAnalysisManager.h"
36 #include "AliAODEvent.h"
37 #include "AliAODVertex.h"
38 #include "AliAODHandler.h"
39 #include "AliAODTrack.h"
40 #include "AliAODJet.h"
41 //#include "AliMCEvent.h"
42
43 #include "AliAnalysisHelperJetTasks.h"
44
45 //
46 //
47 // corrections to jet energy by sona 
48 // 
49 //
50
51 ClassImp(AliAnalysisTaskJetCorrections)
52
53   AliAnalysisTaskJetCorrections::AliAnalysisTaskJetCorrections() : AliAnalysisTaskSE(),                                                            
54                                                                    fAOD(0x0),
55                                                                   
56                                                                    fBranchRec(""),
57                                                                    fBranchGen(""),
58                                                                    
59                                                                    fUseAODInput(kFALSE),
60
61                                                                    fR(0x0),
62                                                                    fList(0x0),
63                                                                    
64                                                                    fGlobVar(1),
65                                                                    fXsection(1),
66                                                                    
67                                                                    fhEGen(0x0),
68                                                                    fhERec(0x0),
69                                                                    
70                                                                    fhEGenRest(0x0),
71                                                                    fhERecRest(0x0),
72
73                                                                    fhEsumGenRest(0x0),
74                                                                    fhEsumRecRest(0x0),
75                                                                    
76                                                                    fhE2vsE1Gen(0x0),
77                                                                    fhE2vsE1Rec(0x0),
78                                                                    fhE2E1vsEsumGen(0x0),
79                                                                    fhE2E1vsEsumRec(0x0),
80                                                                    fhE2E1vsE1Gen(0x0),
81                                                                    fhE2E1vsE1Rec(0x0),
82                                                                    fhE2E1vsdPhiGen(0x0),
83                                                                    fhE2E1vsdPhiRec(0x0),
84
85                                                                    fhTrackBalance2(0x0),
86                                                                    fhTrackBalance3(0x0),
87
88                                                                    fhEt1Et22(0x0),
89                                                                    fhEt1Et23(0x0)
90
91 {
92   //
93   // ctor
94   //
95   for (Int_t i = 0; i < 3; i++)
96     {
97       fhECorrJet10[i] = 0;    
98       fhECorrJet05[i] = 0;    
99       fhECorrJet01[i] = 0;    
100       fhECorrJet001[i] = 0;
101       fhdEvsErec10[i] = 0;
102       fhdEvsErec05[i] = 0;
103       fhdEvsErec01[i] = 0;
104       fhdEvsErec001[i] = 0;
105       fhdPhidEta10[i] = 0;
106       fhdPhidEta05[i] = 0;
107       fhdPhidEta01[i] = 0;
108       fhdPhidEta001[i] = 0;
109       fhdPhidEtaPt10[i] = 0;
110       fhdPhidEtaPt05[i] = 0;
111       fhdPhidEtaPt01[i] = 0;
112       fhdPhidEtaPt001[i] = 0;
113     }
114 }
115
116 AliAnalysisTaskJetCorrections::AliAnalysisTaskJetCorrections(const char * name):
117   AliAnalysisTaskSE(name),
118   
119   fAOD(0x0),
120   
121   fBranchRec(""),
122   fBranchGen(""),
123   
124   fUseAODInput(kFALSE),
125
126   fR(0x0),
127   fList(0x0),
128
129   fGlobVar(1),
130   fXsection(1),
131
132   fhEGen(0x0),
133   fhERec(0x0),
134
135   fhEGenRest(0x0),
136   fhERecRest(0x0),
137   
138   fhEsumGenRest(0x0),
139   fhEsumRecRest(0x0),
140   
141   fhE2vsE1Gen(0x0),
142   fhE2vsE1Rec(0x0),
143   fhE2E1vsEsumGen(0x0),
144   fhE2E1vsEsumRec(0x0),
145   fhE2E1vsE1Gen(0x0),
146   fhE2E1vsE1Rec(0x0),
147   fhE2E1vsdPhiGen(0x0),
148   fhE2E1vsdPhiRec(0x0),
149
150   fhTrackBalance2(0x0),
151   fhTrackBalance3(0x0),
152
153   fhEt1Et22(0x0),
154   fhEt1Et23(0x0)
155 {
156   //
157   // ctor
158   //
159   for (Int_t i = 0; i < 3; i++)
160     {
161       fhECorrJet10[i] = 0;    
162       fhECorrJet05[i] = 0;    
163       fhECorrJet01[i] = 0;    
164       fhECorrJet001[i] = 0;
165       fhdEvsErec10[i] = 0;
166       fhdEvsErec05[i] = 0;
167       fhdEvsErec01[i] = 0;
168       fhdEvsErec001[i] = 0;
169       fhdPhidEta10[i] = 0;
170       fhdPhidEta05[i] = 0;
171       fhdPhidEta01[i] = 0;
172       fhdPhidEta001[i] = 0;
173       fhdPhidEtaPt10[i] = 0;
174       fhdPhidEtaPt05[i] = 0;
175       fhdPhidEtaPt01[i] = 0;
176       fhdPhidEtaPt001[i] = 0;
177     }
178   DefineOutput(1, TList::Class());
179 }
180
181
182
183 Bool_t AliAnalysisTaskJetCorrections::Notify()
184 {
185   //
186   // Implemented Notify() to read the cross sections
187   // and number of trials from pyxsec.root
188   // 
189
190
191   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
192   Float_t xsection = 0;
193   Float_t ftrials  = 1;
194
195   Float_t fAvgTrials = 1;
196   if(tree){
197     TFile *curfile = tree->GetCurrentFile();
198     if (!curfile) {
199       Error("Notify","No current file");
200       return kFALSE;
201     }
202     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
203     // construct a poor man average trials 
204     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
205     if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; // CKB take this into account for normalisation
206   }  
207
208   if(xsection>0)fXsection  = xsection;
209
210   return kTRUE;
211
212 }
213
214
215 //___________________________________________________________________________________________________________________________________
216 void AliAnalysisTaskJetCorrections::UserCreateOutputObjects()
217 {
218   //
219   // Create the output container
220   //
221   //  Printf("Analysing event  %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
222
223   if(fUseAODInput){
224     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
225     if(!fAOD){
226       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
227       return;
228     }    
229   }
230   else{
231     //  assume that the AOD is in the general output...
232     fAOD  = AODEvent();
233     if(!fAOD){
234       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
235       return;
236     }    
237   }
238
239   printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
240
241   fList = new TList();
242
243   fhEGen = new TH1F("EGen", "", 100, 0, 200);
244   fhEGen->Sumw2();
245   fList->Add(fhEGen);
246
247   fhERec = new TH1F("ERec", "", 100, 0, 200);
248   fhERec->Sumw2();
249   fList->Add(fhERec);
250
251   fhEGenRest = new TH1F("EGenRest", "", 100, 0, 200);
252   fhEGenRest->Sumw2();
253   fList->Add(fhEGenRest);
254
255   fhERecRest = new TH1F("ERecRest", "", 100, 0, 200);
256   fhERecRest->Sumw2();
257   fList->Add(fhERecRest);
258
259   fhEsumGenRest = new TH1F("EsumGenRest", "", 100, 0, 200);
260   fhEsumGenRest->Sumw2();
261   fList->Add(fhEsumGenRest);
262
263   fhEsumRecRest = new TH1F("EsumRecRest", "", 100, 0, 200);
264   fhEsumRecRest->Sumw2();
265   fList->Add(fhEsumRecRest);
266
267   fhE2vsE1Gen = new TH2F("E2vsE1Gen", "", 100, 0, 200, 100, 0, 200);
268   fhE2vsE1Gen->Sumw2();
269   fList->Add(fhE2vsE1Gen);
270
271   fhE2vsE1Rec = new TH2F("E2vsE1Rec", "", 100, 0, 200, 100, 0, 200);
272   fhE2vsE1Rec->Sumw2();
273   fList->Add(fhE2vsE1Rec);
274
275   fhE2E1vsEsumGen = new TH2F("E2E1vsEsumGen", "", 100, 0, 200, 25, 0, 1);
276   fhE2E1vsEsumGen->Sumw2();
277   fList->Add(fhE2E1vsEsumGen);
278
279   fhE2E1vsEsumRec = new TH2F("E2E1vsEsumRec", "", 100, 0, 200, 25, 0, 1);
280   fhE2E1vsEsumRec->Sumw2();
281   fList->Add(fhE2E1vsEsumRec);
282
283   fhE2E1vsE1Gen = new TH2F("E2E1vsE1Gen", "", 100, 0, 200, 25, 0, 1);
284   fhE2E1vsE1Gen->Sumw2();
285   fList->Add(fhE2E1vsE1Gen);
286
287   fhE2E1vsE1Rec = new TH2F("E2E1vsE1Rec", "", 100, 0, 200, 25, 0, 1);
288   fhE2E1vsE1Rec->Sumw2();
289   fList->Add(fhE2E1vsE1Rec);
290
291   fhE2E1vsdPhiGen =  new TH2F("E2E1vsdPhiGen", "", 64, -3.20, 3.20, 25, 0, 1);
292   fList->Add(fhE2E1vsdPhiGen);
293
294   fhE2E1vsdPhiRec =  new TH2F("E2E1vsdPhiRec", "", 64, -3.20, 3.20, 25, 0, 1);
295   fList->Add(fhE2E1vsdPhiRec);
296   
297   fhTrackBalance2 = new TH2F("TrackBalance2", "", 60, 0, 30, 60, 0, 30);
298   fhTrackBalance2->Sumw2();
299   fList->Add(fhTrackBalance2);
300   
301   fhTrackBalance3 = new TH2F("TrackBalance3", "", 60, 0, 30, 60, 0, 30);
302   fhTrackBalance3->Sumw2();
303   fList->Add(fhTrackBalance3);
304
305   fhEt1Et22 = new TH2F("Et1Et22", "", 100, 0, 50, 100, 0, 50);
306   fhEt1Et22->Sumw2();
307   fList->Add(fhEt1Et22);
308
309   fhEt1Et23 = new TH2F("Et1Et23", "", 100, 0, 50, 100, 0, 50);
310   fhEt1Et23->Sumw2();
311   fList->Add(fhEt1Et23);
312
313   for(Int_t i = 0; i < 3; i++)
314     {
315       fhECorrJet10[i] = new TProfile(Form("ECorrJet10%d", i+1), "", 100, 0, 200, 0, 10);
316       fhECorrJet10[i]->SetXTitle("E_{rec} [GeV]");
317       fhECorrJet10[i]->SetYTitle("C=E_{gen}/E_{rec}");
318       fhECorrJet10[i]->Sumw2();
319
320       fhECorrJet05[i] = new TProfile(Form("ECorrJet05%d", i+1), "", 100, 0, 200, 0, 10);
321       fhECorrJet05[i]->SetXTitle("E_{rec} [GeV]");
322       fhECorrJet05[i]->SetYTitle("C=E_{gen}/E_{rec}");
323       fhECorrJet05[i]->Sumw2();
324
325       fhECorrJet01[i] = new TProfile(Form("ECorrJet01%d", i+1), "", 100, 0, 200, 0, 10);
326       fhECorrJet01[i]->SetXTitle("E_{rec} [GeV]");
327       fhECorrJet01[i]->SetYTitle("C=E_{gen}/E_{rec}");
328       fhECorrJet01[i]->Sumw2();
329
330       fhECorrJet001[i] = new TProfile(Form("ECorrJet001%d", i+1), "", 100, 0, 200, 0, 10);
331       fhECorrJet001[i]->SetXTitle("E_{rec} [GeV]");
332       fhECorrJet001[i]->SetYTitle("C=E_{gen}/E_{rec}");
333       fhECorrJet001[i]->Sumw2();
334
335       fhdEvsErec10[i] = new TProfile(Form("dEvsErec10_%d", i+1),"", 100, 0, 200, -1, 10);
336       fhdEvsErec10[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
337       fhdEvsErec10[i]->SetXTitle("E_{rec} [GeV]");
338       fhdEvsErec10[i]->Sumw2();
339
340       fhdEvsErec05[i] = new TProfile(Form("dEvsErec05_%d", i+1),"", 100, 0, 200, -1, 10);
341       fhdEvsErec05[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
342       fhdEvsErec05[i]->SetXTitle("E_{rec} [GeV]");
343       fhdEvsErec05[i]->Sumw2();
344
345       fhdEvsErec01[i] = new TProfile(Form("dEvsErec01_%d", i+1),"", 100, 0, 200, -1, 10);
346       fhdEvsErec01[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
347       fhdEvsErec01[i]->SetXTitle("E_{rec} [GeV]");
348       fhdEvsErec01[i]->Sumw2();
349
350       fhdEvsErec001[i] = new TProfile(Form("dEvsErec001_%d", i+1),"", 100, 0, 200, -1, 10);
351       fhdEvsErec001[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
352       fhdEvsErec001[i]->SetXTitle("E_{rec} [GeV]");
353       fhdEvsErec001[i]->Sumw2();
354
355       fhdPhidEta10[i] = new TH2F(Form("dPhidEta10_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
356       fhdPhidEta10[i]->SetXTitle("#phi [rad]");
357       fhdPhidEta10[i]->SetYTitle("#eta");
358       fhdPhidEta10[i]->Sumw2();
359
360       fhdPhidEta05[i] = new TH2F(Form("dPhidEta05_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
361       fhdPhidEta05[i]->SetXTitle("#phi [rad]");
362       fhdPhidEta05[i]->SetYTitle("#eta");
363       fhdPhidEta05[i]->Sumw2();
364
365       fhdPhidEta01[i] = new TH2F(Form("dPhidEta01_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
366       fhdPhidEta01[i]->SetXTitle("#phi [rad]");
367       fhdPhidEta01[i]->SetYTitle("#eta");
368       fhdPhidEta01[i]->Sumw2();
369
370       fhdPhidEta001[i] = new TH2F(Form("dPhidEta001_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
371       fhdPhidEta001[i]->SetXTitle("#phi [rad]");
372       fhdPhidEta001[i]->SetYTitle("#eta");
373       fhdPhidEta001[i]->Sumw2();
374
375       fhdPhidEtaPt10[i] = new TH2F(Form("dPhidEtaPt10_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
376       fhdPhidEtaPt10[i]->SetXTitle("#phi [rad]");
377       fhdPhidEtaPt10[i]->SetYTitle("#eta");
378       fhdPhidEtaPt10[i]->Sumw2();
379
380       fhdPhidEtaPt05[i] = new TH2F(Form("dPhidEtaPt05_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
381       fhdPhidEtaPt05[i]->SetXTitle("#phi [rad]");
382       fhdPhidEtaPt05[i]->SetYTitle("#eta");
383       fhdPhidEtaPt05[i]->Sumw2();
384
385       fhdPhidEtaPt01[i] = new TH2F(Form("dPhidEtaPt01_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
386       fhdPhidEtaPt01[i]->SetXTitle("#phi [rad]");
387       fhdPhidEtaPt01[i]->SetYTitle("#eta");
388       fhdPhidEtaPt01[i]->Sumw2();
389
390       fhdPhidEtaPt001[i] = new TH2F(Form("dPhidEtaPt001_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
391       fhdPhidEtaPt001[i]->SetXTitle("#phi [rad]");
392       fhdPhidEtaPt001[i]->SetYTitle("#eta");
393       fhdPhidEtaPt001[i]->Sumw2();
394
395       fList->Add(fhECorrJet10[i]);
396       fList->Add(fhECorrJet05[i]);
397       fList->Add(fhECorrJet01[i]);
398       fList->Add(fhECorrJet001[i]);
399       fList->Add(fhdEvsErec10[i]);
400       fList->Add(fhdEvsErec05[i]);
401       fList->Add(fhdEvsErec01[i]);
402       fList->Add(fhdEvsErec001[i]);
403       fList->Add(fhdPhidEta10[i]);
404       fList->Add(fhdPhidEta05[i]);
405       fList->Add(fhdPhidEta01[i]);
406       fList->Add(fhdPhidEta001[i]);
407       fList->Add(fhdPhidEtaPt10[i]);
408       fList->Add(fhdPhidEtaPt05[i]);
409       fList->Add(fhdPhidEtaPt01[i]);
410       fList->Add(fhdPhidEtaPt001[i]);
411     }
412     
413   Printf("UserCreateOutputObjects finished\n");
414 }
415
416 //__________________________________________________________________________________________________________________________________________
417 void AliAnalysisTaskJetCorrections::Init()
418 {
419   printf("AliAnalysisJetCut::Init() \n");
420 }
421
422 //____________________________________________________________________________________________________________________________________________
423 void AliAnalysisTaskJetCorrections::UserExec(Option_t * )
424 {
425 //  if (fDebug > 1) printf("Analysing event # %5d\n", (Int_t) fEntry);
426
427  
428   //create an AOD handler
429   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
430   
431   if(!aodH)
432     {
433       Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
434       return;
435     }
436   
437 //   AliMCEvent* mcEvent =MCEvent();
438 //   if(!mcEvent){
439 //     Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
440 //     return;
441 //   }
442   
443   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
444
445  //primary vertex
446   AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
447   if(!pvtx) return;
448
449  AliAODJet genJets[kMaxJets];
450  Int_t nGenJets = 0;
451  
452  AliAODJet recJets[kMaxJets];
453  Int_t nRecJets = 0;
454
455   //array of reconstructed jets from the AOD input
456  TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
457  if(!aodRecJets){
458    Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
459    return;
460  }
461  
462  // reconstructed jets
463  nRecJets = aodRecJets->GetEntries(); 
464  nRecJets = TMath::Min(nRecJets, kMaxJets);
465  
466  for(int ir = 0;ir < nRecJets;++ir)
467    {
468      AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
469      if(!tmp)continue;
470      recJets[ir] = *tmp;
471     }
472  
473  // If we set a second branch for the input jets fetch this
474  TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
475  
476  if(!aodGenJets)
477    {
478      printf("NO MC jets branch with name %s  Found \n",fBranchGen.Data());
479      return;
480    }
481  
482  //   //Generated jets
483  nGenJets = aodGenJets->GetEntries();
484  nGenJets = TMath::Min(nGenJets, kMaxJets);
485  
486  for(Int_t ig =0 ; ig < nGenJets; ++ig)
487    {
488      AliAODJet * tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
489      if(!tmp)continue;
490      genJets[ig] = * tmp;
491    }
492
493  Double_t eRec[kMaxJets];
494  Double_t eGen[kMaxJets];
495  
496  Double_t eRecRest[kMaxJets];
497  Double_t eGenRest[kMaxJets];
498
499 // AliAODJet jetRec[kMaxJets];
500  AliAODJet jetGen[kMaxJets];
501  
502  Int_t idxRec[kMaxJets];
503  Int_t idxGen[kMaxJets];
504
505 // Double_t EsumRec = 0;
506  // Double_t EsumGen =0;
507  
508  TLorentzVector vRec[kMaxJets];
509  TLorentzVector vGen[kMaxJets];
510
511  TLorentzVector vsumRec;
512  TLorentzVector vsumGen;
513  
514  TVector3 pRec[kMaxJets];
515  TVector3 pGen[kMaxJets];
516  
517  // HISTOS FOR MC
518  Int_t nGenSel = 0;
519  Int_t counter = 0;
520  Int_t tag = 0;
521  
522  AliAODJet selJets[kMaxJets];
523  
524  // loop for applying the separation cut
525  for(Int_t i = 0; i < nGenJets; i++)
526    {
527      if(nGenJets == 1)
528        {
529          selJets[nGenSel] = genJets[i];
530          nGenSel++;
531        }
532      else
533        {
534          counter = 0;
535          tag = 0;
536          for(Int_t j = 0; j < nGenJets; j++)
537            {
538              if(i!=j)
539                {
540                  Double_t dRij = genJets[i].DeltaR(&genJets[j]);
541                  counter++;
542                  if(dRij > 2*fR) tag++;
543                }
544            }
545          if(counter!=0)
546            {
547              if(tag/counter == 1)
548                {
549                  selJets[nGenSel] = genJets[i];
550                  nGenSel++;
551                }
552            }
553        }
554    }
555  
556  for (Int_t gj = 0; gj < nGenSel; gj++)
557    {
558      eGen[gj] = selJets[gj].E();
559      fhEGen->Fill(eGen[gj], fXsection);
560    }
561  
562  TMath::Sort(nGenSel, eGen, idxGen);
563  for (Int_t ig = 0; ig < nGenSel; ig++)
564    jetGen[ig] = selJets[idxGen[ig]];
565  
566  //rest frame MC jets
567  for (Int_t i = 0; i < nGenSel; ++i)
568    {
569      vGen[i].SetPxPyPzE(jetGen[i].Px(), jetGen[i].Py(), jetGen[i].Pz(), jetGen[i].E());
570      pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
571      vsumGen += vGen[i];
572    }
573  
574  if(nGenSel > 1 && pGen[0].DeltaPhi(pGen[1]) > 2.8)
575    {
576      fhE2vsE1Gen->Fill(jetGen[0].E(), jetGen[1].E(), fXsection);
577      fhE2E1vsEsumGen->Fill(jetGen[0].E()+jetGen[1].E(), TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
578      fhE2E1vsE1Gen->Fill(jetGen[0].E(),  TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
579      Double_t deltaPhi = (jetGen[0].Phi()-jetGen[1].Phi());
580      if(deltaPhi > TMath::Pi()) deltaPhi = deltaPhi - 2.*TMath::Pi();
581      if(deltaPhi < (-1.*TMath::Pi())) deltaPhi = deltaPhi + 2.*TMath::Pi();
582      fhE2E1vsdPhiGen->Fill(deltaPhi,  TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
583    }
584      
585  Double_t fPxGen = vsumGen.Px();
586  Double_t fPyGen = vsumGen.Py();
587  Double_t fPzGen = vsumGen.Pz();
588  Double_t fEGen = vsumGen.E();
589  
590  Double_t eSumGenRest = 0; 
591  for (Int_t j = 0; j < nGenSel; j++)
592    {
593      vGen[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
594      eGenRest[j] = vGen[j].E();
595      if(nGenSel > 1)
596        fhEGenRest->Fill(eGenRest[j], fXsection);
597      eSumGenRest += eGenRest[j];
598    }
599  
600  if(nGenSel > 1)    
601    fhEsumGenRest->Fill(eSumGenRest, fXsection);
602  
603  //END VARIABLES FOR MC JETS ---------------
604  //   }
605  
606  //AOD JET VARIABLES
607  Int_t nRecSel = 0;
608  Int_t counter1 = 0;
609  Int_t tag1 = 0;
610  
611  AliAODJet recSelJets[kMaxJets];
612  
613  for(Int_t i = 0; i < nRecJets; i++)
614    {
615      if(nRecJets == 1)
616        {
617          recSelJets[nRecSel] = recJets[i];
618          nRecSel++;
619        }
620      else
621        {
622          counter1 = 0;
623          tag1 = 0;
624          for(Int_t j = 0; j < nRecJets; j++)
625            {
626              if(i!=j)
627                {
628                  Double_t dRij = recJets[i].DeltaR(&recJets[j]);
629                  counter1++;
630                  if(dRij > 2*fR) tag1++;
631                }
632            }
633          if(counter1!=0)
634            {
635              if(tag1/counter1 == 1)
636                {
637                  recSelJets[nRecSel] = recJets[i];
638                  nRecSel++;
639                }
640            }
641        }
642    } 
643   
644  if(nRecSel == 0) return;
645  Printf("******NUMBER OF JETS AFTER DELTA R CUT : %d **********\n", nRecSel);
646  //sort rec/gen jets by energy in C.M.S
647  AliAODJet jetRecTmp[kMaxJets];
648  Int_t nAccJets = 0;
649  Double_t jetTrackPt[kTracks];
650  TLorentzVector jetTrackTmp[kTracks];
651  Int_t nTracks = 0;
652  for (Int_t rj = 0; rj < nRecSel; rj++)
653    {
654      TRefArray * jetTracksAOD = dynamic_cast<TRefArray*>(recSelJets[rj].GetRefTracks());
655      if(!jetTracksAOD) continue;
656      if(jetTracksAOD->GetEntries() < 3) continue;
657      Int_t nJetTracks = 0;
658      for(Int_t j = 0; j < jetTracksAOD->GetEntries(); j++)
659        {
660          AliAODTrack * track = dynamic_cast<AliAODTrack*>(jetTracksAOD->At(j));
661          if(!track) continue;
662          Double_t cv[21];
663          track->GetCovarianceXYZPxPyPz(cv);
664          if(cv[14] > 1000.) continue;
665          jetTrackPt[nTracks] = track->Pt();
666          jetTrackTmp[nTracks].SetPxPyPzE(track->Px(),track->Py(),track->Pz(),track->E());
667          nTracks++;
668          nJetTracks++;
669        }
670      if(nJetTracks < 4) continue;
671      jetRecTmp[nAccJets] = recSelJets[rj];
672      eRec[nAccJets] = recSelJets[rj].E();
673      fhERec->Fill(eRec[nAccJets], fXsection);
674      nAccJets++;
675    }
676
677  if(nAccJets == 0) return;
678  if(nTracks == 0) return;
679
680  Printf(" ************ Number of accepted jets : %d ************ \n", nAccJets);
681
682  AliAODJet jetRecAcc[kMaxJets];
683  TMath::Sort(nAccJets, eRec, idxRec);
684  for (Int_t rj = 0; rj < nAccJets; rj++)
685    jetRecAcc[rj] = jetRecTmp[idxRec[rj]];
686  
687  //rest frame for reconstructed jets
688  for (Int_t i = 0; i < nAccJets; i++)
689    {
690      vRec[i].SetPxPyPzE(jetRecAcc[i].Px(), jetRecAcc[i].Py(), jetRecAcc[i].Pz(), jetRecAcc[i].E());
691      pRec[i].SetXYZ(vRec[i].Px(), vRec[i].Py(), vRec[i].Pz());
692      vsumRec += vRec[i];
693    }
694
695  //check balance of two leading hadrons, deltaPhi > 2.
696  Int_t idxTrack[kTracks];
697  TMath::Sort(nTracks, jetTrackPt, idxTrack);
698
699  TLorentzVector jetTrack[kTracks];
700  for(Int_t iTr = 0; iTr < nTracks; iTr++)
701    jetTrack[iTr] = jetTrackTmp[idxTrack[iTr]];
702  
703  Int_t n = 1;
704  while(jetTrack[0].DeltaPhi(jetTrack[n]) < 2.8)
705    n++;
706
707  Double_t et1 = 0;
708  Double_t et2 = 0;
709  for(Int_t iTr = 0; iTr < nTracks; iTr++)
710    {
711      if(TMath::Abs(jetTrack[0].DeltaPhi(jetTrack[iTr]) < 1.) && iTr != 0)
712        et1 += jetTrack[iTr].Et();
713
714      if(TMath::Abs(jetTrack[n].DeltaPhi(jetTrack[iTr]) < 1.) && iTr != n)
715        et2 += jetTrack[iTr].Et();
716    }
717
718  if(nAccJets == 2)
719    {
720      fhTrackBalance2->Fill(jetTrack[0].Et(), jetTrack[n].Et());
721      fhEt1Et22->Fill(et1, et2);
722    }
723  if(nAccJets == 3)
724    {
725      fhTrackBalance3->Fill(jetTrack[0].Et(), jetTrack[n].Et());
726      fhEt1Et23->Fill(et1, et2);
727    }
728     
729  if(nAccJets > 1 && pRec[0].DeltaPhi(pRec[1]) > 2.8)
730    {
731      fhE2vsE1Rec->Fill(jetRecAcc[0].E(), jetRecAcc[1].E(), fXsection);
732      fhE2E1vsEsumRec->Fill(jetRecAcc[0].E()+jetRecAcc[1].E(), TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
733      fhE2E1vsE1Rec->Fill(jetRecAcc[0].E(), TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
734      Double_t deltaPhi = (jetRecAcc[0].Phi()-jetRecAcc[1].Phi());
735      if(deltaPhi > TMath::Pi()) deltaPhi = deltaPhi - 2.*TMath::Pi();
736      if(deltaPhi < (-1.*TMath::Pi())) deltaPhi = deltaPhi + 2.*TMath::Pi();
737      fhE2E1vsdPhiRec->Fill(-1*deltaPhi, TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
738    }
739  
740  Double_t fPx = vsumRec.Px();
741  Double_t fPy = vsumRec.Py();
742  Double_t fPz = vsumRec.Pz();
743  Double_t fE = vsumRec.E();
744
745  Double_t eSumRecRest = 0;
746  for (Int_t j = 0; j < nAccJets; j++)
747    {
748      vRec[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
749      eRecRest[j] = vRec[j].E();
750      if(nAccJets > 1)
751        fhERecRest->Fill(eRecRest[j], fXsection);
752      eSumRecRest += eRecRest[j];
753    }
754  if(nAccJets > 1)
755    fhEsumRecRest->Fill(eSumRecRest, fXsection);
756  
757  // Relate the jets
758  Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
759  Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
760  
761  for(int i = 0;i<kMaxJets;++i){
762    iGenIndex[i] = iRecIndex[i] = -1;
763  }
764  Double_t dR = 1.4;
765  if(nAccJets == nGenSel)
766    {
767     AliAnalysisHelperJetTasks::GetClosestJets(jetGen,nGenSel,jetRecAcc,nAccJets,
768                                                  iGenIndex,iRecIndex,0);
769       
770      for(int ir = 0;ir < nAccJets;ir++){
771        Int_t ig = iGenIndex[ir];
772        if(ig>=0&&ig<nGenSel){
773          Double_t dPhi = TMath::Abs(jetRecAcc[ir].Phi()-jetGen[ig].Phi());
774          if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi();
775          if(dPhi < (-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi();
776          Double_t sigma = TMath::Abs(jetGen[ig].E()-jetRecAcc[ir].E())/jetGen[ig].E();
777          dR = jetRecAcc[ir].DeltaR(&jetGen[ig]);
778          if(dR < 2*fR && dR >= fR) 
779            {
780              switch(nAccJets)
781                {
782                case 1:
783                  {
784                    fhdEvsErec10[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
785                    fhECorrJet10[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
786                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
787                      {
788                        fhdPhidEtaPt10[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
789                        fhdPhidEta10[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
790                      }
791                  }
792                  break;
793                case 2:
794                  {
795                    fhdEvsErec10[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
796                    fhECorrJet10[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
797                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
798                      {
799                        fhdPhidEtaPt10[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
800                        fhdPhidEta10[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
801                      }
802                  }
803                  break;
804                case 3:
805                  {
806                    fhdEvsErec10[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
807                    fhECorrJet10[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
808                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
809                      {
810                        fhdPhidEtaPt10[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
811                        fhdPhidEta10[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
812                      }
813                  }
814                  break;
815                }
816            }
817          if(dR < fR && dR >= 0.1) 
818            {
819              switch(nAccJets)
820                {
821                case 1:
822                  {
823                    fhdEvsErec05[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
824                    fhECorrJet05[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
825                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
826                      {
827                        fhdPhidEtaPt05[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
828                        fhdPhidEta05[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
829                      }
830                  }
831                  break;
832                case 2:
833                  {
834                    fhdEvsErec05[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
835                    fhECorrJet05[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
836                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
837                      {
838                        fhdPhidEtaPt05[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
839                        fhdPhidEta05[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
840                      }
841                  }
842                  break;
843                case 3:
844                  {
845                    fhdEvsErec05[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
846                    fhECorrJet05[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
847                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
848                      {
849                        fhdPhidEtaPt05[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
850                        fhdPhidEta05[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
851                      }
852                  }
853                  break;
854                }
855            }
856          if(dR < 0.1 && dR >= 0.01)
857            {
858              switch(nAccJets)
859                {
860                case 1:
861                  {
862                    fhdEvsErec01[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
863                    fhECorrJet01[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
864                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
865                      {
866                        fhdPhidEtaPt01[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
867                        fhdPhidEta01[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
868                      }
869                  }
870                  break;
871                case 2:
872                  {
873                    fhdEvsErec01[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
874                    fhECorrJet01[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
875                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
876                      {
877                        fhdPhidEtaPt01[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
878                        fhdPhidEta01[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
879                      }
880                  }
881                  break;
882                case 3:
883                  {
884                    fhdEvsErec01[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
885                    fhECorrJet01[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
886                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
887                      {
888                        fhdPhidEtaPt01[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
889                        fhdPhidEta01[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
890                      }
891                  }
892                  break;
893                }
894            }
895          if(dR > 0.01) continue;
896          switch(nAccJets)
897            {
898            case 1:
899              {
900                fhECorrJet001[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
901                fhdEvsErec001[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
902                for(Int_t iTr = 0; iTr < nTracks; iTr++)
903                  {
904                    fhdPhidEtaPt001[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
905                    fhdPhidEta001[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
906                  }
907              }
908              break;
909            case 2:
910              {
911                fhECorrJet001[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
912                fhdEvsErec001[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
913                for(Int_t iTr = 0; iTr < nTracks; iTr++)
914                  {
915                    fhdPhidEtaPt001[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
916                    fhdPhidEta001[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
917                  }
918              }
919              break;
920            case 3:
921              {
922                fhECorrJet001[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
923                fhdEvsErec001[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
924                for(Int_t iTr = 0; iTr < nTracks; iTr++)
925                  {
926                    fhdPhidEtaPt001[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
927                    fhdPhidEta001[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
928                  }
929              }
930              break;
931            }
932        }
933      }
934      // loop over reconstructed jets
935    }
936
937  Printf("%s:%d",(char*)__FILE__,__LINE__);
938  
939  PostData(1, fList);
940  
941  Printf("%s:%d Data Posted",(char*)__FILE__,__LINE__);
942   
943 }
944
945
946 //__________________________________________________________________________________________________________________________________________________
947 void AliAnalysisTaskJetCorrections::Terminate(Option_t *)
948 {
949   printf("AnalysisJetCorrections::Terminate()");
950
951
952
953 //_______________________________________User defined functions_____________________________________________________________________________________
954
955