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