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