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