]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetCorrections.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliACORDEQAChecker::Check fixed
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetCorrections.cxx
CommitLineData
95392a57 1#include <TROOT.h>
2#include <TSystem.h>
3#include <TInterpreter.h>
4#include <TChain.h>
5#include <TFile.h>
6#include <TH1F.h>
95392a57 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>
95392a57 24
25#include "AliAnalysisTaskJetCorrections.h"
26#include "AliAnalysisManager.h"
95392a57 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"
95392a57 37#include "AliGenCocktailEventHeader.h"
95392a57 38
95392a57 39
40#include "AliAnalysisHelperJetTasks.h"
41
393afcc4 42//
43//
44// corrections to jet energy by sona
45//
46//
95392a57 47
48ClassImp(AliAnalysisTaskJetCorrections)
49
393afcc4 50 AliAnalysisTaskJetCorrections::AliAnalysisTaskJetCorrections() : AliAnalysisTaskSE(),
95392a57 51 fAOD(0x0),
393afcc4 52
95392a57 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{
393afcc4 89 //
90 // ctor
91 //
95392a57 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
113AliAnalysisTaskJetCorrections::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{
393afcc4 153 //
154 // ctor
155 //
95392a57 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
180Bool_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//___________________________________________________________________________________________________________________________________
226void 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//__________________________________________________________________________________________________________________________________________
427void AliAnalysisTaskJetCorrections::Init()
428{
429 printf("AliAnalysisJetCut::Init() \n");
430}
431
432//____________________________________________________________________________________________________________________________________________
433void 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//__________________________________________________________________________________________________________________________________________________
957void AliAnalysisTaskJetCorrections::Terminate(Option_t *)
958{
959 printf("AnalysisJetCorrelation::Terminate()");
960
961}
962
963//_______________________________________User defined functions_____________________________________________________________________________________
393afcc4 964void AliAnalysisTaskJetCorrections::GetThrustAxis(TVector3 &n01, TVector3 * pTrack,const Int_t &nTracks)
95392a57 965{
393afcc4 966 //
967 // Getthrust axis
968 //
95392a57 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