New book-keeping class
[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 "AliAODMCParticle.h"
48 //#include "AliGenPythiaEventHeader.h"
49 //#include "AliMCEvent.h"
50 //#include "AliStack.h"
51
52
53 #include "AliAnalysisHelperJetTasks.h"
54
55
56 ClassImp(AliAnalysisTaskThreeJets)
57
58   AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets() : AliAnalysisTaskSE(),
59                                                    
60                                                          fAOD(0x0),
61                                                          fUseMC(0x0),
62                                                  
63                                                          fBranchRec(""),
64                                                          fBranchGen(""),
65                                                 
66                                                          fUseAODInput(kFALSE),
67
68                                                          fR(0x0),
69                                                          fList(0x0),
70
71                                                          fGlobVar(1),
72                                                          fXsection(1),
73
74                                                          fhX3X4Rec(0x0),
75                                                          fhX3X4Gen(0x0),
76                                                          
77                                                          fhMu35Rec(0x0),
78                                                          fhMu34Rec(0x0),
79                                                          fhMu45Rec(0x0),
80     
81                                                          fhMu35Gen(0x0),
82                                                          fhMu34Gen(0x0),
83                                                          fhMu45Gen(0x0),
84         
85                                                          fhInOut(0x0),
86
87                                                          fhThrustRec2(0x0),
88                                                          fhThrustRec3(0x0),
89
90                                                          fhThrustGen2(0x0),
91                                                          fhThrustGen3(0x0),
92
93                                                          fhCGen2(0x0),
94                                                          fhCGen3(0x0),
95
96                                                          fhSGen2(0x0),
97                                                          fhSGen3(0x0),
98
99                                                          fhAGen2(0x0),
100                                                          fhAGen3(0x0),
101
102                                                          fhCRec2(0x0),
103                                                          fhCRec3(0x0),
104
105                                                          fhSRec2(0x0),
106                                                          fhSRec3(0x0),
107
108                                                          fhARec2(0x0),
109                                                          fhARec3(0x0),
110
111                                                          fhX3(0x0),
112                                                          fhX4(0x0),
113                                                          fhX5(0x0),
114
115                                                          fhXSec(0x0),
116
117                                                          fhX3X4Rec60(0x0),
118                                                          fhX3X4Rec60100(0x0),
119                                                          fhX3X4Rec100(0x0),
120
121                                                          fhX3X4Gen60(0x0),
122                                                          fhX3X4Gen60100(0x0),
123                                                          fhX3X4Gen100(0x0),
124
125                                                          fhdPhiThrustGen(0x0),
126                                                          fhdPhiThrustGenALL(0x0),
127
128                                                          fhdPhiThrustRec(0x0),
129                                                          fhdPhiThrustRecALL(0x0)
130
131 {
132
133 }
134
135 AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets(const char * name):
136   AliAnalysisTaskSE(name),
137   
138   fAOD(0x0),
139   fUseMC(0x0),
140   
141   fBranchRec(""),
142   fBranchGen(""),
143   
144   fUseAODInput(kFALSE),
145
146   fR(0x0),
147   fList(0x0),
148
149   fGlobVar(1),
150   fXsection(1),
151
152   fhX3X4Rec(0x0),
153   fhX3X4Gen(0x0),
154   
155   fhMu35Rec(0x0),
156   fhMu34Rec(0x0),
157   fhMu45Rec(0x0),
158   
159   fhMu35Gen(0x0),
160   fhMu34Gen(0x0),
161   fhMu45Gen(0x0),
162   
163   fhInOut(0x0),
164   
165   fhThrustRec2(0x0),
166   fhThrustRec3(0x0),
167   
168   fhThrustGen2(0x0),
169   fhThrustGen3(0x0),
170   
171   fhCGen2(0x0),
172   fhCGen3(0x0),
173   
174   fhSGen2(0x0),
175   fhSGen3(0x0),
176
177   fhAGen2(0x0),
178   fhAGen3(0x0),
179   
180   fhCRec2(0x0),
181   fhCRec3(0x0),
182   
183   fhSRec2(0x0),
184   fhSRec3(0x0),
185   
186   fhARec2(0x0),
187   fhARec3(0x0),
188   
189   fhX3(0x0),
190   fhX4(0x0),
191   fhX5(0x0),
192   
193   fhXSec(0x0),
194   
195   fhX3X4Rec60(0x0),
196   fhX3X4Rec60100(0x0),
197   fhX3X4Rec100(0x0),
198   
199   fhX3X4Gen60(0x0),
200   fhX3X4Gen60100(0x0),
201   fhX3X4Gen100(0x0),
202   
203   fhdPhiThrustGen(0x0),
204   fhdPhiThrustGenALL(0x0),
205   
206   fhdPhiThrustRec(0x0),
207   fhdPhiThrustRecALL(0x0)
208 {
209   DefineOutput(1, TList::Class());
210 }
211
212
213
214 Bool_t AliAnalysisTaskThreeJets::Notify()
215 {
216   //
217   // Implemented Notify() to read the cross sections
218   // and number of trials from pyxsec.root
219   // 
220
221
222   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
223   Float_t xsection = 0;
224   Float_t ftrials  = 1;
225
226   Float_t fAvgTrials = 1;
227   if(tree){
228     TFile *curfile = tree->GetCurrentFile();
229     if (!curfile) {
230       Error("Notify","No current file");
231       return kFALSE;
232     }
233     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
234     // construct a poor man average trials 
235     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
236     if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; // CKB take this into account for normalisation
237   }  
238
239   if(xsection>0)fXsection  = xsection;
240
241   return kTRUE;
242
243 }
244
245
246 //___________________________________________________________________________________________________________________________________
247 void AliAnalysisTaskThreeJets::UserCreateOutputObjects()
248 {
249   //
250   // Create the output container
251   //
252   //  Printf("Analysing event  %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
253
254   printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
255
256   fList = new TList();
257
258   fhX3X4Gen = new TH2F("X3vsX4Gen", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
259   fhX3X4Gen->SetXTitle("X_{3}");
260   fhX3X4Gen->SetYTitle("X_{4}");
261   fhX3X4Gen->Sumw2();
262   fList->Add(fhX3X4Gen);
263      
264   fhX3X4Rec = new TH2F("X3vsX4Rec", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
265   fhX3X4Rec->SetXTitle("X_{3}");
266   fhX3X4Rec->SetYTitle("X_{4}");
267   fhX3X4Rec->Sumw2();
268   fList->Add(fhX3X4Rec);
269      
270   fhMu35Rec = new TH1F("Mu35Rec", "", 20,0.1,0.8);
271   fhMu35Rec->Sumw2();
272   fhMu35Rec->SetXTitle("#mu_{35}");
273   fhMu35Rec->SetYTitle("#frac{dN}{d#mu_{35Rec}}");  
274   fList->Add(fhMu35Rec);
275           
276   fhMu34Rec = new TH1F("Mu34Rec", "", 20,0.5,1);
277   fhMu34Rec->Sumw2();
278   fhMu34Rec->SetXTitle("#mu_{34}");
279   fhMu34Rec->SetYTitle("#frac{dN}{d#mu_{34}}");
280   fList->Add(fhMu34Rec);
281   
282   fhMu45Rec = new TH1F("Mu45Rec", "", 20,0,0.65);
283   fhMu45Rec->Sumw2();
284   fhMu45Rec->SetXTitle("#mu_{45}");
285   fhMu45Rec->SetYTitle("#frac{dN}{d#mu_{45}}");
286   fList->Add(fhMu45Rec);
287      
288   fhMu35Gen = new TH1F("Mu35Gen", "", 20,0.1,0.8);
289   fhMu35Gen->Sumw2();
290   fhMu35Gen->SetXTitle("#mu_{35Gen}");
291   fhMu35Gen->SetYTitle("#frac{dN}{d#mu_{35Gen}}");  
292   fList->Add(fhMu35Gen);
293           
294   fhMu34Gen = new TH1F("Mu34Gen", "", 20,0.5,1);
295   fhMu34Gen->Sumw2();
296   fhMu34Gen->SetXTitle("#mu_{34Gen}");
297   fhMu34Gen->SetYTitle("#frac{dN}{d#mu_{34Gen}}");
298   fList->Add(fhMu34Gen);
299   
300   fhMu45Gen = new TH1F("Mu45Gen", "", 20,0,0.65);
301   fhMu45Gen->Sumw2();
302   fhMu45Gen->SetXTitle("#mu_{45Gen}");
303   fhMu45Gen->SetYTitle("#frac{dN}{d#mu_{45}}");
304   fList->Add(fhMu45Gen);
305
306   fhInOut = new TH1I("InOut", "", 6, 0, 6);
307   fhInOut->SetXTitle("#RecJets_{GenJets=3}");
308   fhInOut->SetYTitle("#Entries");
309   fList->Add(fhInOut);
310
311   fhThrustGen2 = new TH1F("ThrustGen2", "", 50, 0.5, 1);
312   fhThrustGen2->Sumw2();
313   fList->Add(fhThrustGen2);
314
315   fhThrustGen3 = new TH1F("ThrustGen3", "", 50, 0.5, 1);
316   fhThrustGen3->Sumw2();
317   fList->Add(fhThrustGen3);
318
319   fhThrustRec2 = new TH1F("ThrustRec2", "", 50, 0.5, 1);
320   fhThrustRec2->Sumw2();
321   fList->Add(fhThrustRec2);
322
323   fhThrustRec3 = new TH1F("ThrustRec3", "", 50, 0.5, 1);
324   fhThrustRec3->Sumw2();
325   fList->Add(fhThrustRec3);
326
327   fhCGen2 = new TH1F("CGen2", "", 100, 0, 1);
328   fhCGen2->Sumw2();
329   fList->Add(fhCGen2);
330
331   fhCGen3 = new TH1F("CGen3", "", 100, 0, 1);
332   fhCGen3->Sumw2();
333   fList->Add(fhCGen3);
334
335   fhCRec2 = new TH1F("CRec2", "", 100, 0, 1);
336   fhCRec2->Sumw2();
337   fList->Add(fhCRec2);
338
339   fhCRec3 = new TH1F("CRec3", "", 100, 0, 1);
340   fhCRec3->Sumw2();
341   fList->Add(fhCRec3);
342
343   fhSGen2 = new TH1F("SGen2", "", 100, 0, 1);
344   fList->Add(fhSGen2);
345
346   fhSGen3 = new TH1F("SGen3", "", 100, 0, 1);
347   fList->Add(fhSGen3);
348
349   fhSRec2 = new TH1F("SRec2", "", 100, 0, 1);
350   fList->Add(fhSRec2);
351
352   fhSRec3 = new TH1F("SRec3", "", 100, 0, 1);
353   fList->Add(fhSRec3);
354
355   fhAGen2 = new TH1F("AGen2", "", 50, 0, 0.5);
356   fList->Add(fhAGen2);
357
358   fhAGen3 = new TH1F("AGen3", "", 50, 0, 0.5);
359   fList->Add(fhAGen3);
360
361   fhARec2 = new TH1F("ARec2", "", 50, 0, 0.5);
362   fList->Add(fhARec2);
363
364   fhARec3 = new TH1F("ARec3", "", 50, 0, 0.5);
365   fList->Add(fhARec3);
366
367   fhX3 = new TH2F("X3", "", 22, 0.6, 1.02, 100, 0, 1);
368   fhX3->SetYTitle("|X_{3}^{MC} - X_{3}^{AOD}|/X_{3}^{MC}");
369   fhX3->SetXTitle("X_{3}");
370   fhX3->Sumw2();
371   fList->Add(fhX3);
372
373   fhX4 = new TH2F("X4", "",33, 0.4, 1.02, 100, 0, 1);
374   fhX4->SetYTitle("|X_{4}^{MC} - X_{4}^{AOD}|/X_{4}^{MC}");
375   fhX4->SetXTitle("X_{4}");
376   fhX4->Sumw2();
377   fList->Add(fhX4);
378
379   fhX5 = new TH2F("X5", "",100, 0., 1., 100, 0, 1);
380   fhX5->SetYTitle("|X_{5}^{MC} - X_{5}^{AOD}|/X_{5}^{MC}");
381   fhX5->SetXTitle("X_{5}");
382   fhX5->Sumw2();
383   fList->Add(fhX5);
384
385   fhXSec = new TProfile("XSec", "", 200, 0, 200, 0, 1);
386   fList->Add(fhXSec);
387
388   fhX3X4Rec60 = new TH2F("X3vsX4Rec60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
389   fhX3X4Rec60->SetXTitle("X_{3}");
390   fhX3X4Rec60->SetYTitle("X_{4}");
391   fhX3X4Rec60->Sumw2();
392   fList->Add(fhX3X4Rec60);
393
394   fhX3X4Rec60100 = new TH2F("X3vsX4Rec60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
395   fhX3X4Rec60100->SetXTitle("X_{3}");
396   fhX3X4Rec60100->SetYTitle("X_{4}");
397   fhX3X4Rec60100->Sumw2();
398   fList->Add(fhX3X4Rec60100);
399
400   fhX3X4Rec100 = new TH2F("X3vsX4Rec100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
401   fhX3X4Rec100->SetXTitle("X_{3}");
402   fhX3X4Rec100->SetYTitle("X_{4}");
403   fhX3X4Rec100->Sumw2();
404   fList->Add(fhX3X4Rec100);
405
406   fhX3X4Gen60 = new TH2F("X3vsX4Gen60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
407   fhX3X4Gen60->SetXTitle("X_{3}");
408   fhX3X4Gen60->SetYTitle("X_{4}");
409   fhX3X4Gen60->Sumw2();
410   fList->Add(fhX3X4Gen60);
411
412   fhX3X4Gen60100 = new TH2F("X3vsX4Gen60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
413   fhX3X4Gen60100->SetXTitle("X_{3}");
414   fhX3X4Gen60100->SetYTitle("X_{4}");
415   fhX3X4Gen60100->Sumw2();
416   fList->Add(fhX3X4Gen60100);
417
418   fhX3X4Gen100 = new TH2F("X3vsX4Gen100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
419   fhX3X4Gen100->SetXTitle("X_{3}");
420   fhX3X4Gen100->SetYTitle("X_{4}");
421   fhX3X4Gen100->Sumw2();
422   fList->Add(fhX3X4Gen100);
423
424   fhdPhiThrustGen = new TH2F("dPhiThrustGen", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
425   fhdPhiThrustGen->Sumw2();
426   fList->Add(fhdPhiThrustGen);
427
428   fhdPhiThrustGenALL = new TH2F("dPhiThrustGenALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0,  150);
429   fhdPhiThrustGenALL->Sumw2();
430   fList->Add(fhdPhiThrustGenALL);
431
432   fhdPhiThrustRec = new TH2F("dPhiThrustRec", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
433   fhdPhiThrustRec->Sumw2();
434   fList->Add(fhdPhiThrustRec);
435
436   fhdPhiThrustRecALL = new TH2F("dPhiThrustRecALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
437   fhdPhiThrustRecALL->Sumw2();
438   fList->Add(fhdPhiThrustRecALL);
439     
440   if(fDebug)Printf("UserCreateOutputObjects finished\n");
441 }
442
443 //__________________________________________________________________________________________________________________________________________
444 void AliAnalysisTaskThreeJets::Init()
445 {
446   printf("AliAnalysisJetCut::Init() \n");
447 }
448
449 //____________________________________________________________________________________________________________________________________________
450 void AliAnalysisTaskThreeJets::UserExec(Option_t * )
451 {
452   if (fDebug > 1) printf("AliAnlysisTaskThreeJets::Analysing event # %5d\n", (Int_t) fEntry);
453
454   if(fUseAODInput){
455     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
456     if(!fAOD){
457       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
458       return;
459     }    
460   }
461   else{
462     //  assume that the AOD is in the general output...
463     fAOD  = AODEvent();
464     if(!fAOD){
465       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
466       return;
467     }    
468   }
469
470 //   AliMCEvent* mcEvent =MCEvent();
471 //   if(!mcEvent){
472 //     Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
473 //     return;
474 //   }
475   
476   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
477   
478   //primary vertex
479   AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
480   if(!pvtx){
481     if (fDebug > 1)     Printf("%s:%d AOD Vertex found",(char*)__FILE__,__LINE__);
482     return;
483   }
484   
485 //   AliAODJet genJetsPythia[kMaxJets];
486 //   Int_t nPythiaGenJets = 0;
487
488   AliAODJet recJets[kMaxJets];
489   Int_t nRecJets = 0;
490   
491   //array of reconstructed jets from the AOD input
492   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
493   if(!aodRecJets){
494     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
495     return;
496   }
497   
498   // reconstructed jets
499   nRecJets = aodRecJets->GetEntries(); 
500   nRecJets = TMath::Min(nRecJets, kMaxJets);
501   
502   for(int ir = 0;ir < nRecJets;++ir)
503     {
504      AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
505      if(!tmp)continue;
506      recJets[ir] = *tmp;
507     }
508    
509   AliAODJet genJets[kMaxJets];
510   Int_t nGenJets = 0;
511   if(fUseMC){
512   // If we set a second branch for the input jets fetch this
513   TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
514     
515     if(!aodGenJets)
516       {
517         printf("NO MC jets Found\n");
518         return;
519       }
520     
521     //   //Generated jets
522     nGenJets = aodGenJets->GetEntries();
523     nGenJets = TMath::Min(nGenJets, kMaxJets);
524
525     for(Int_t ig =0 ; ig < nGenJets; ++ig)
526       {
527         AliAODJet * tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
528         if(!tmp)continue;
529         genJets[ig] = * tmp;
530       }
531   }
532   //   AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
533 //   if(!pythiaGenHeader){
534 //     Printf("!!!NO GEN HEADER AVALABLE!!!");
535 //     return;
536 //   }
537   
538 //   // Int_t ProcessType = pythiaGenHeader->ProcessType();
539 //   // if(ProcessType != 28) return;
540 //   nPythiaGenJets = pythiaGenHeader->NTriggerJets();
541 //   nPythiaGenJets = TMath::Min(nPythiaGenJets, kMaxJets);
542    
543 //  fXsection = 1;
544
545   Double_t eRec[kMaxJets];
546   Double_t eGen[kMaxJets];
547  
548   Double_t eJetRec[kMaxJets];
549   //  Double_t EJetGen[kMaxJets];
550
551   AliAODJet jetRec[kMaxJets];
552   AliAODJet jetGen[kMaxJets];
553
554   Int_t idxRec[kMaxJets];
555   Int_t idxGen[kMaxJets];
556    
557   Double_t xRec[kMaxJets]; 
558   Double_t xGen[kMaxJets];
559
560   Double_t eSumRec = 0;
561   Double_t eSumGen = 0;
562   
563   TLorentzVector vRec[kMaxJets];
564   TLorentzVector vRestRec[kMaxJets];
565   
566   TLorentzVector vGen[kMaxJets];
567   TLorentzVector vRestGen[kMaxJets];
568
569   TLorentzVector vsumRec;
570   TLorentzVector vsumGen;
571
572   TVector3 pRec[kMaxJets];
573   TVector3 pGen[kMaxJets];
574
575   TVector3 pTrack[kTracks];
576
577   TVector3 pRestRec[kMaxJets];
578   TVector3 pRestGen[kMaxJets];
579
580   Double_t psumRestRec = 0;
581   //  Double_t psumRestGen = 0;
582 //   //Pythia_________________________________________________________________________________________________________________
583
584 //   for(int ip = 0;ip < nPythiaGenJets;++ip)
585 //     {
586 //       if(ip>=kMaxJets)continue;
587 //       Float_t p[4];
588 //       pythiaGenHeader->TriggerJet(ip,p);
589 //       genJetsPythia[ip].SetPxPyPzE(p[0],p[1],p[2],p[3]);
590 //     }
591   
592 //_________________________________________________________________________________________________________________________
593  
594  
595 //________histos for MC___________________________________________________________________________________________________________
596   Int_t nGenSel = 0;
597   if(fUseMC){
598     Int_t counter = 0;
599     Int_t tag = 0;
600     
601     AliAODJet selJets[kMaxJets];
602     
603     for(Int_t i = 0; i < nGenJets; i++)
604       {
605         if(nGenJets == 1)
606           {
607             selJets[nGenSel] = genJets[i];
608             nGenSel++;
609           }
610         else
611           {
612             counter = 0;
613             tag = 0;
614             for(Int_t j = 0; j < nGenJets; j++)
615               {
616                 if(i!=j)
617                   {
618                     Double_t dRij = genJets[i].DeltaR(&genJets[j]);
619                     counter++;
620                     if(dRij > 2*fR) tag++;
621                   }
622               }
623             if(counter!=0)
624               {
625                 if(tag/counter == 1)
626                   {
627                     selJets[nGenSel] = genJets[i];
628                     nGenSel++;
629                   }
630               }
631           }
632       }
633
634     if(nGenSel == 0) return;
635     
636     for (Int_t gj = 0; gj < nGenSel; gj++)
637       {
638         eGen[gj] = selJets[gj].E();
639       }
640     
641     TMath::Sort(nGenSel, eGen, idxGen);
642     for (Int_t ig = 0; ig < nGenSel; ig++)
643       {
644         jetGen[ig] = selJets[idxGen[ig]];
645       }
646     
647     fhXSec->Fill(jetGen[0].Pt(), fXsection);
648     //  AliStack * stack = mcEvent->Stack();
649     
650     Int_t nMCtracks = 0;
651     Double_t eTracksMC[kTracks];
652     Double_t pTracksMC[kTracks];
653     Int_t idxTracksMC[kTracks];
654     TLorentzVector jetTracksMC[kTracks];
655     TLorentzVector jetTracksSortMC[kTracks];
656     TVector3 pTrackMC[kTracks];
657     TLorentzVector vTrackMCAll[kTracks];
658     Double_t pTrackMCAll[kTracks];
659     TLorentzVector vTrackMC[kTracks];
660     TVector3 pTrackMCBoost[kTracks];
661     Double_t eventShapes[4];
662     
663     Int_t nAccTr = 0;
664     Int_t nInJet[kMaxJets];
665     TLorentzVector inJetPartV[kMaxJets][kTracks];
666     Int_t nAllTracksMC = 0;
667     TVector3 n01MC;
668
669     TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
670     if(!tca){
671       if(fDebug)Printf("NO Ref Tracks\n");
672       tca = 0;
673     }
674     else{
675       nMCtracks = tca->GetEntries();
676       for(Int_t iTrack = 0; iTrack < nMCtracks; iTrack++)
677         {
678           //      TParticle * part = (TParticle*)stack->Particle(iTrack);
679           AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(iTrack));
680           if (!part) continue;
681           if(!part->IsPhysicalPrimary())continue;      
682           Double_t fEta = part->Eta();
683           if(TMath::Abs(fEta) > .9) continue;
684           vTrackMCAll[nAllTracksMC].SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
685           pTrackMCAll[nAllTracksMC] = part->Pt();
686           nAllTracksMC++;
687         }
688       if(nAllTracksMC == 0) return;
689       for(Int_t iJet = 0; iJet < nGenSel; iJet++)
690         {
691           Int_t nJetTracks = 0;
692           for(Int_t i = 0; i < nAllTracksMC; i++)
693             {
694               Double_t dPhi = (jetGen[iJet].Phi()-vTrackMCAll[i].Phi());
695               if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi();
696               if(dPhi < (-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi();
697               Double_t dEta = (jetGen[iJet].Eta()-vTrackMCAll[i].Eta());
698               Double_t deltaR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
699               if(deltaR < fR && vTrackMCAll[i].Pt() > 1.5)
700                 {
701                   jetTracksMC[nAccTr] = vTrackMCAll[i];
702                   eTracksMC[nAccTr] = vTrackMCAll[i].E();
703                   pTracksMC[nAccTr] = vTrackMCAll[i].Pt();
704                   inJetPartV[iJet][nJetTracks].SetPxPyPzE(vTrackMCAll[i].Px(), vTrackMCAll[i].Py(), vTrackMCAll[i].Pz(),vTrackMCAll[i].E());
705                   nAccTr++;
706                   nJetTracks++;
707                 }
708             }
709           nInJet[iJet] = nJetTracks;
710         }
711       
712       if(nAccTr == 0) return;
713       if(fDebug)Printf("*********** Number of Jets : %d ***************\n", nGenSel);  
714       Double_t pTav[kMaxJets];
715       for(Int_t i = 0; i < nGenSel; i++)
716         {
717           Double_t pTsum = 0;
718           if(fDebug)Printf("*********** Number of particles in Jet %d = %d *******************\n", i+3, nInJet[i]);
719           for(Int_t iT = 0; iT < nInJet[i]; iT++)
720             {
721               Double_t pt = inJetPartV[i][iT].Pt();
722               pTsum += pt;
723             }
724           pTav[i] = pTsum/nInJet[i]; 
725         }
726       
727       TMath::Sort(nAllTracksMC, pTrackMCAll, idxTracksMC);
728       for(Int_t i = 0; i < nAllTracksMC; i++)
729         {
730           jetTracksSortMC[i] = vTrackMCAll[idxTracksMC[i]];
731           pTrackMC[i].SetXYZ(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz());
732           vTrackMC[i].SetPxPyPzE(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz(), jetTracksSortMC[i].E());
733         }
734       
735       n01MC = pTrackMC[0].Unit();
736       n01MC.SetZ(0.);
737       
738       //Thrust calculation, iterative method
739       if(nGenSel > 1)
740         { 
741           //       if(fGlobVar == 1)
742           //    {
743           if(fDebug)Printf("**************Shapes for MC*************");
744           AliAnalysisHelperJetTasks::GetEventShapes(n01MC, pTrackMC, nAllTracksMC, eventShapes);
745           //    }
746           if(eventShapes[0] < 2/TMath::Pi()){
747             Double_t eventShapesTest[4];
748             TVector3 n01Test;
749             Int_t rnd_max = nAllTracksMC;
750             Int_t k = (rand()%rnd_max)+3;
751             while(TMath::Abs(pTrackMC[k].X()) < 10e-5 && TMath::Abs(pTrackMC[k].Y()) < 10e-5){
752               k--;
753             }
754             n01Test = pTrackMC[k].Unit();
755             n01Test.SetZ(0.);
756             AliAnalysisHelperJetTasks::GetEventShapes(n01Test, pTrackMC, nAllTracksMC, eventShapesTest);
757             eventShapes[0] = TMath::Max(eventShapes[0], eventShapesTest[0]);
758             if(TMath::Abs(eventShapes[0]-eventShapesTest[0]) < 10e-7) n01MC = n01Test;
759           }
760         
761           Double_t s = eventShapes[1];
762           Double_t a = eventShapes[2];
763           Double_t c = eventShapes[3];
764           
765           switch(nGenSel)
766             {
767             case 2:
768               {
769                 fhAGen2->Fill(a);
770                 fhSGen2->Fill(s);
771                 fhCGen2->Fill(c);
772               }
773               break;
774             case 3:
775               {
776                 fhAGen3->Fill(a);
777                 fhSGen3->Fill(s);
778                 fhCGen3->Fill(c);
779               }
780               break;
781             }
782           Double_t thrust01MC = eventShapes[0];
783           
784           switch(nGenSel)
785             {
786             case 2:
787               fhThrustGen2->Fill(thrust01MC, fXsection);
788               break;
789             case 3:
790               fhThrustGen3->Fill(thrust01MC, fXsection);
791               break;
792             }
793         }
794     }
795   
796     
797     //rest frame MC jets
798     for (Int_t i = 0; i < nGenSel; ++i)
799       {
800         vGen[i].SetPxPyPzE(jetGen[i].Px(), jetGen[i].Py(), jetGen[i].Pz(), jetGen[i].E());
801         pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
802         vsumGen += vGen[i];
803       }
804     if(tca){
805       if(eventShapes[0] > 0.8 && nGenSel > 1)
806         {
807           for(Int_t i = 0; i < nGenSel; i++)
808             fhdPhiThrustGen->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
809         }
810  
811       if(eventShapes[0] <= 0.8 && nGenSel > 1)   
812         {
813           for(Int_t i = 0; i < nGenSel; i++)
814             fhdPhiThrustGenALL->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
815         }
816     }
817     
818     Double_t fPxGen = vsumGen.Px();
819     Double_t fPyGen = vsumGen.Py();
820     Double_t fPzGen = vsumGen.Pz();
821     Double_t fEGen = vsumGen.E();
822     
823     Double_t eRestGen[kMaxJets];
824     for (Int_t j = 0; j < nGenSel; j++)
825       {
826         vGen[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
827         eRestGen[j] = vGen[j].E();
828       }
829     
830     for (Int_t j = 0; j < nAccTr; j++)
831       {
832         vTrackMC[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
833         pTrackMCBoost[j].SetXYZ(vTrackMC[j].Px(),vTrackMC[j].Py(),vTrackMC[j].Pz());
834       }
835     
836     Int_t idxRestGen[kMaxJets];
837     TMath::Sort(nGenSel, eRestGen, idxRestGen);
838     for(Int_t j = 0; j < nGenSel; j++)
839       {
840         vRestGen[j] = vGen[idxRestGen[j]];
841         eSumGen += vRestGen[j].E();
842       }
843
844     if (nGenSel == 3)
845       {
846         //      if(nInJet[0] < 3 || nInJet[1] < 3 || nInJet[2] < 3) return;
847         //       if(pRestGen[1].DeltaPhi(pRestGen[2]) > 0.95 && pRestGen[1].DeltaPhi(pRestGen[2]) < 1.15)
848         //      {
849       
850         for(Int_t i = 0; i < nGenSel; i++)
851           {
852             xGen[i] = 2*vRestGen[i].E()/eSumGen;
853         }
854       
855         if(fDebug)Printf("***************** Values of Dalitz variables are : %f, %f, %f ****************\n", xGen[0], xGen[1], xGen[2]);
856         
857         if(fDebug)Printf("***************** fXSection = %f ******************\n", fXsection);
858         if(eSumGen <= 60)
859           fhX3X4Gen60->Fill(xGen[0], xGen[1], fXsection);
860         
861         if(eSumGen > 60 && eSumGen <= 100)
862           fhX3X4Gen60100->Fill(xGen[0], xGen[1], fXsection);
863         
864         if(eSumGen > 100)
865           fhX3X4Gen100->Fill(xGen[0], xGen[1], fXsection);
866         
867         FillTopology(fhX3X4Gen, fhMu34Gen, fhMu45Gen, fhMu35Gen, xGen, pRestGen, fXsection);
868       }
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           if(fUseMC) 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