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