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