]> 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
Add the efficiency estimation to common FF task.
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetCorrections.cxx
CommitLineData
47dcf196 1// **************************************
2// Task used for estimating a charged to neutral correction
6f3f79de 3// sona.pochybova@cern.ch
47dcf196 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
95392a57 23#include <TFile.h>
24#include <TH1F.h>
95392a57 25#include <TH2F.h>
95392a57 26#include <TList.h>
27#include <TTree.h>
95392a57 28#include <TLorentzVector.h>
29#include <TClonesArray.h>
95392a57 30#include <TRefArray.h>
95392a57 31#include <TVector3.h>
95392a57 32#include <TProfile.h>
95392a57 33
34#include "AliAnalysisTaskJetCorrections.h"
35#include "AliAnalysisManager.h"
95392a57 36#include "AliAODEvent.h"
37#include "AliAODVertex.h"
38#include "AliAODHandler.h"
39#include "AliAODTrack.h"
40#include "AliAODJet.h"
e946cd3a 41//#include "AliMCEvent.h"
95392a57 42
43#include "AliAnalysisHelperJetTasks.h"
44
393afcc4 45//
46//
47// corrections to jet energy by sona
48//
49//
95392a57 50
51ClassImp(AliAnalysisTaskJetCorrections)
52
393afcc4 53 AliAnalysisTaskJetCorrections::AliAnalysisTaskJetCorrections() : AliAnalysisTaskSE(),
95392a57 54 fAOD(0x0),
393afcc4 55
95392a57 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{
393afcc4 92 //
93 // ctor
94 //
95392a57 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
116AliAnalysisTaskJetCorrections::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{
393afcc4 156 //
157 // ctor
158 //
95392a57 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
183Bool_t AliAnalysisTaskJetCorrections::Notify()
184{
185 //
186 // Implemented Notify() to read the cross sections
187 // and number of trials from pyxsec.root
188 //
e946cd3a 189
190
95392a57 191 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
e946cd3a 192 Float_t xsection = 0;
193 Float_t ftrials = 1;
194
195 Float_t fAvgTrials = 1;
95392a57 196 if(tree){
197 TFile *curfile = tree->GetCurrentFile();
198 if (!curfile) {
199 Error("Notify","No current file");
200 return kFALSE;
201 }
e946cd3a 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;
95392a57 209
95392a57 210 return kTRUE;
e946cd3a 211
95392a57 212}
213
214
215//___________________________________________________________________________________________________________________________________
216void AliAnalysisTaskJetCorrections::UserCreateOutputObjects()
217{
218 //
219 // Create the output container
220 //
47dcf196 221 // Printf("Analysing event %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
95392a57 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//__________________________________________________________________________________________________________________________________________
417void AliAnalysisTaskJetCorrections::Init()
418{
419 printf("AliAnalysisJetCut::Init() \n");
420}
421
422//____________________________________________________________________________________________________________________________________________
423void 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
e946cd3a 437// AliMCEvent* mcEvent =MCEvent();
438// if(!mcEvent){
439// Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
440// return;
441// }
95392a57 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//__________________________________________________________________________________________________________________________________________________
947void AliAnalysisTaskJetCorrections::Terminate(Option_t *)
948{
e946cd3a 949 printf("AnalysisJetCorrections::Terminate()");
95392a57 950
951}
952
953//_______________________________________User defined functions_____________________________________________________________________________________
47dcf196 954
95392a57 955