AliACORDEQAChecker::Check fixed
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskThreeJets.cxx
CommitLineData
95392a57 1#include <TSystem.h>
95392a57 2#include <TFile.h>
3#include <TH1F.h>
95392a57 4#include <TH2F.h>
5#include <TH3F.h>
6#include <TList.h>
7#include <TTree.h>
95392a57 8#include <TLorentzVector.h>
9#include <TClonesArray.h>
95392a57 10#include <TRefArray.h>
95392a57 11#include <TVector3.h>
12#include <TVectorD.h>
13#include <TMatrixDSym.h>
14#include <TMatrixDSymEigen.h>
95392a57 15#include <TProfile.h>
95392a57 16
17#include "AliAnalysisTaskThreeJets.h"
18#include "AliAnalysisManager.h"
95392a57 19#include "AliAODEvent.h"
20#include "AliAODVertex.h"
21#include "AliAODHandler.h"
22#include "AliAODTrack.h"
23#include "AliAODJet.h"
393afcc4 24#include "AliGenPythiaEventHeader.h"
95392a57 25#include "AliMCEventHandler.h"
26#include "AliMCEvent.h"
27#include "AliStack.h"
95392a57 28
95392a57 29
30#include "AliAnalysisHelperJetTasks.h"
31
393afcc4 32//
33//
34// Trhreee jet task by sona
35//
36//
95392a57 37
38ClassImp(AliAnalysisTaskThreeJets)
39
40 AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets() : AliAnalysisTaskSE(),
41
42 fAOD(0x0),
43
44 fBranchRec(""),
45 fBranchGen(""),
46
47 fUseAODInput(kFALSE),
48
49 fR(0x0),
50 fList(0x0),
51
52 fGlobVar(1),
53 fXsection(1),
54
55 fhX3X4Rec(0x0),
56 fhX3X4Gen(0x0),
57
58 fhMu35Rec(0x0),
59 fhMu34Rec(0x0),
60 fhMu45Rec(0x0),
61
62 fhMu35Gen(0x0),
63 fhMu34Gen(0x0),
64 fhMu45Gen(0x0),
65
66 fhInOut(0x0),
67
68 fhThrustRec2(0x0),
69 fhThrustRec3(0x0),
70
71 fhThrustGen2(0x0),
72 fhThrustGen3(0x0),
73
74 fhCGen2(0x0),
75 fhCGen3(0x0),
76
77 fhSGen2(0x0),
78 fhSGen3(0x0),
79
80 fhAGen2(0x0),
81 fhAGen3(0x0),
82
83 fhCRec2(0x0),
84 fhCRec3(0x0),
85
86 fhSRec2(0x0),
87 fhSRec3(0x0),
88
89 fhARec2(0x0),
90 fhARec3(0x0),
91
92 fhX3(0x0),
93 fhX4(0x0),
94 fhX5(0x0),
95
96 fhXSec(0x0),
97
98 fhX3X4Rec60(0x0),
99 fhX3X4Rec60100(0x0),
100 fhX3X4Rec100(0x0),
101
102 fhX3X4Gen60(0x0),
103 fhX3X4Gen60100(0x0),
104 fhX3X4Gen100(0x0),
105
106 fhdPhiThrustGen(0x0),
107 fhdPhiThrustGenALL(0x0),
108
109 fhdPhiThrustRec(0x0),
110 fhdPhiThrustRecALL(0x0)
111{
112
113}
114
115AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets(const char * name):
116 AliAnalysisTaskSE(name),
117
118 fAOD(0x0),
119
120 fBranchRec(""),
121 fBranchGen(""),
122
123 fUseAODInput(kFALSE),
124
125 fR(0x0),
126 fList(0x0),
127
128 fGlobVar(1),
129 fXsection(1),
130
131 fhX3X4Rec(0x0),
132 fhX3X4Gen(0x0),
133
134 fhMu35Rec(0x0),
135 fhMu34Rec(0x0),
136 fhMu45Rec(0x0),
137
138 fhMu35Gen(0x0),
139 fhMu34Gen(0x0),
140 fhMu45Gen(0x0),
141
142 fhInOut(0x0),
143
144 fhThrustRec2(0x0),
145 fhThrustRec3(0x0),
146
147 fhThrustGen2(0x0),
148 fhThrustGen3(0x0),
149
150 fhCGen2(0x0),
151 fhCGen3(0x0),
152
153 fhSGen2(0x0),
154 fhSGen3(0x0),
155
156 fhAGen2(0x0),
157 fhAGen3(0x0),
158
159 fhCRec2(0x0),
160 fhCRec3(0x0),
161
162 fhSRec2(0x0),
163 fhSRec3(0x0),
164
165 fhARec2(0x0),
166 fhARec3(0x0),
167
168 fhX3(0x0),
169 fhX4(0x0),
170 fhX5(0x0),
171
172 fhXSec(0x0),
173
174 fhX3X4Rec60(0x0),
175 fhX3X4Rec60100(0x0),
176 fhX3X4Rec100(0x0),
177
178 fhX3X4Gen60(0x0),
179 fhX3X4Gen60100(0x0),
180 fhX3X4Gen100(0x0),
181
182 fhdPhiThrustGen(0x0),
183 fhdPhiThrustGenALL(0x0),
184
185 fhdPhiThrustRec(0x0),
186 fhdPhiThrustRecALL(0x0)
187{
188 DefineOutput(1, TList::Class());
189}
190
191
192
193Bool_t AliAnalysisTaskThreeJets::Notify()
194{
195 //
196 // Implemented Notify() to read the cross sections
197 // and number of trials from pyxsec.root
198 //
199 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
200 UInt_t ntrials = 0;
201 if(tree){
202 TFile *curfile = tree->GetCurrentFile();
203 if (!curfile) {
204 Error("Notify","No current file");
205 return kFALSE;
206 }
207
208 TString fileName(curfile->GetName());
209 if(fileName.Contains("AliESDs.root")){
210 fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
211 }
212 else if(fileName.Contains("AliAOD.root")){
213 fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
214 }
215 else if(fileName.Contains("galice.root")){
216 // for running with galice and kinematics alone...
217 fileName.ReplaceAll("galice.root", "pyxsec.root");
218 }
219 TFile *fxsec = TFile::Open(fileName.Data());
220 if(!fxsec){
221 Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
222 // no a severe condition
223 return kTRUE;
224 }
225 TTree *xtree = (TTree*)fxsec->Get("Xsection");
226 if(!xtree){
227 Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
228 }
229 xtree->SetBranchAddress("xsection",&fXsection);
230 xtree->SetBranchAddress("ntrials",&ntrials);
231 xtree->GetEntry(0);
232 }
233
234 return kTRUE;
235}
236
237
238//___________________________________________________________________________________________________________________________________
239void AliAnalysisTaskThreeJets::UserCreateOutputObjects()
240{
241 //
242 // Create the output container
243 //
244 Printf("Analysing event %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
245
246 if(fUseAODInput){
247 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
248 if(!fAOD){
249 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
250 return;
251 }
252 }
253 else{
254 // assume that the AOD is in the general output...
255 fAOD = AODEvent();
256 if(!fAOD){
257 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
258 return;
259 }
260 }
261
262 printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
263
264 fList = new TList();
265
266 fhX3X4Gen = new TH2F("X3vsX4Gen", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
267 fhX3X4Gen->SetXTitle("X_{3}");
268 fhX3X4Gen->SetYTitle("X_{4}");
269 fhX3X4Gen->Sumw2();
270 fList->Add(fhX3X4Gen);
271
272 fhX3X4Rec = new TH2F("X3vsX4Rec", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
273 fhX3X4Rec->SetXTitle("X_{3}");
274 fhX3X4Rec->SetYTitle("X_{4}");
275 fhX3X4Rec->Sumw2();
276 fList->Add(fhX3X4Rec);
277
278 fhMu35Rec = new TH1F("Mu35Rec", "", 20,0.1,0.8);
279 fhMu35Rec->Sumw2();
280 fhMu35Rec->SetXTitle("#mu_{35}");
281 fhMu35Rec->SetYTitle("#frac{dN}{d#mu_{35Rec}}");
282 fList->Add(fhMu35Rec);
283
284 fhMu34Rec = new TH1F("Mu34Rec", "", 20,0.5,1);
285 fhMu34Rec->Sumw2();
286 fhMu34Rec->SetXTitle("#mu_{34}");
287 fhMu34Rec->SetYTitle("#frac{dN}{d#mu_{34}}");
288 fList->Add(fhMu34Rec);
289
290 fhMu45Rec = new TH1F("Mu45Rec", "", 20,0,0.65);
291 fhMu45Rec->Sumw2();
292 fhMu45Rec->SetXTitle("#mu_{45}");
293 fhMu45Rec->SetYTitle("#frac{dN}{d#mu_{45}}");
294 fList->Add(fhMu45Rec);
295
296 fhMu35Gen = new TH1F("Mu35Gen", "", 20,0.1,0.8);
297 fhMu35Gen->Sumw2();
298 fhMu35Gen->SetXTitle("#mu_{35Gen}");
299 fhMu35Gen->SetYTitle("#frac{dN}{d#mu_{35Gen}}");
300 fList->Add(fhMu35Gen);
301
302 fhMu34Gen = new TH1F("Mu34Gen", "", 20,0.5,1);
303 fhMu34Gen->Sumw2();
304 fhMu34Gen->SetXTitle("#mu_{34Gen}");
305 fhMu34Gen->SetYTitle("#frac{dN}{d#mu_{34Gen}}");
306 fList->Add(fhMu34Gen);
307
308 fhMu45Gen = new TH1F("Mu45Gen", "", 20,0,0.65);
309 fhMu45Gen->Sumw2();
310 fhMu45Gen->SetXTitle("#mu_{45Gen}");
311 fhMu45Gen->SetYTitle("#frac{dN}{d#mu_{45}}");
312 fList->Add(fhMu45Gen);
313
314 fhInOut = new TH1I("InOut", "", 6, 0, 6);
315 fhInOut->SetXTitle("#RecJets_{GenJets=3}");
316 fhInOut->SetYTitle("#Entries");
317 fList->Add(fhInOut);
318
319 fhThrustGen2 = new TH1F("ThrustGen2", "", 50, 0.5, 1);
320 fhThrustGen2->Sumw2();
321 fList->Add(fhThrustGen2);
322
323 fhThrustGen3 = new TH1F("ThrustGen3", "", 50, 0.5, 1);
324 fhThrustGen3->Sumw2();
325 fList->Add(fhThrustGen3);
326
327 fhThrustRec2 = new TH1F("ThrustRec2", "", 50, 0.5, 1);
328 fhThrustRec2->Sumw2();
329 fList->Add(fhThrustRec2);
330
331 fhThrustRec3 = new TH1F("ThrustRec3", "", 50, 0.5, 1);
332 fhThrustRec3->Sumw2();
333 fList->Add(fhThrustRec3);
334
335 fhCGen2 = new TH1F("CGen2", "", 100, 0, 1);
336 fhCGen2->Sumw2();
337 fList->Add(fhCGen2);
338
339 fhCGen3 = new TH1F("CGen3", "", 100, 0, 1);
340 fhCGen3->Sumw2();
341 fList->Add(fhCGen3);
342
343 fhCRec2 = new TH1F("CRec2", "", 100, 0, 1);
344 fhCRec2->Sumw2();
345 fList->Add(fhCRec2);
346
347 fhCRec3 = new TH1F("CRec3", "", 100, 0, 1);
348 fhCRec3->Sumw2();
349 fList->Add(fhCRec3);
350
351 fhSGen2 = new TH1F("SGen2", "", 100, 0, 1);
352 fList->Add(fhSGen2);
353
354 fhSGen3 = new TH1F("SGen3", "", 100, 0, 1);
355 fList->Add(fhSGen3);
356
357 fhSRec2 = new TH1F("SRec2", "", 100, 0, 1);
358 fList->Add(fhSRec2);
359
360 fhSRec3 = new TH1F("SRec3", "", 100, 0, 1);
361 fList->Add(fhSRec3);
362
363 fhAGen2 = new TH1F("AGen2", "", 50, 0, 0.5);
364 fList->Add(fhAGen2);
365
366 fhAGen3 = new TH1F("AGen3", "", 50, 0, 0.5);
367 fList->Add(fhAGen3);
368
369 fhARec2 = new TH1F("ARec2", "", 50, 0, 0.5);
370 fList->Add(fhARec2);
371
372 fhARec3 = new TH1F("ARec3", "", 50, 0, 0.5);
373 fList->Add(fhARec3);
374
375 fhX3 = new TH2F("X3", "", 22, 0.6, 1.02, 100, 0, 1);
376 fhX3->SetYTitle("|X_{3}^{MC} - X_{3}^{AOD}|/X_{3}^{MC}");
377 fhX3->SetXTitle("X_{3}");
378 fhX3->Sumw2();
379 fList->Add(fhX3);
380
381 fhX4 = new TH2F("X4", "",33, 0.4, 1.02, 100, 0, 1);
382 fhX4->SetYTitle("|X_{4}^{MC} - X_{4}^{AOD}|/X_{4}^{MC}");
383 fhX4->SetXTitle("X_{4}");
384 fhX4->Sumw2();
385 fList->Add(fhX4);
386
387 fhX5 = new TH2F("X5", "",100, 0., 1., 100, 0, 1);
388 fhX5->SetYTitle("|X_{5}^{MC} - X_{5}^{AOD}|/X_{5}^{MC}");
389 fhX5->SetXTitle("X_{5}");
390 fhX5->Sumw2();
391 fList->Add(fhX5);
392
393 fhXSec = new TProfile("XSec", "", 200, 0, 200, 0, 1);
394 fList->Add(fhXSec);
395
396 fhX3X4Rec60 = new TH2F("X3vsX4Rec60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
397 fhX3X4Rec60->SetXTitle("X_{3}");
398 fhX3X4Rec60->SetYTitle("X_{4}");
399 fhX3X4Rec60->Sumw2();
400 fList->Add(fhX3X4Rec60);
401
402 fhX3X4Rec60100 = new TH2F("X3vsX4Rec60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
403 fhX3X4Rec60100->SetXTitle("X_{3}");
404 fhX3X4Rec60100->SetYTitle("X_{4}");
405 fhX3X4Rec60100->Sumw2();
406 fList->Add(fhX3X4Rec60100);
407
408 fhX3X4Rec100 = new TH2F("X3vsX4Rec100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
409 fhX3X4Rec100->SetXTitle("X_{3}");
410 fhX3X4Rec100->SetYTitle("X_{4}");
411 fhX3X4Rec100->Sumw2();
412 fList->Add(fhX3X4Rec100);
413
414 fhX3X4Gen60 = new TH2F("X3vsX4Gen60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
415 fhX3X4Gen60->SetXTitle("X_{3}");
416 fhX3X4Gen60->SetYTitle("X_{4}");
417 fhX3X4Gen60->Sumw2();
418 fList->Add(fhX3X4Gen60);
419
420 fhX3X4Gen60100 = new TH2F("X3vsX4Gen60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
421 fhX3X4Gen60100->SetXTitle("X_{3}");
422 fhX3X4Gen60100->SetYTitle("X_{4}");
423 fhX3X4Gen60100->Sumw2();
424 fList->Add(fhX3X4Gen60100);
425
426 fhX3X4Gen100 = new TH2F("X3vsX4Gen100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
427 fhX3X4Gen100->SetXTitle("X_{3}");
428 fhX3X4Gen100->SetYTitle("X_{4}");
429 fhX3X4Gen100->Sumw2();
430 fList->Add(fhX3X4Gen100);
431
432 fhdPhiThrustGen = new TH2F("dPhiThrustGen", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
433 fhdPhiThrustGen->Sumw2();
434 fList->Add(fhdPhiThrustGen);
435
436 fhdPhiThrustGenALL = new TH2F("dPhiThrustGenALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
437 fhdPhiThrustGenALL->Sumw2();
438 fList->Add(fhdPhiThrustGenALL);
439
440 fhdPhiThrustRec = new TH2F("dPhiThrustRec", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
441 fhdPhiThrustRec->Sumw2();
442 fList->Add(fhdPhiThrustRec);
443
444 fhdPhiThrustRecALL = new TH2F("dPhiThrustRecALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
445 fhdPhiThrustRecALL->Sumw2();
446 fList->Add(fhdPhiThrustRecALL);
447
448 Printf("UserCreateOutputObjects finished\n");
449}
450
451//__________________________________________________________________________________________________________________________________________
452void AliAnalysisTaskThreeJets::Init()
453{
454 printf("AliAnalysisJetCut::Init() \n");
455}
456
457//____________________________________________________________________________________________________________________________________________
458void AliAnalysisTaskThreeJets::UserExec(Option_t * )
459{
460 // if (fDebug > 1) printf("Analysing event # %5d\n", (Int_t) fEntry);
461
462
463 //create an AOD handler
464 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
465
466 if(!aodH)
467 {
468 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
469 return;
470 }
471
472 AliMCEvent* mcEvent =MCEvent();
473 if(!mcEvent){
474 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
475 return;
476 }
477
478 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
479
480 //primary vertex
481 AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
482 if(!pvtx) return;
483
484 AliAODJet genJetsPythia[kMaxJets];
485 Int_t nPythiaGenJets = 0;
486
487 AliAODJet genJets[kMaxJets];
488 Int_t nGenJets = 0;
489
490 AliAODJet recJets[kMaxJets];
491 Int_t nRecJets = 0;
492
493 //array of reconstructed jets from the AOD input
494 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
495 if(!aodRecJets){
496 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
497 return;
498 }
499
500 // reconstructed jets
501 nRecJets = aodRecJets->GetEntries();
502 nRecJets = TMath::Min(nRecJets, kMaxJets);
503
504 for(int ir = 0;ir < nRecJets;++ir)
505 {
506 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
507 if(!tmp)continue;
508 recJets[ir] = *tmp;
509 }
510
511 // If we set a second branch for the input jets fetch this
512 TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
513
514 if(!aodGenJets)
515 {
516 printf("NO MC jets Found\n");
517 return;
518 }
519
520 // //Generated jets
521 nGenJets = aodGenJets->GetEntries();
522 nGenJets = TMath::Min(nGenJets, kMaxJets);
523
524 for(Int_t ig =0 ; ig < nGenJets; ++ig)
525 {
526 AliAODJet * tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
527 if(!tmp)continue;
528 genJets[ig] = * tmp;
529 }
530
531 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
532 if(!pythiaGenHeader){
533 Printf("!!!NO GEN HEADER AVALABLE!!!");
534 return;
535 }
536
537 // Int_t ProcessType = pythiaGenHeader->ProcessType();
538 // if(ProcessType != 28) return;
539 fhXSec->Fill(pythiaGenHeader->GetPtHard(), fXsection);
540 nPythiaGenJets = pythiaGenHeader->NTriggerJets();
541 nPythiaGenJets = TMath::Min(nPythiaGenJets, kMaxJets);
542
543 Double_t eRec[kMaxJets];
544 Double_t eGen[kMaxJets];
545
546 Double_t eJetRec[kMaxJets];
547 // Double_t EJetGen[kMaxJets];
548
549 AliAODJet jetRec[kMaxJets];
550 AliAODJet jetGen[kMaxJets];
551
552 Int_t idxRec[kMaxJets];
553 Int_t idxGen[kMaxJets];
554
555 Double_t xRec[kMaxJets];
556 Double_t xGen[kMaxJets];
557
558 Double_t eSumRec = 0;
559 Double_t eSumGen = 0;
560
561 TLorentzVector vRec[kMaxJets];
562 TLorentzVector vRestRec[kMaxJets];
563
564 TLorentzVector vGen[kMaxJets];
565 TLorentzVector vRestGen[kMaxJets];
566
567 TLorentzVector vsumRec;
568 TLorentzVector vsumGen;
569
570 TVector3 pRec[kMaxJets];
571 TVector3 pGen[kMaxJets];
572
573 TVector3 pTrack[kTracks];
574
575 TVector3 pRestRec[kMaxJets];
576 TVector3 pRestGen[kMaxJets];
577
578 Double_t psumRestRec = 0;
579 // Double_t psumRestGen = 0;
580 //Pythia_________________________________________________________________________________________________________________
581
582 for(int ip = 0;ip < nPythiaGenJets;++ip)
583 {
584 if(ip>=kMaxJets)continue;
585 Float_t p[4];
586 pythiaGenHeader->TriggerJet(ip,p);
587 genJetsPythia[ip].SetPxPyPzE(p[0],p[1],p[2],p[3]);
588 }
589
590//_________________________________________________________________________________________________________________________
591
592
593//________histos for MC___________________________________________________________________________________________________________
594
595 Int_t nGenSel = 0;
596 Int_t counter = 0;
597 Int_t tag = 0;
598
599 AliAODJet selJets[kMaxJets];
600
601 for(Int_t i = 0; i < nGenJets; i++)
602 {
603 if(nGenJets == 1)
604 {
605 selJets[nGenSel] = genJets[i];
606 nGenSel++;
607 }
608 else
609 {
610 counter = 0;
611 tag = 0;
612 for(Int_t j = 0; j < nGenJets; j++)
613 {
614 if(i!=j)
615 {
616 Double_t dRij = genJets[i].DeltaR(&genJets[j]);
617 counter++;
618 if(dRij > 2*fR) tag++;
619 }
620 }
621 if(counter!=0)
622 {
623 if(tag/counter == 1)
624 {
625 selJets[nGenSel] = genJets[i];
626 nGenSel++;
627 }
628 }
629 }
630 }
631
632 if(nGenSel == 0) return;
633
634 for (Int_t gj = 0; gj < nGenSel; gj++)
635 {
636 eGen[gj] = selJets[gj].E();
637 }
638
639 TMath::Sort(nGenSel, eGen, idxGen);
640 for (Int_t ig = 0; ig < nGenSel; ig++)
641 {
642 jetGen[ig] = selJets[idxGen[ig]];
643 }
644 AliStack * stack = mcEvent->Stack();
645 Int_t nMCtracks = stack->GetNprimary();
646 Double_t * eTracksMC = new Double_t[kTracks];
647 Double_t pTracksMC[kTracks];
648 Int_t * idxTracksMC = new Int_t[kTracks];
649 TLorentzVector jetTracksMC[kTracks];
650 TLorentzVector jetTracksSortMC[kTracks];
651 TVector3 pTrackMC[kTracks];
652 TLorentzVector vTrackMCAll[kTracks];
653 Double_t pTrackMCAll[kTracks];
654 TLorentzVector vTrackMC[kTracks];
655 TVector3 pTrackMCBoost[kTracks];
656 Double_t eventShapes[4];
657
658 Int_t nAccTr = 0;
659 Int_t nInJet[kMaxJets];
660 TLorentzVector inJetPartV[kMaxJets][kTracks];
661 Int_t nAllTracksMC = 0;
662
663 for(Int_t iTrack = 0; iTrack < nMCtracks; iTrack++)
664 {
665 TParticle * part = (TParticle*)stack->Particle(iTrack);
666 if (!part) continue;
667 Double_t fEta = part->Eta();
668 if(TMath::Abs(fEta) > .9) continue;
669 if(!IsPrimChar(part, nMCtracks, 0)) continue;
670 vTrackMCAll[nAllTracksMC].SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
671 pTrackMCAll[nAllTracksMC] = part->Pt();
672 nAllTracksMC++;
673 }
674 if(nAllTracksMC == 0) return;
675 for(Int_t iJet = 0; iJet < nGenSel; iJet++)
676 {
677 Int_t nJetTracks = 0;
678 for(Int_t i = 0; i < nAllTracksMC; i++)
679 {
680 Double_t dPhi = (jetGen[iJet].Phi()-vTrackMCAll[i].Phi());
681 if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi();
682 if(dPhi < (-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi();
683 Double_t dEta = (jetGen[iJet].Eta()-vTrackMCAll[i].Eta());
684 Double_t deltaR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
685 if(deltaR < fR && vTrackMCAll[i].Pt() > 1.5)
686 {
687 jetTracksMC[nAccTr] = vTrackMCAll[i];
688 eTracksMC[nAccTr] = vTrackMCAll[i].E();
689 pTracksMC[nAccTr] = vTrackMCAll[i].Pt();
690 inJetPartV[iJet][nJetTracks].SetPxPyPzE(vTrackMCAll[i].Px(), vTrackMCAll[i].Py(), vTrackMCAll[i].Pz(),vTrackMCAll[i].E());
691 nAccTr++;
692 nJetTracks++;
693 }
694 }
695 nInJet[iJet] = nJetTracks;
696 }
697
698 if(nAccTr == 0) return;
699 Printf("*********** Number of Jets : %d ***************\n", nGenSel);
700 Double_t pTav[kMaxJets];
701 for(Int_t i = 0; i < nGenSel; i++)
702 {
703 Double_t pTsum = 0;
704 Printf("*********** Number of particles in Jet %d = %d *******************\n", i+3, nInJet[i]);
705 for(Int_t iT = 0; iT < nInJet[i]; iT++)
706 {
707 Double_t pt = inJetPartV[i][iT].Pt();
708 pTsum += pt;
709 }
710 pTav[i] = pTsum/nInJet[i];
711 }
712
713 TMath::Sort(nAccTr, pTracksMC, idxTracksMC);
714 for(Int_t i = 0; i < nAccTr; i++)
715 {
716 jetTracksSortMC[i] = jetTracksMC[idxTracksMC[i]];
717 pTrackMC[i].SetXYZ(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz());
718 vTrackMC[i].SetPxPyPzE(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz(), jetTracksSortMC[i].E());
719 }
720
721 TVector3 n01MC = pTrackMC[0].Unit();
722 //Thrust calculation, iterative method
723 if(nGenSel > 1)
724 {
725// if(fGlobVar == 1)
726// {
727 GetEventShapes(n01MC, pTrackMC, nAccTr, eventShapes);
728// }
729
730 Double_t s = eventShapes[1];
731 Double_t a = eventShapes[2];
732 Double_t c = eventShapes[3];
733
734 switch(nGenSel)
735 {
736 case 2:
737 {
738 fhAGen2->Fill(a);
739 fhSGen2->Fill(s);
740 fhCGen2->Fill(c);
741 }
742 break;
743 case 3:
744 {
745 fhAGen3->Fill(a);
746 fhSGen3->Fill(s);
747 fhCGen3->Fill(c);
748 }
749 break;
750 }
751 Double_t thrust01MC = eventShapes[0];
752
753 switch(nGenSel)
754 {
755 case 2:
756 fhThrustGen2->Fill(thrust01MC, fXsection);
757 break;
758 case 3:
759 fhThrustGen3->Fill(thrust01MC, fXsection);
760 break;
761 }
762 }
763
764
765 //rest frame MC jets
766 for (Int_t i = 0; i < nGenSel; ++i)
767 {
768 vGen[i].SetPxPyPzE(jetGen[i].Px(), jetGen[i].Py(), jetGen[i].Pz(), jetGen[i].E());
769 pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
770 vsumGen += vGen[i];
771 }
772 if(eventShapes[0] >0.8 )
773 {
774 for(Int_t i = 0; i < nGenSel; i++)
775 fhdPhiThrustGen->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
776 }
777
778 if(eventShapes[0] <= 0.8)
779 {
780 for(Int_t i = 0; i < nGenSel; i++)
781 fhdPhiThrustGenALL->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
782 }
783
784 Double_t fPxGen = vsumGen.Px();
785 Double_t fPyGen = vsumGen.Py();
786 Double_t fPzGen = vsumGen.Pz();
787 Double_t fEGen = vsumGen.E();
788
789 Double_t eRestGen[kMaxJets];
790 for (Int_t j = 0; j < nGenSel; j++)
791 {
792 vGen[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
793 eRestGen[j] = vGen[j].E();
794 }
795
796 for (Int_t j = 0; j < nAccTr; j++)
797 {
798 vTrackMC[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
799 pTrackMCBoost[j].SetXYZ(vTrackMC[j].Px(),vTrackMC[j].Py(),vTrackMC[j].Pz());
800 }
801
802 Int_t idxRestGen[kMaxJets];
803 TMath::Sort(nGenSel, eRestGen, idxRestGen);
804 for(Int_t j = 0; j < nGenSel; j++)
805 {
806 vRestGen[j] = vGen[idxRestGen[j]];
807 eSumGen += vRestGen[j].E();
808 }
809
810 if (nGenSel == 3)
811 {
812 // if(nInJet[0] < 3 || nInJet[1] < 3 || nInJet[2] < 3) return;
813 // if(pRestGen[1].DeltaPhi(pRestGen[2]) > 0.95 && pRestGen[1].DeltaPhi(pRestGen[2]) < 1.15)
814 // {
815
816 for(Int_t i = 0; i < nGenSel; i++)
817 {
818 xGen[i] = 2*vRestGen[i].E()/eSumGen;
819 }
820
821 Printf("***************** Values of Dalitz variables are : %f, %f, %f ****************\n", xGen[0], xGen[1], xGen[2]);
822
823 Printf("***************** fXSection = %f ******************\n", fXsection);
824 if(eSumGen <= 60)
825 fhX3X4Gen60->Fill(xGen[0], xGen[1], fXsection);
826
827 if(eSumGen > 60 && eSumGen <= 100)
828 fhX3X4Gen60100->Fill(xGen[0], xGen[1], fXsection);
829
830 if(eSumGen > 100)
831 fhX3X4Gen100->Fill(xGen[0], xGen[1], fXsection);
832
833 FillTopology(fhX3X4Gen, fhMu34Gen, fhMu45Gen, fhMu35Gen, xGen, pRestGen, fXsection);
834 }
835
836
837
838//_______________________________________________histos for MC_____________________________________________________
839
840
841//_______________________________________________histos AOD________________________________________________________
842
843// Printf("Event Number : %d, Number of gen jets : %d ", fEntry, nGenJets);
844
845 Int_t nRecSel = 0;
846 Int_t counter1 = 0;
847 Int_t tag1 = 0;
848
849 AliAODJet recSelJets[kMaxJets];
850
851 for(Int_t i = 0; i < nRecJets; i++)
852 {
853 if(nRecJets == 1)
854 {
855 recSelJets[nRecSel] = recJets[i];
856 nRecSel++;
857 }
858 else
859 {
860 counter1 = 0;
861 tag1 = 0;
862 for(Int_t j = 0; j < nRecJets; j++)
863 {
864 if(i!=j)
865 {
866 Double_t dRij = recJets[i].DeltaR(&recJets[j]);
867 counter1++;
868 if(dRij > 2*fR) tag1++;
869 }
870 }
871 if(counter1!=0)
872 {
873 if(tag1/counter1 == 1)
874 {
875 recSelJets[nRecSel] = recJets[i];
876 nRecSel++;
877 }
878 }
879 }
880 }
881
882 if(nRecSel == 0) return;
883
884 //sort rec/gen jets by energy in C.M.S
885 for (Int_t rj = 0; rj < nRecSel; rj++)
886 {
887 eRec[rj] = recSelJets[rj].E();
888 }
889
890 // Int_t nAODtracks = fAOD->GetNumberOfTracks();
891 Int_t nTracks = 0; //tracks accepted in the whole event
892 // Int_t nTracksALL = 0;
893 TLorentzVector jetTracks[kTracks];
894 TLorentzVector jetTracksSort[kTracks];
895 Double_t * eTracks = new Double_t[kTracks];
896 Double_t pTracks[kTracks];
897 Int_t * idxTracks = new Int_t[kTracks];
898 Double_t eventShapesRec[4];
899 Int_t jetMult[kMaxJets];
900 // TLorentzVector vTracksAll[kTracks];
901 // Double_t pTracksAll[kTracks];
902 Int_t nAccJets = 0;
903 AliAODJet jetRecAcc[kMaxJets];
904 Int_t nJetTracks = 0;
905
906 AliAODTrack jetTrack[kTracks];
907 Double_t * cv = new Double_t[21];
908 TMath::Sort(nRecSel, eRec, idxRec);
909 for (Int_t rj = 0; rj < nRecSel; rj++)
910 {
911 nJetTracks = 0;
912 eJetRec[rj] = eRec[idxRec[rj]];
913 jetRec[rj] = recSelJets[idxRec[rj]];
914 TRefArray * jetTracksAOD = dynamic_cast<TRefArray*>(jetRec[rj].GetRefTracks());
915 if(!jetTracksAOD) continue;
916 if(jetTracksAOD->GetEntries() < 3) continue;
917 for(Int_t i = 0; i < jetTracksAOD->GetEntries(); i++)
918 {
919 AliAODTrack * track = (AliAODTrack*)jetTracksAOD->At(i);
920 track->GetCovarianceXYZPxPyPz(cv);
921 if(cv[14] > 1000.) continue;
922 jetTrack[nTracks] = *track;
923 jetTracks[nTracks].SetPxPyPzE(jetTrack[nTracks].Px(), jetTrack[nTracks].Py(), jetTrack[nTracks].Pz(), jetTrack[nTracks].E());
924 eTracks[nTracks] = jetTracks[nTracks].E();
925 pTracks[nTracks] = jetTracks[nTracks].Pt();
926 nTracks++;
927 nJetTracks++;
928 }
929 if(nJetTracks < 3) continue;
930 jetRecAcc[nAccJets] = jetRec[rj];
931 jetMult[nAccJets] = jetTracksAOD->GetEntries();
932 nAccJets++;
933 }
934
935 if (nAccJets == 0) return;
936
937 TLorentzVector vTrack[kTracks];
938 TMath::Sort(nTracks, pTracks, idxTracks);
939 for(Int_t i = 0; i < nTracks; i++)
940 {
941 jetTracksSort[i] = jetTracks[idxTracks[i]];
942 pTrack[i].SetXYZ(jetTracksSort[i].Px(), jetTracksSort[i].Py(), jetTracksSort[i].Pz());
943 vTrack[i].SetPxPyPzE(jetTracksSort[i].Px(), jetTracksSort[i].Py(), jetTracksSort[i].Pz(), jetTracksSort[i].E());
944 }
945
946 for (Int_t i = 0; i < nAccJets; ++i)
947 {
948 vRec[i].SetPxPyPzE(jetRecAcc[i].Px(), jetRecAcc[i].Py(), jetRecAcc[i].Pz(), jetRecAcc[i].E());
949 pRec[i].SetXYZ(vRec[i].Px(), vRec[i].Py(), vRec[i].Pz());
950 vsumRec += vRec[i];
951 }
952
953 //Thrust, iterative method, AODs
954 TVector3 n01 = pTrack[0].Unit();
955 if(nAccJets > 1)
956 {
957// if(fGlobVar == 1)
958// {
959 GetEventShapes(n01, pTrack, nTracks, eventShapesRec);
960// }
961// fGlobVar = 0;
962// Double_t Max3 = TMath::Max(eventShapesRec0[0],eventShapesRec1[0]);
963// Double_t Max4 = TMath::Max(eventShapesRec3[0],eventShapesRec2[0]);
964
965 Double_t thrust = eventShapesRec[0];
966
967 if(eventShapesRec[0] > 0.8)
968 {
969 for(Int_t i = 0; i < nAccJets; i++)
970 fhdPhiThrustRec->Fill(n01.DeltaPhi(pRec[i]), jetRecAcc[i].E());
971
972 }
973
974 if(eventShapesRec[0] <= 0.8)
975 {
976 for(Int_t i = 0; i < nAccJets; i++)
977 fhdPhiThrustRecALL->Fill(n01.DeltaPhi(pRec[i]), jetRecAcc[i].E());
978 }
979
980 switch(nAccJets)
981 {
982 case 2:
983 fhThrustRec2->Fill(thrust, fXsection);
984 break;
985 case 3:
986 fhThrustRec3->Fill(thrust, fXsection);
987 break;
988 }
989
990 switch(nAccJets)
991 {
992 case 2:
993 {
994 fhARec2->Fill(eventShapesRec[2], fXsection);
995 fhSRec2->Fill(eventShapesRec[1], fXsection);
996 fhCRec2->Fill(eventShapesRec[3], fXsection);
997 }
998 break;
999 case 3:
1000 {
1001 fhARec3->Fill(eventShapesRec[2], fXsection);
1002 fhSRec3->Fill(eventShapesRec[1], fXsection);
1003 fhCRec3->Fill(eventShapesRec[3], fXsection);
1004 }
1005 break;
1006 }
1007
1008 }
1009
1010 //rest frame for reconstructed jets
1011 Double_t fPx = vsumRec.Px();
1012 Double_t fPy = vsumRec.Py();
1013 Double_t fPz = vsumRec.Pz();
1014 Double_t fE = vsumRec.E();
1015
1016 TVector3 pTrackRest[kTracks];
1017 for(Int_t j = 0; j < nTracks; j++)
1018 {
1019 vTrack[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
1020 pTrackRest[j].SetXYZ(vTrack[j].Px(), vTrack[j].Py(),vTrack[j].Pz());
1021 }
1022
1023 Double_t eRestRec[kMaxJets];
1024 Int_t idxRestRec[kMaxJets];
1025 for (Int_t j = 0; j < nAccJets; j++)
1026 {
1027 vRec[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
1028 eRestRec[j] = vRec[j].E();
1029 eSumRec += vRec[j].E();
1030 }
1031
1032 TMath::Sort(nAccJets, eRestRec, idxRestRec);
1033 for (Int_t i = 0; i < nAccJets; i++)
1034 {
1035 vRestRec[i] = vRec[idxRestRec[i]];
1036 pRestRec[i].SetXYZ(vRestRec[i].Px(), vRestRec[i].Py(), vRestRec[i].Pz());
1037 psumRestRec += pRestRec[i].Perp();
1038 }
1039
1040 if(nAccJets == 3)
1041 {
1042// if(pRest[1].DeltaPhi(pRest[2]) > 0.95 && pRest[1].DeltaPhi(pRest[2]) < 1.15)
1043// {
1044 fhInOut->Fill(nGenSel);
1045// for(Int_t j = 0; j < nTracksALL; j++)
1046// {
1047// vTracksAll[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
1048// pTracksAll[j].SetXYZ(vTracksAll[j].Px(), vTracksAll[j].Py(),vTracksAll[j].Pz());
1049// fhdPhiRec->Fill(pRest[0].DeltaPhi(pTracksAll[j]), pTracksAll[j].Perp(), fXsection);
1050// }
1051 //and the Dalitz variables and Energy distributions in the rest frame
1052 for (Int_t i = 0; i < nAccJets; i++)
1053 xRec[i] = 2*vRestRec[i].E()/eSumRec;
1054
1055 if(eSumRec <= 60)
1056 fhX3X4Rec60->Fill(xRec[0], xRec[1], fXsection);
1057
1058 if(eSumRec > 60 && eSumRec <= 100)
1059 fhX3X4Rec60100->Fill(xRec[0], xRec[1], fXsection);
1060
1061 if(eSumRec > 100)
1062 fhX3X4Rec100->Fill(xRec[0], xRec[1], fXsection);
1063
1064 if(nAccJets == 3 && nAccJets == nGenJets)
1065 {
1066 fhX3->Fill(xGen[0], TMath::Abs(xGen[0]-xRec[0])/xGen[0], fXsection);
1067 fhX4->Fill(xGen[1], TMath::Abs(xGen[1]-xRec[1])/xGen[1], fXsection);
1068 fhX5->Fill(xGen[2], TMath::Abs(xGen[2]-xRec[2])/xGen[2], fXsection);
1069 }
1070
1071 FillTopology(fhX3X4Rec, fhMu34Rec, fhMu45Rec, fhMu35Rec, xRec, pRestRec, fXsection);
1072 }
1073 Printf("%s:%d",(char*)__FILE__,__LINE__);
1074
1075 PostData(1, fList);
1076
1077 Printf("%s:%d Data Posted",(char*)__FILE__,__LINE__);
1078
1079}
1080
1081//__________________________________________________________________________________________________________________________________________________
1082void AliAnalysisTaskThreeJets::Terminate(Option_t *)
1083{
1084 printf("AnalysisJetCorrelation::Terminate()");
1085
1086}
1087
1088//_______________________________________User defined functions_____________________________________________________________________________________
1089void AliAnalysisTaskThreeJets::FillTopology(TH2F * Dalitz, TH1F * fhMu34, TH1F * fhMu45, TH1F * fhMu35, Double_t * x, TVector3 * pRest, Double_t xsection)
1090{
393afcc4 1091 //
1092 // fill the topology histos
1093 //
95392a57 1094 Dalitz->Fill(x[0], x[1], xsection);
1095 fhMu35->Fill(TMath::Sqrt(x[0]*x[2]*(1-(pRest[0].Unit()).Dot(pRest[2].Unit()))/2), xsection);
1096 fhMu34->Fill(TMath::Sqrt(x[0]*x[1]*(1-(pRest[0].Unit()).Dot(pRest[1].Unit()))/2), xsection);
1097 fhMu45->Fill(TMath::Sqrt(x[1]*x[2]*(1-(pRest[1].Unit()).Dot(pRest[2].Unit()))/2), xsection);
1098}
1099
1100//_____________________________________________________________________________________________________________________________
1101
1102Bool_t AliAnalysisTaskThreeJets::IsPrimChar(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
1103{
1104 //
1105 // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
1106 // shall be counted as a primary particle
1107 //
1108 // This function or a equivalent should be available in some common place of AliRoot
1109 //
1110 // WARNING: Call this function only for particles that are among the particles from the event generator!
1111 // --> stack->Particle(id) with id < stack->GetNprimary()
1112
1113 // if the particle has a daughter primary, we do not want to count it
1114 if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
1115 {
1116 if (adebug)
1117 printf("Dropping particle because it has a daughter among the primaries.\n");
1118 return kFALSE;
1119 }
1120
1121 Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
1122
1123
1124 // skip quarks and gluon
1125 if (pdgCode <= 10 || pdgCode == 21)
1126 {
1127 if (adebug)
1128 printf("Dropping particle because it is a quark or gluon.\n");
1129 return kFALSE;
1130 }
1131
1132 Int_t status = aParticle->GetStatusCode();
1133 // skip non final state particles..
1134 if(status!=1){
1135 if (adebug)
1136 printf("Dropping particle because it is not a final state particle.\n");
1137 return kFALSE;
1138 }
1139
1140 if (strcmp(aParticle->GetName(),"XXX") == 0)
1141 {
1142 Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode);
1143 return kFALSE;
1144 }
1145
1146 TParticlePDG* pdgPart = aParticle->GetPDG();
1147
1148 if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
1149 {
1150 Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode);
1151 return kFALSE;
1152 }
1153
1154 if (pdgPart->Charge() == 0)
1155 {
1156 if (adebug)
1157 printf("Dropping particle because it is not charged.\n");
1158 return kFALSE;
1159 }
1160
1161 return kTRUE;
1162}
1163
1164//______________________________________________________________________________________________________
1165
393afcc4 1166void AliAnalysisTaskThreeJets::GetThrustAxis(TVector3 &n01, TVector3 * pTrack, const Int_t &nTracks)
95392a57 1167{
393afcc4 1168 //
1169 // fetch the thrust axis
1170 //
95392a57 1171 TVector3 psum;
1172 Double_t psum1 = 0;
1173 Double_t psum2 = 0;
1174 Double_t thrust[kTracks];
1175 Double_t th = -3;
1176
1177 for(Int_t j = 0; j < nTracks; j++)
1178 {
1179 psum.SetXYZ(0., 0., 0.);
1180 psum1 = 0;
1181 psum2 = 0;
1182 for(Int_t i = 0; i < nTracks; i++)
1183 {
1184 psum1 += (TMath::Abs(n01.Dot(pTrack[i])));
1185 psum2 += pTrack[i].Mag();
1186
1187 if (n01.Dot(pTrack[i]) > 0) psum += pTrack[i];
1188 if (n01.Dot(pTrack[i]) < 0) psum -= pTrack[i];
1189 }
1190
1191 thrust[j] = psum1/psum2;
1192
1193 if(th == thrust[j])
1194 break;
1195
1196 th = thrust[j];
1197
1198 n01 = psum.Unit();
1199 }
1200}
1201
1202//___________________________________________________________________________________________________________
1203
1204void AliAnalysisTaskThreeJets::GetEventShapes(TVector3 &n01, TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes)
1205{
393afcc4 1206 //
1207 // get the event shape
1208 //
95392a57 1209 TVector3 psum;
1210 Double_t psum1 = 0;
1211 Double_t psum2 = 0;
1212 Double_t thrust[kTracks];
1213 Double_t th = -3;
1214
1215 //Sphericity calculation
1216 TMatrixDSym m(3);
1217 Double_t s00 = 0;
1218 Double_t s01 = 0;
1219 Double_t s02 = 0;
1220
1221 Double_t s10 = 0;
1222 Double_t s11 = 0;
1223 Double_t s12 = 0;
1224
1225 Double_t s20 = 0;
1226 Double_t s21 = 0;
1227 Double_t s22 = 0;
1228
1229 Double_t ptot = 0;
1230
1231 Double_t c = -10;
1232
1233 for(Int_t j = 0; j < nTracks; j++)
1234 {
1235 psum.SetXYZ(0., 0., 0.);
1236 psum1 = 0;
1237 psum2 = 0;
1238 for(Int_t i = 0; i < nTracks; i++)
1239 {
1240 psum1 += (TMath::Abs(n01.Dot(pTrack[i])));
1241 psum2 += pTrack[i].Mag();
1242
1243 if (n01.Dot(pTrack[i]) > 0) psum += pTrack[i];
1244 if (n01.Dot(pTrack[i]) < 0) psum -= pTrack[i];
1245 }
1246
1247 thrust[j] = psum1/psum2;
1248 if(th == thrust[j])
1249 break;
1250
1251 th = thrust[j];
1252
1253 n01 = psum.Unit();
1254 }
1255
1256 eventShapes[0] = th;
1257
1258 for(Int_t j = 0; j < nTracks; j++)
1259 {
1260 s00 = s00 + (pTrack[j].Px()*pTrack[j].Px())/pTrack[j].Mag();
1261 s01 = s01 + (pTrack[j].Px()*pTrack[j].Py())/pTrack[j].Mag();
1262 s02 = s02 + (pTrack[j].Px()*pTrack[j].Pz())/pTrack[j].Mag();
1263
1264 s10 = s10 + (pTrack[j].Py()*pTrack[j].Px())/pTrack[j].Mag();
1265 s11 = s11 + (pTrack[j].Py()*pTrack[j].Py())/pTrack[j].Mag();
1266 s12 = s12 + (pTrack[j].Py()*pTrack[j].Pz())/pTrack[j].Mag();
1267
1268 s20 = s20 + (pTrack[j].Pz()*pTrack[j].Px())/pTrack[j].Mag();
1269 s21 = s21 + (pTrack[j].Pz()*pTrack[j].Py())/pTrack[j].Mag();
1270 s22 = s22 + (pTrack[j].Pz()*pTrack[j].Pz())/pTrack[j].Mag();
1271
1272 ptot += pTrack[j].Mag();
1273 }
1274
1275 if(ptot !=0)
1276 {
1277 m(0,0) = s00/ptot;
1278 m(0,1) = s01/ptot;
1279 m(0,2) = s02/ptot;
1280
1281 m(1,0) = s10/ptot;
1282 m(1,1) = s11/ptot;
1283 m(1,2) = s12/ptot;
1284
1285 m(2,0) = s20/ptot;
1286 m(2,1) = s21/ptot;
1287 m(2,2) = s22/ptot;
1288
1289 TMatrixDSymEigen eigen(m);
1290 TVectorD eigenVal = eigen.GetEigenValues();
1291
1292 Double_t sphericity = (3/2)*(eigenVal(2)+eigenVal(1));
1293 eventShapes[1] = sphericity;
1294
1295 Double_t aplaarity = (3/2)*(eigenVal(2));
1296 eventShapes[2] = aplaarity;
1297
1298 c = 3*(eigenVal(0)*eigenVal(1)+eigenVal(0)*eigenVal(2)+eigenVal(1)*eigenVal(2));
1299 eventShapes[3] = c;
1300 }
1301}
1302
1303//__________________________________________________________________________________________________________________________
1304
393afcc4 1305Double_t AliAnalysisTaskThreeJets::Exponent(Double_t x,const Double_t * const par) const
95392a57 1306{
1307 return par[0]*TMath::Power(1/TMath::E(), TMath::Power(par[1]/x, par[2])+0.5*TMath::Power((x-par[3])/par[0], 2))+par[4]*x;
1308}
1309
393afcc4 1310Double_t AliAnalysisTaskThreeJets::Exponent2(Double_t x,const Double_t * const par) const
95392a57 1311{
1312 return par[0]*TMath::Power(1/TMath::E(), TMath::Power(par[1]/x, par[2]))+par[3]*x;
1313}
1314
393afcc4 1315Double_t AliAnalysisTaskThreeJets::Gauss(Double_t x,const Double_t * const par) const
95392a57 1316{
1317 return 1/(par[1])*TMath::Power(1/TMath::E(), 0.5*(x-par[0])*(x-par[0])/(par[1]*par[1]));
1318}
1319
393afcc4 1320Double_t AliAnalysisTaskThreeJets::Total(Double_t x,const Double_t * const par) const
95392a57 1321{
1322 return Exponent(x, par)+Gauss(x, &par[4]);
1323}
1324
1325