git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
account for limits of thrust variable (Sona)
[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
20d01a2f 24#include <cstdlib>
25#include <iostream>
95392a57 26#include <TFile.h>
27#include <TH1F.h>
95392a57 28#include <TH2F.h>
95392a57 29#include <TList.h>
30#include <TTree.h>
95392a57 31#include <TLorentzVector.h>
32#include <TClonesArray.h>
95392a57 33#include <TRefArray.h>
95392a57 34#include <TVector3.h>
35#include <TVectorD.h>
36#include <TMatrixDSym.h>
37#include <TMatrixDSymEigen.h>
95392a57 38#include <TProfile.h>
95392a57 39
40#include "AliAnalysisTaskThreeJets.h"
41#include "AliAnalysisManager.h"
95392a57 42#include "AliAODEvent.h"
43#include "AliAODVertex.h"
44#include "AliAODHandler.h"
45#include "AliAODTrack.h"
46#include "AliAODJet.h"
393afcc4 47#include "AliGenPythiaEventHeader.h"
95392a57 48#include "AliMCEvent.h"
49#include "AliStack.h"
95392a57 50
95392a57 51
52#include "AliAnalysisHelperJetTasks.h"
53
54
55ClassImp(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)
20d01a2f 128
95392a57 129{
130
131}
132
133AliAnalysisTaskThreeJets::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
211Bool_t AliAnalysisTaskThreeJets::Notify()
212{
213 //
214 // Implemented Notify() to read the cross sections
215 // and number of trials from pyxsec.root
216 //
217 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
218 UInt_t ntrials = 0;
219 if(tree){
220 TFile *curfile = tree->GetCurrentFile();
221 if (!curfile) {
222 Error("Notify","No current file");
223 return kFALSE;
224 }
225
226 TString fileName(curfile->GetName());
227 if(fileName.Contains("AliESDs.root")){
228 fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
229 }
230 else if(fileName.Contains("AliAOD.root")){
231 fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
232 }
233 else if(fileName.Contains("galice.root")){
234 // for running with galice and kinematics alone...
235 fileName.ReplaceAll("galice.root", "pyxsec.root");
236 }
237 TFile *fxsec = TFile::Open(fileName.Data());
238 if(!fxsec){
239 Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
240 // no a severe condition
241 return kTRUE;
242 }
243 TTree *xtree = (TTree*)fxsec->Get("Xsection");
244 if(!xtree){
245 Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
246 }
247 xtree->SetBranchAddress("xsection",&fXsection);
248 xtree->SetBranchAddress("ntrials",&ntrials);
249 xtree->GetEntry(0);
250 }
251
252 return kTRUE;
253}
254
255
256//___________________________________________________________________________________________________________________________________
257void AliAnalysisTaskThreeJets::UserCreateOutputObjects()
258{
259 //
260 // Create the output container
261 //
47dcf196 262 // Printf("Analysing event %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
95392a57 263
264 if(fUseAODInput){
265 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
266 if(!fAOD){
267 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
268 return;
269 }
270 }
271 else{
272 // assume that the AOD is in the general output...
273 fAOD = AODEvent();
274 if(!fAOD){
275 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
276 return;
277 }
278 }
279
280 printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
281
282 fList = new TList();
283
284 fhX3X4Gen = new TH2F("X3vsX4Gen", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
285 fhX3X4Gen->SetXTitle("X_{3}");
286 fhX3X4Gen->SetYTitle("X_{4}");
287 fhX3X4Gen->Sumw2();
288 fList->Add(fhX3X4Gen);
289
290 fhX3X4Rec = new TH2F("X3vsX4Rec", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
291 fhX3X4Rec->SetXTitle("X_{3}");
292 fhX3X4Rec->SetYTitle("X_{4}");
293 fhX3X4Rec->Sumw2();
294 fList->Add(fhX3X4Rec);
295
296 fhMu35Rec = new TH1F("Mu35Rec", "", 20,0.1,0.8);
297 fhMu35Rec->Sumw2();
298 fhMu35Rec->SetXTitle("#mu_{35}");
299 fhMu35Rec->SetYTitle("#frac{dN}{d#mu_{35Rec}}");
300 fList->Add(fhMu35Rec);
301
302 fhMu34Rec = new TH1F("Mu34Rec", "", 20,0.5,1);
303 fhMu34Rec->Sumw2();
304 fhMu34Rec->SetXTitle("#mu_{34}");
305 fhMu34Rec->SetYTitle("#frac{dN}{d#mu_{34}}");
306 fList->Add(fhMu34Rec);
307
308 fhMu45Rec = new TH1F("Mu45Rec", "", 20,0,0.65);
309 fhMu45Rec->Sumw2();
310 fhMu45Rec->SetXTitle("#mu_{45}");
311 fhMu45Rec->SetYTitle("#frac{dN}{d#mu_{45}}");
312 fList->Add(fhMu45Rec);
313
314 fhMu35Gen = new TH1F("Mu35Gen", "", 20,0.1,0.8);
315 fhMu35Gen->Sumw2();
316 fhMu35Gen->SetXTitle("#mu_{35Gen}");
317 fhMu35Gen->SetYTitle("#frac{dN}{d#mu_{35Gen}}");
318 fList->Add(fhMu35Gen);
319
320 fhMu34Gen = new TH1F("Mu34Gen", "", 20,0.5,1);
321 fhMu34Gen->Sumw2();
322 fhMu34Gen->SetXTitle("#mu_{34Gen}");
323 fhMu34Gen->SetYTitle("#frac{dN}{d#mu_{34Gen}}");
324 fList->Add(fhMu34Gen);
325
326 fhMu45Gen = new TH1F("Mu45Gen", "", 20,0,0.65);
327 fhMu45Gen->Sumw2();
328 fhMu45Gen->SetXTitle("#mu_{45Gen}");
329 fhMu45Gen->SetYTitle("#frac{dN}{d#mu_{45}}");
330 fList->Add(fhMu45Gen);
331
332 fhInOut = new TH1I("InOut", "", 6, 0, 6);
333 fhInOut->SetXTitle("#RecJets_{GenJets=3}");
334 fhInOut->SetYTitle("#Entries");
335 fList->Add(fhInOut);
336
337 fhThrustGen2 = new TH1F("ThrustGen2", "", 50, 0.5, 1);
338 fhThrustGen2->Sumw2();
339 fList->Add(fhThrustGen2);
340
341 fhThrustGen3 = new TH1F("ThrustGen3", "", 50, 0.5, 1);
342 fhThrustGen3->Sumw2();
343 fList->Add(fhThrustGen3);
344
345 fhThrustRec2 = new TH1F("ThrustRec2", "", 50, 0.5, 1);
346 fhThrustRec2->Sumw2();
347 fList->Add(fhThrustRec2);
348
349 fhThrustRec3 = new TH1F("ThrustRec3", "", 50, 0.5, 1);
350 fhThrustRec3->Sumw2();
351 fList->Add(fhThrustRec3);
352
353 fhCGen2 = new TH1F("CGen2", "", 100, 0, 1);
354 fhCGen2->Sumw2();
355 fList->Add(fhCGen2);
356
357 fhCGen3 = new TH1F("CGen3", "", 100, 0, 1);
358 fhCGen3->Sumw2();
359 fList->Add(fhCGen3);
360
361 fhCRec2 = new TH1F("CRec2", "", 100, 0, 1);
362 fhCRec2->Sumw2();
363 fList->Add(fhCRec2);
364
365 fhCRec3 = new TH1F("CRec3", "", 100, 0, 1);
366 fhCRec3->Sumw2();
367 fList->Add(fhCRec3);
368
369 fhSGen2 = new TH1F("SGen2", "", 100, 0, 1);
370 fList->Add(fhSGen2);
371
372 fhSGen3 = new TH1F("SGen3", "", 100, 0, 1);
373 fList->Add(fhSGen3);
374
375 fhSRec2 = new TH1F("SRec2", "", 100, 0, 1);
376 fList->Add(fhSRec2);
377
378 fhSRec3 = new TH1F("SRec3", "", 100, 0, 1);
379 fList->Add(fhSRec3);
380
381 fhAGen2 = new TH1F("AGen2", "", 50, 0, 0.5);
382 fList->Add(fhAGen2);
383
384 fhAGen3 = new TH1F("AGen3", "", 50, 0, 0.5);
385 fList->Add(fhAGen3);
386
387 fhARec2 = new TH1F("ARec2", "", 50, 0, 0.5);
388 fList->Add(fhARec2);
389
390 fhARec3 = new TH1F("ARec3", "", 50, 0, 0.5);
391 fList->Add(fhARec3);
392
393 fhX3 = new TH2F("X3", "", 22, 0.6, 1.02, 100, 0, 1);
394 fhX3->SetYTitle("|X_{3}^{MC} - X_{3}^{AOD}|/X_{3}^{MC}");
395 fhX3->SetXTitle("X_{3}");
396 fhX3->Sumw2();
397 fList->Add(fhX3);
398
399 fhX4 = new TH2F("X4", "",33, 0.4, 1.02, 100, 0, 1);
400 fhX4->SetYTitle("|X_{4}^{MC} - X_{4}^{AOD}|/X_{4}^{MC}");
401 fhX4->SetXTitle("X_{4}");
402 fhX4->Sumw2();
403 fList->Add(fhX4);
404
405 fhX5 = new TH2F("X5", "",100, 0., 1., 100, 0, 1);
406 fhX5->SetYTitle("|X_{5}^{MC} - X_{5}^{AOD}|/X_{5}^{MC}");
407 fhX5->SetXTitle("X_{5}");
408 fhX5->Sumw2();
409 fList->Add(fhX5);
410
411 fhXSec = new TProfile("XSec", "", 200, 0, 200, 0, 1);
412 fList->Add(fhXSec);
413
414 fhX3X4Rec60 = new TH2F("X3vsX4Rec60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
415 fhX3X4Rec60->SetXTitle("X_{3}");
416 fhX3X4Rec60->SetYTitle("X_{4}");
417 fhX3X4Rec60->Sumw2();
418 fList->Add(fhX3X4Rec60);
419
420 fhX3X4Rec60100 = new TH2F("X3vsX4Rec60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
421 fhX3X4Rec60100->SetXTitle("X_{3}");
422 fhX3X4Rec60100->SetYTitle("X_{4}");
423 fhX3X4Rec60100->Sumw2();
424 fList->Add(fhX3X4Rec60100);
425
426 fhX3X4Rec100 = new TH2F("X3vsX4Rec100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
427 fhX3X4Rec100->SetXTitle("X_{3}");
428 fhX3X4Rec100->SetYTitle("X_{4}");
429 fhX3X4Rec100->Sumw2();
430 fList->Add(fhX3X4Rec100);
431
432 fhX3X4Gen60 = new TH2F("X3vsX4Gen60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
433 fhX3X4Gen60->SetXTitle("X_{3}");
434 fhX3X4Gen60->SetYTitle("X_{4}");
435 fhX3X4Gen60->Sumw2();
436 fList->Add(fhX3X4Gen60);
437
438 fhX3X4Gen60100 = new TH2F("X3vsX4Gen60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
439 fhX3X4Gen60100->SetXTitle("X_{3}");
440 fhX3X4Gen60100->SetYTitle("X_{4}");
441 fhX3X4Gen60100->Sumw2();
442 fList->Add(fhX3X4Gen60100);
443
444 fhX3X4Gen100 = new TH2F("X3vsX4Gen100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
445 fhX3X4Gen100->SetXTitle("X_{3}");
446 fhX3X4Gen100->SetYTitle("X_{4}");
447 fhX3X4Gen100->Sumw2();
448 fList->Add(fhX3X4Gen100);
449
450 fhdPhiThrustGen = new TH2F("dPhiThrustGen", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
451 fhdPhiThrustGen->Sumw2();
452 fList->Add(fhdPhiThrustGen);
453
454 fhdPhiThrustGenALL = new TH2F("dPhiThrustGenALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
455 fhdPhiThrustGenALL->Sumw2();
456 fList->Add(fhdPhiThrustGenALL);
457
458 fhdPhiThrustRec = new TH2F("dPhiThrustRec", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
459 fhdPhiThrustRec->Sumw2();
460 fList->Add(fhdPhiThrustRec);
461
462 fhdPhiThrustRecALL = new TH2F("dPhiThrustRecALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
463 fhdPhiThrustRecALL->Sumw2();
464 fList->Add(fhdPhiThrustRecALL);
465
466 Printf("UserCreateOutputObjects finished\n");
467}
468
469//__________________________________________________________________________________________________________________________________________
470void AliAnalysisTaskThreeJets::Init()
471{
472 printf("AliAnalysisJetCut::Init() \n");
473}
474
475//____________________________________________________________________________________________________________________________________________
476void AliAnalysisTaskThreeJets::UserExec(Option_t * )
477{
478 // if (fDebug > 1) printf("Analysing event # %5d\n", (Int_t) fEntry);
479
480
481 //create an AOD handler
482 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
483
484 if(!aodH)
485 {
486 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
487 return;
488 }
489
490 AliMCEvent* mcEvent =MCEvent();
491 if(!mcEvent){
492 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
493 return;
494 }
495
496 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
497
498 //primary vertex
499 AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
500 if(!pvtx) return;
501
502 AliAODJet genJetsPythia[kMaxJets];
503 Int_t nPythiaGenJets = 0;
504
505 AliAODJet genJets[kMaxJets];
506 Int_t nGenJets = 0;
507
508 AliAODJet recJets[kMaxJets];
509 Int_t nRecJets = 0;
510
511 //array of reconstructed jets from the AOD input
512 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
513 if(!aodRecJets){
514 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
515 return;
516 }
517
518 // reconstructed jets
519 nRecJets = aodRecJets->GetEntries();
520 nRecJets = TMath::Min(nRecJets, kMaxJets);
521
522 for(int ir = 0;ir < nRecJets;++ir)
523 {
524 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
525 if(!tmp)continue;
526 recJets[ir] = *tmp;
527 }
528
529 // If we set a second branch for the input jets fetch this
530 TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
531
532 if(!aodGenJets)
533 {
534 printf("NO MC jets Found\n");
535 return;
536 }
537
538 // //Generated jets
539 nGenJets = aodGenJets->GetEntries();
540 nGenJets = TMath::Min(nGenJets, kMaxJets);
541
542 for(Int_t ig =0 ; ig < nGenJets; ++ig)
543 {
544 AliAODJet * tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
545 if(!tmp)continue;
546 genJets[ig] = * tmp;
547 }
548
549 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
550 if(!pythiaGenHeader){
551 Printf("!!!NO GEN HEADER AVALABLE!!!");
552 return;
553 }
554
555 // Int_t ProcessType = pythiaGenHeader->ProcessType();
556 // if(ProcessType != 28) return;
557 fhXSec->Fill(pythiaGenHeader->GetPtHard(), fXsection);
558 nPythiaGenJets = pythiaGenHeader->NTriggerJets();
559 nPythiaGenJets = TMath::Min(nPythiaGenJets, kMaxJets);
560
561 Double_t eRec[kMaxJets];
562 Double_t eGen[kMaxJets];
563
564 Double_t eJetRec[kMaxJets];
565 // Double_t EJetGen[kMaxJets];
566
567 AliAODJet jetRec[kMaxJets];
568 AliAODJet jetGen[kMaxJets];
569
570 Int_t idxRec[kMaxJets];
571 Int_t idxGen[kMaxJets];
572
573 Double_t xRec[kMaxJets];
574 Double_t xGen[kMaxJets];
575
576 Double_t eSumRec = 0;
577 Double_t eSumGen = 0;
578
579 TLorentzVector vRec[kMaxJets];
580 TLorentzVector vRestRec[kMaxJets];
581
582 TLorentzVector vGen[kMaxJets];
583 TLorentzVector vRestGen[kMaxJets];
584
585 TLorentzVector vsumRec;
586 TLorentzVector vsumGen;
587
588 TVector3 pRec[kMaxJets];
589 TVector3 pGen[kMaxJets];
590
591 TVector3 pTrack[kTracks];
592
593 TVector3 pRestRec[kMaxJets];
594 TVector3 pRestGen[kMaxJets];
595
596 Double_t psumRestRec = 0;
597 // Double_t psumRestGen = 0;
598 //Pythia_________________________________________________________________________________________________________________
599
600 for(int ip = 0;ip < nPythiaGenJets;++ip)
601 {
602 if(ip>=kMaxJets)continue;
603 Float_t p[4];
604 pythiaGenHeader->TriggerJet(ip,p);
605 genJetsPythia[ip].SetPxPyPzE(p[0],p[1],p[2],p[3]);
606 }
607
608//_________________________________________________________________________________________________________________________
609
610
611//________histos for MC___________________________________________________________________________________________________________
612
613 Int_t nGenSel = 0;
614 Int_t counter = 0;
615 Int_t tag = 0;
616
617 AliAODJet selJets[kMaxJets];
618
619 for(Int_t i = 0; i < nGenJets; i++)
620 {
621 if(nGenJets == 1)
622 {
623 selJets[nGenSel] = genJets[i];
624 nGenSel++;
625 }
626 else
627 {
628 counter = 0;
629 tag = 0;
630 for(Int_t j = 0; j < nGenJets; j++)
631 {
632 if(i!=j)
633 {
634 Double_t dRij = genJets[i].DeltaR(&genJets[j]);
635 counter++;
636 if(dRij > 2*fR) tag++;
637 }
638 }
639 if(counter!=0)
640 {
641 if(tag/counter == 1)
642 {
643 selJets[nGenSel] = genJets[i];
644 nGenSel++;
645 }
646 }
647 }
648 }
649
650 if(nGenSel == 0) return;
651
652 for (Int_t gj = 0; gj < nGenSel; gj++)
653 {
654 eGen[gj] = selJets[gj].E();
655 }
656
657 TMath::Sort(nGenSel, eGen, idxGen);
658 for (Int_t ig = 0; ig < nGenSel; ig++)
659 {
660 jetGen[ig] = selJets[idxGen[ig]];
661 }
662 AliStack * stack = mcEvent->Stack();
663 Int_t nMCtracks = stack->GetNprimary();
664 Double_t * eTracksMC = new Double_t[kTracks];
665 Double_t pTracksMC[kTracks];
666 Int_t * idxTracksMC = new Int_t[kTracks];
667 TLorentzVector jetTracksMC[kTracks];
668 TLorentzVector jetTracksSortMC[kTracks];
669 TVector3 pTrackMC[kTracks];
670 TLorentzVector vTrackMCAll[kTracks];
671 Double_t pTrackMCAll[kTracks];
672 TLorentzVector vTrackMC[kTracks];
673 TVector3 pTrackMCBoost[kTracks];
674 Double_t eventShapes[4];
675
676 Int_t nAccTr = 0;
677 Int_t nInJet[kMaxJets];
678 TLorentzVector inJetPartV[kMaxJets][kTracks];
679 Int_t nAllTracksMC = 0;
680
681 for(Int_t iTrack = 0; iTrack < nMCtracks; iTrack++)
682 {
683 TParticle * part = (TParticle*)stack->Particle(iTrack);
684 if (!part) continue;
685 Double_t fEta = part->Eta();
686 if(TMath::Abs(fEta) > .9) continue;
687 if(!IsPrimChar(part, nMCtracks, 0)) continue;
688 vTrackMCAll[nAllTracksMC].SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
689 pTrackMCAll[nAllTracksMC] = part->Pt();
690 nAllTracksMC++;
691 }
692 if(nAllTracksMC == 0) return;
693 for(Int_t iJet = 0; iJet < nGenSel; iJet++)
694 {
695 Int_t nJetTracks = 0;
696 for(Int_t i = 0; i < nAllTracksMC; i++)
697 {
698 Double_t dPhi = (jetGen[iJet].Phi()-vTrackMCAll[i].Phi());
699 if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi();
700 if(dPhi < (-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi();
701 Double_t dEta = (jetGen[iJet].Eta()-vTrackMCAll[i].Eta());
702 Double_t deltaR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
703 if(deltaR < fR && vTrackMCAll[i].Pt() > 1.5)
704 {
705 jetTracksMC[nAccTr] = vTrackMCAll[i];
706 eTracksMC[nAccTr] = vTrackMCAll[i].E();
707 pTracksMC[nAccTr] = vTrackMCAll[i].Pt();
708 inJetPartV[iJet][nJetTracks].SetPxPyPzE(vTrackMCAll[i].Px(), vTrackMCAll[i].Py(), vTrackMCAll[i].Pz(),vTrackMCAll[i].E());
709 nAccTr++;
710 nJetTracks++;
711 }
712 }
713 nInJet[iJet] = nJetTracks;
714 }
715
716 if(nAccTr == 0) return;
717 Printf("*********** Number of Jets : %d ***************\n", nGenSel);
718 Double_t pTav[kMaxJets];
719 for(Int_t i = 0; i < nGenSel; i++)
720 {
721 Double_t pTsum = 0;
722 Printf("*********** Number of particles in Jet %d = %d *******************\n", i+3, nInJet[i]);
723 for(Int_t iT = 0; iT < nInJet[i]; iT++)
724 {
725 Double_t pt = inJetPartV[i][iT].Pt();
726 pTsum += pt;
727 }
728 pTav[i] = pTsum/nInJet[i];
729 }
730
20d01a2f 731 TMath::Sort(nAllTracksMC, pTrackMCAll, idxTracksMC);
732 for(Int_t i = 0; i < nAllTracksMC; i++)
95392a57 733 {
20d01a2f 734 jetTracksSortMC[i] = vTrackMCAll[idxTracksMC[i]];
95392a57 735 pTrackMC[i].SetXYZ(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz());
736 vTrackMC[i].SetPxPyPzE(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz(), jetTracksSortMC[i].E());
737 }
738
739 TVector3 n01MC = pTrackMC[0].Unit();
20d01a2f 740 n01MC.SetZ(0.);
741
95392a57 742 //Thrust calculation, iterative method
743 if(nGenSel > 1)
744 {
745// if(fGlobVar == 1)
746// {
20d01a2f 747 Printf("**************Shapes for MC*************");
748 AliAnalysisHelperJetTasks::GetEventShapes(n01MC, pTrackMC, nAllTracksMC, eventShapes);
95392a57 749// }
20d01a2f 750 if(eventShapes[0] < 2/TMath::Pi()){
751 Double_t eventShapesTest[4];
752 TVector3 n01Test;
753 Int_t rnd_max = nAllTracksMC;
754 Int_t k = (rand()%rnd_max)+3;
755 while(TMath::Abs(pTrackMC[k].X()) < 10e-5 && TMath::Abs(pTrackMC[k].Y()) < 10e-5){
756 k--;
757 }
758 n01Test = pTrackMC[k].Unit();
759 n01Test.SetZ(0.);
760 AliAnalysisHelperJetTasks::GetEventShapes(n01Test, pTrackMC, nAllTracksMC, eventShapesTest);
761 eventShapes[0] = TMath::Max(eventShapes[0], eventShapesTest[0]);
762 if(TMath::Abs(eventShapes[0]-eventShapesTest[0]) < 10e-7) n01MC = n01Test;
763 }
764
95392a57 765 Double_t s = eventShapes[1];
766 Double_t a = eventShapes[2];
767 Double_t c = eventShapes[3];
768
769 switch(nGenSel)
770 {
771 case 2:
772 {
773 fhAGen2->Fill(a);
774 fhSGen2->Fill(s);
775 fhCGen2->Fill(c);
776 }
777 break;
778 case 3:
779 {
780 fhAGen3->Fill(a);
781 fhSGen3->Fill(s);
782 fhCGen3->Fill(c);
783 }
784 break;
785 }
786 Double_t thrust01MC = eventShapes[0];
787
788 switch(nGenSel)
789 {
790 case 2:
791 fhThrustGen2->Fill(thrust01MC, fXsection);
792 break;
793 case 3:
794 fhThrustGen3->Fill(thrust01MC, fXsection);
795 break;
796 }
797 }
798
799
800 //rest frame MC jets
801 for (Int_t i = 0; i < nGenSel; ++i)
802 {
803 vGen[i].SetPxPyPzE(jetGen[i].Px(), jetGen[i].Py(), jetGen[i].Pz(), jetGen[i].E());
804 pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
805 vsumGen += vGen[i];
806 }
20d01a2f 807 if(eventShapes[0] > 0.8 && nGenSel > 1)
95392a57 808 {
809 for(Int_t i = 0; i < nGenSel; i++)
810 fhdPhiThrustGen->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
811 }
812
20d01a2f 813 if(eventShapes[0] <= 0.8 && nGenSel > 1)
95392a57 814 {
815 for(Int_t i = 0; i < nGenSel; i++)
816 fhdPhiThrustGenALL->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
817 }
818
819 Double_t fPxGen = vsumGen.Px();
820 Double_t fPyGen = vsumGen.Py();
821 Double_t fPzGen = vsumGen.Pz();
822 Double_t fEGen = vsumGen.E();
823
824 Double_t eRestGen[kMaxJets];
825 for (Int_t j = 0; j < nGenSel; j++)
826 {
827 vGen[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
828 eRestGen[j] = vGen[j].E();
829 }
830
831 for (Int_t j = 0; j < nAccTr; j++)
832 {
833 vTrackMC[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
834 pTrackMCBoost[j].SetXYZ(vTrackMC[j].Px(),vTrackMC[j].Py(),vTrackMC[j].Pz());
835 }
836
837 Int_t idxRestGen[kMaxJets];
838 TMath::Sort(nGenSel, eRestGen, idxRestGen);
839 for(Int_t j = 0; j < nGenSel; j++)
840 {
841 vRestGen[j] = vGen[idxRestGen[j]];
842 eSumGen += vRestGen[j].E();
843 }
844
845 if (nGenSel == 3)
846 {
847 // if(nInJet[0] < 3 || nInJet[1] < 3 || nInJet[2] < 3) return;
848 // if(pRestGen[1].DeltaPhi(pRestGen[2]) > 0.95 && pRestGen[1].DeltaPhi(pRestGen[2]) < 1.15)
849 // {
850
851 for(Int_t i = 0; i < nGenSel; i++)
852 {
853 xGen[i] = 2*vRestGen[i].E()/eSumGen;
854 }
855
856 Printf("***************** Values of Dalitz variables are : %f, %f, %f ****************\n", xGen[0], xGen[1], xGen[2]);
857
858 Printf("***************** fXSection = %f ******************\n", fXsection);
859 if(eSumGen <= 60)
860 fhX3X4Gen60->Fill(xGen[0], xGen[1], fXsection);
861
862 if(eSumGen > 60 && eSumGen <= 100)
863 fhX3X4Gen60100->Fill(xGen[0], xGen[1], fXsection);
864
865 if(eSumGen > 100)
866 fhX3X4Gen100->Fill(xGen[0], xGen[1], fXsection);
867
868 FillTopology(fhX3X4Gen, fhMu34Gen, fhMu45Gen, fhMu35Gen, xGen, pRestGen, fXsection);
869 }
870
871
872
873//_______________________________________________histos for MC_____________________________________________________
874
875
876//_______________________________________________histos AOD________________________________________________________
877
878// Printf("Event Number : %d, Number of gen jets : %d ", fEntry, nGenJets);
879
880 Int_t nRecSel = 0;
881 Int_t counter1 = 0;
882 Int_t tag1 = 0;
883
884 AliAODJet recSelJets[kMaxJets];
885
886 for(Int_t i = 0; i < nRecJets; i++)
887 {
888 if(nRecJets == 1)
889 {
890 recSelJets[nRecSel] = recJets[i];
891 nRecSel++;
892 }
893 else
894 {
895 counter1 = 0;
896 tag1 = 0;
897 for(Int_t j = 0; j < nRecJets; j++)
898 {
899 if(i!=j)
900 {
901 Double_t dRij = recJets[i].DeltaR(&recJets[j]);
902 counter1++;
903 if(dRij > 2*fR) tag1++;
904 }
905 }
906 if(counter1!=0)
907 {
908 if(tag1/counter1 == 1)
909 {
910 recSelJets[nRecSel] = recJets[i];
911 nRecSel++;
912 }
913 }
914 }
915 }
916
917 if(nRecSel == 0) return;
918
919 //sort rec/gen jets by energy in C.M.S
920 for (Int_t rj = 0; rj < nRecSel; rj++)
921 {
922 eRec[rj] = recSelJets[rj].E();
923 }
924
20d01a2f 925 Int_t nAODtracks = fAOD->GetNumberOfTracks();
95392a57 926 Int_t nTracks = 0; //tracks accepted in the whole event
20d01a2f 927 Int_t nTracksALL = 0;
95392a57 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;
20d01a2f 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();
95392a57 961 nTracks++;
962 nJetTracks++;
963 }
964 if(nJetTracks < 3) continue;
965 jetRecAcc[nAccJets] = jetRec[rj];
966 jetMult[nAccJets] = jetTracksAOD->GetEntries();
967 nAccJets++;
968 }
20d01a2f 969
95392a57 970 if (nAccJets == 0) return;
971
20d01a2f 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
95392a57 986 TLorentzVector vTrack[kTracks];
20d01a2f 987 TMath::Sort(nTracksALL, pTracks, idxTracks);
988 for(Int_t i = 0; i < nTracksALL; i++)
95392a57 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();
20d01a2f 1004 n01.SetZ(0.);
95392a57 1005 if(nAccJets > 1)
1006 {
1007// if(fGlobVar == 1)
1008// {
20d01a2f 1009 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
95392a57 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 }
95392a57 1075 }
1076
1077 //rest frame for reconstructed jets
1078 Double_t fPx = vsumRec.Px();
1079 Double_t fPy = vsumRec.Py();
1080 Double_t fPz = vsumRec.Pz();
1081 Double_t fE = vsumRec.E();
1082
1083 TVector3 pTrackRest[kTracks];
1084 for(Int_t j = 0; j < nTracks; j++)
1085 {
1086 vTrack[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
1087 pTrackRest[j].SetXYZ(vTrack[j].Px(), vTrack[j].Py(),vTrack[j].Pz());
1088 }
1089
1090 Double_t eRestRec[kMaxJets];
1091 Int_t idxRestRec[kMaxJets];
1092 for (Int_t j = 0; j < nAccJets; j++)
1093 {
1094 vRec[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
1095 eRestRec[j] = vRec[j].E();
1096 eSumRec += vRec[j].E();
1097 }
1098
1099 TMath::Sort(nAccJets, eRestRec, idxRestRec);
1100 for (Int_t i = 0; i < nAccJets; i++)
1101 {
1102 vRestRec[i] = vRec[idxRestRec[i]];
1103 pRestRec[i].SetXYZ(vRestRec[i].Px(), vRestRec[i].Py(), vRestRec[i].Pz());
1104 psumRestRec += pRestRec[i].Perp();
1105 }
1106
1107 if(nAccJets == 3)
1108 {
1109// if(pRest[1].DeltaPhi(pRest[2]) > 0.95 && pRest[1].DeltaPhi(pRest[2]) < 1.15)
1110// {
1111 fhInOut->Fill(nGenSel);
1112// for(Int_t j = 0; j < nTracksALL; j++)
1113// {
1114// vTracksAll[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
1115// pTracksAll[j].SetXYZ(vTracksAll[j].Px(), vTracksAll[j].Py(),vTracksAll[j].Pz());
1116// fhdPhiRec->Fill(pRest[0].DeltaPhi(pTracksAll[j]), pTracksAll[j].Perp(), fXsection);
1117// }
1118 //and the Dalitz variables and Energy distributions in the rest frame
1119 for (Int_t i = 0; i < nAccJets; i++)
1120 xRec[i] = 2*vRestRec[i].E()/eSumRec;
1121
1122 if(eSumRec <= 60)
1123 fhX3X4Rec60->Fill(xRec[0], xRec[1], fXsection);
1124
1125 if(eSumRec > 60 && eSumRec <= 100)
1126 fhX3X4Rec60100->Fill(xRec[0], xRec[1], fXsection);
1127
1128 if(eSumRec > 100)
1129 fhX3X4Rec100->Fill(xRec[0], xRec[1], fXsection);
1130
1131 if(nAccJets == 3 && nAccJets == nGenJets)
1132 {
1133 fhX3->Fill(xGen[0], TMath::Abs(xGen[0]-xRec[0])/xGen[0], fXsection);
1134 fhX4->Fill(xGen[1], TMath::Abs(xGen[1]-xRec[1])/xGen[1], fXsection);
1135 fhX5->Fill(xGen[2], TMath::Abs(xGen[2]-xRec[2])/xGen[2], fXsection);
1136 }
1137
1138 FillTopology(fhX3X4Rec, fhMu34Rec, fhMu45Rec, fhMu35Rec, xRec, pRestRec, fXsection);
1139 }
1140 Printf("%s:%d",(char*)__FILE__,__LINE__);
1141
1142 PostData(1, fList);
1143
1144 Printf("%s:%d Data Posted",(char*)__FILE__,__LINE__);
1145
1146}
1147
1148//__________________________________________________________________________________________________________________________________________________
1149void AliAnalysisTaskThreeJets::Terminate(Option_t *)
1150{
1151 printf("AnalysisJetCorrelation::Terminate()");
1152
1153}
1154
1155//_______________________________________User defined functions_____________________________________________________________________________________
1156void AliAnalysisTaskThreeJets::FillTopology(TH2F * Dalitz, TH1F * fhMu34, TH1F * fhMu45, TH1F * fhMu35, Double_t * x, TVector3 * pRest, Double_t xsection)
1157{
393afcc4 1158 //
1159 // fill the topology histos
1160 //
95392a57 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
1169Bool_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");
20d01a2f 1225 return kTRUE;
95392a57 1226 }
1227
1228 return kTRUE;
1229}
1230
1231//______________________________________________________________________________________________________
1232
47dcf196 1233
95392a57 1234//__________________________________________________________________________________________________________________________
1235
393afcc4 1236Double_t AliAnalysisTaskThreeJets::Exponent(Double_t x,const Double_t * const par) const
95392a57 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
393afcc4 1241Double_t AliAnalysisTaskThreeJets::Exponent2(Double_t x,const Double_t * const par) const
95392a57 1242{
1243 return par[0]*TMath::Power(1/TMath::E(), TMath::Power(par[1]/x, par[2]))+par[3]*x;
1244}
1245
393afcc4 1246Double_t AliAnalysisTaskThreeJets::Gauss(Double_t x,const Double_t * const par) const
95392a57 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
393afcc4 1251Double_t AliAnalysisTaskThreeJets::Total(Double_t x,const Double_t * const par) const
95392a57 1252{
1253 return Exponent(x, par)+Gauss(x, &par[4]);
1254}
1255
1256