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