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