]> 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
coverity fixes
[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>
0eb7227a 39#include <TRandom.h>
95392a57 40
41#include "AliAnalysisTaskThreeJets.h"
42#include "AliAnalysisManager.h"
95392a57 43#include "AliAODEvent.h"
44#include "AliAODVertex.h"
45#include "AliAODHandler.h"
46#include "AliAODTrack.h"
47#include "AliAODJet.h"
e946cd3a 48#include "AliAODMCParticle.h"
49//#include "AliGenPythiaEventHeader.h"
50//#include "AliMCEvent.h"
51//#include "AliStack.h"
95392a57 52
95392a57 53
54#include "AliAnalysisHelperJetTasks.h"
55
56
57ClassImp(AliAnalysisTaskThreeJets)
58
59 AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets() : AliAnalysisTaskSE(),
60
61 fAOD(0x0),
7fa8b2da 62 fUseMC(0x0),
95392a57 63
64 fBranchRec(""),
65 fBranchGen(""),
66
67 fUseAODInput(kFALSE),
68
69 fR(0x0),
70 fList(0x0),
71
216601f0 72 fhStopHisto(0x0),
73
95392a57 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)
20d01a2f 133
95392a57 134{
135
136}
137
138AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets(const char * name):
139 AliAnalysisTaskSE(name),
140
141 fAOD(0x0),
7fa8b2da 142 fUseMC(0x0),
95392a57 143
144 fBranchRec(""),
145 fBranchGen(""),
146
147 fUseAODInput(kFALSE),
148
149 fR(0x0),
150 fList(0x0),
151
216601f0 152 fhStopHisto(0x0),
153
95392a57 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 //
5010a3f7 225
226
95392a57 227 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
5010a3f7 228 Float_t xsection = 0;
229 Float_t ftrials = 1;
230
231 Float_t fAvgTrials = 1;
95392a57 232 if(tree){
233 TFile *curfile = tree->GetCurrentFile();
234 if (!curfile) {
235 Error("Notify","No current file");
236 return kFALSE;
237 }
5010a3f7 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;
95392a57 245
95392a57 246 return kTRUE;
5010a3f7 247
95392a57 248}
249
250
251//___________________________________________________________________________________________________________________________________
252void AliAnalysisTaskThreeJets::UserCreateOutputObjects()
253{
254 //
255 // Create the output container
256 //
47dcf196 257 // Printf("Analysing event %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
95392a57 258
95392a57 259 printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
260
261 fList = new TList();
262
216601f0 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
95392a57 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
ec00a06a 450 if(fDebug)Printf("UserCreateOutputObjects finished\n");
95392a57 451}
452
453//__________________________________________________________________________________________________________________________________________
454void AliAnalysisTaskThreeJets::Init()
455{
456 printf("AliAnalysisJetCut::Init() \n");
457}
458
459//____________________________________________________________________________________________________________________________________________
460void AliAnalysisTaskThreeJets::UserExec(Option_t * )
461{
5010a3f7 462 if (fDebug > 1) printf("AliAnlysisTaskThreeJets::Analysing event # %5d\n", (Int_t) fEntry);
95392a57 463
5010a3f7 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);
95392a57 468 return;
5010a3f7 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
e946cd3a 480// AliMCEvent* mcEvent =MCEvent();
481// if(!mcEvent){
482// Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
483// return;
484// }
95392a57 485
486 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
487
488 //primary vertex
489 AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
5010a3f7 490 if(!pvtx){
216601f0 491 //return #1
492 fhStopHisto->Fill(0.5);
493 PostData(1, fList);
5010a3f7 494 return;
e946cd3a 495 }
496
497// AliAODJet genJetsPythia[kMaxJets];
498// Int_t nPythiaGenJets = 0;
7fa8b2da 499
95392a57 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){
216601f0 506 //return #2
507 fhStopHisto->Fill(1.5);
508 PostData(1, fList);
95392a57 509 return;
510 }
511
512 // reconstructed jets
513 nRecJets = aodRecJets->GetEntries();
a923bd34 514 if(fDebug)Printf("--- Jets found in bRec: %d", nRecJets);
95392a57 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 }
7fa8b2da 523
524 AliAODJet genJets[kMaxJets];
525 Int_t nGenJets = 0;
526 if(fUseMC){
95392a57 527 // If we set a second branch for the input jets fetch this
528 TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
7fa8b2da 529
530 if(!aodGenJets)
531 {
216601f0 532 printf("NO MC jets Found\n");
533 //return #3
534 fhStopHisto->Fill(2.5);
535 PostData(1, fList);
7fa8b2da 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);
e946cd3a 551// if(!pythiaGenHeader){
552// Printf("!!!NO GEN HEADER AVALABLE!!!");
553// return;
554// }
95392a57 555
e946cd3a 556// // Int_t ProcessType = pythiaGenHeader->ProcessType();
557// // if(ProcessType != 28) return;
558// nPythiaGenJets = pythiaGenHeader->NTriggerJets();
559// nPythiaGenJets = TMath::Min(nPythiaGenJets, kMaxJets);
95392a57 560
e946cd3a 561// fXsection = 1;
562
95392a57 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;
e946cd3a 600// //Pythia_________________________________________________________________________________________________________________
95392a57 601
e946cd3a 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// }
95392a57 609
610//_________________________________________________________________________________________________________________________
611
612
613//________histos for MC___________________________________________________________________________________________________________
95392a57 614 Int_t nGenSel = 0;
7fa8b2da 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 }
95392a57 651
216601f0 652 if(nGenSel == 0){
653 //return #4
654 fhStopHisto->Fill(3.5);
655 PostData(1, fList);
656 return;
657 }
7fa8b2da 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];
0eb7227a 684 Double_t eventShapes[4] = {0,};
7fa8b2da 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++)
95392a57 700 {
7fa8b2da 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++;
95392a57 710 }
216601f0 711 if(nAllTracksMC == 0){
712 //return #5
713 fhStopHisto->Fill(4.5);
714 PostData(1, fList);
715 return;
716 }
717
7fa8b2da 718 for(Int_t iJet = 0; iJet < nGenSel; iJet++)
95392a57 719 {
7fa8b2da 720 Int_t nJetTracks = 0;
721 for(Int_t i = 0; i < nAllTracksMC; i++)
95392a57 722 {
7fa8b2da 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)
95392a57 729 {
7fa8b2da 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++;
95392a57 736 }
737 }
7fa8b2da 738 nInJet[iJet] = nJetTracks;
739 }
740
216601f0 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);
7fa8b2da 748 Double_t pTav[kMaxJets];
749 for(Int_t i = 0; i < nGenSel; i++)
750 {
751 Double_t pTsum = 0;
216601f0 752 // if(fDebug)Printf("*********** Number of particles in Jet %d = %d *******************\n", i+3, nInJet[i]);
7fa8b2da 753 for(Int_t iT = 0; iT < nInJet[i]; iT++)
95392a57 754 {
7fa8b2da 755 Double_t pt = inJetPartV[i][iT].Pt();
756 pTsum += pt;
95392a57 757 }
7fa8b2da 758 pTav[i] = pTsum/nInJet[i];
95392a57 759 }
7fa8b2da 760
761 TMath::Sort(nAllTracksMC, pTrackMCAll, idxTracksMC);
95392a57 762 for(Int_t i = 0; i < nAllTracksMC; i++)
763 {
7fa8b2da 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()){
0eb7227a 781 Double_t eventShapesTest[4] = {0,};
7fa8b2da 782 TVector3 n01Test;
783 Int_t rnd_max = nAllTracksMC;
0eb7227a 784 Int_t k = (Int_t)(gRandom->Rndm()*rnd_max)+3;
7fa8b2da 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)
95392a57 800 {
7fa8b2da 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;
95392a57 826 }
827 }
95392a57 828 }
829
7fa8b2da 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];
20d01a2f 837 }
7fa8b2da 838 if(tca){
839 if(eventShapes[0] > 0.8 && nGenSel > 1)
95392a57 840 {
7fa8b2da 841 for(Int_t i = 0; i < nGenSel; i++)
842 fhdPhiThrustGen->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
95392a57 843 }
7fa8b2da 844
845 if(eventShapes[0] <= 0.8 && nGenSel > 1)
95392a57 846 {
7fa8b2da 847 for(Int_t i = 0; i < nGenSel; i++)
848 fhdPhiThrustGenALL->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
95392a57 849 }
850 }
7fa8b2da 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 }
95392a57 877
7fa8b2da 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 // {
95392a57 883
7fa8b2da 884 for(Int_t i = 0; i < nGenSel; i++)
885 {
886 xGen[i] = 2*vRestGen[i].E()/eSumGen;
95392a57 887 }
888
7fa8b2da 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
95392a57 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];
a923bd34 919 if(fDebug)Printf("---- Number of reco jets: %d\n",nRecJets);
95392a57 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
216601f0 951 if(nRecSel == 0)
952 {
953 //return #7
954 fhStopHisto->Fill(6.5);
955 PostData(1, fList);
956 return;
957 }
958
95392a57 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
20d01a2f 965 Int_t nAODtracks = fAOD->GetNumberOfTracks();
95392a57 966 Int_t nTracks = 0; //tracks accepted in the whole event
20d01a2f 967 Int_t nTracksALL = 0;
95392a57 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];
0eb7227a 973 Double_t eventShapesRec[4] = {0,};
95392a57 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;
20d01a2f 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();
95392a57 1001 nTracks++;
1002 nJetTracks++;
1003 }
1004 if(nJetTracks < 3) continue;
1005 jetRecAcc[nAccJets] = jetRec[rj];
1006 jetMult[nAccJets] = jetTracksAOD->GetEntries();
1007 nAccJets++;
1008 }
20d01a2f 1009
216601f0 1010 if (nAccJets == 0){
1011 //return #8
1012 fhStopHisto->Fill(7.5);
1013 PostData(1, fList);
1014 return;
1015 }
95392a57 1016
20d01a2f 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
95392a57 1031 TLorentzVector vTrack[kTracks];
20d01a2f 1032 TMath::Sort(nTracksALL, pTracks, idxTracks);
1033 for(Int_t i = 0; i < nTracksALL; i++)
95392a57 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();
20d01a2f 1049 n01.SetZ(0.);
95392a57 1050 if(nAccJets > 1)
1051 {
1052// if(fGlobVar == 1)
1053// {
ec00a06a 1054 if(fDebug)Printf("*********Shapes for AOD********");
20d01a2f 1055 AliAnalysisHelperJetTasks::GetEventShapes(n01, pTrack, nTracksALL, eventShapesRec);
1056 // }
1057 if(eventShapesRec[0] < 2/TMath::Pi()){
0eb7227a 1058 Double_t eventShapesTest[4] = {0,};
20d01a2f 1059 TVector3 n01Test;
1060 Int_t rnd_max = nTracksALL;
0eb7227a 1061 Int_t k = (Int_t)(gRandom->Rndm()*rnd_max)+3;
20d01a2f 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
95392a57 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 }
95392a57 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// {
7fa8b2da 1156 if(fUseMC) fhInOut->Fill(nGenSel);
95392a57 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 }
ec00a06a 1185 if(fDebug)Printf("%s:%d",(char*)__FILE__,__LINE__);
95392a57 1186
1187 PostData(1, fList);
1188
ec00a06a 1189 if(fDebug)Printf("%s:%d Data Posted",(char*)__FILE__,__LINE__);
95392a57 1190
1191}
1192
1193//__________________________________________________________________________________________________________________________________________________
1194void AliAnalysisTaskThreeJets::Terminate(Option_t *)
1195{
ec00a06a 1196 if(fDebug)printf(" AliAnalysisTaskThreeJets::Terminate()");
95392a57 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{
393afcc4 1203 //
1204 // fill the topology histos
1205 //
95392a57 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");
20d01a2f 1270 return kTRUE;
95392a57 1271 }
1272
1273 return kTRUE;
1274}
1275
1276//______________________________________________________________________________________________________
1277
47dcf196 1278
95392a57 1279//__________________________________________________________________________________________________________________________
1280
393afcc4 1281Double_t AliAnalysisTaskThreeJets::Exponent(Double_t x,const Double_t * const par) const
95392a57 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
393afcc4 1286Double_t AliAnalysisTaskThreeJets::Exponent2(Double_t x,const Double_t * const par) const
95392a57 1287{
1288 return par[0]*TMath::Power(1/TMath::E(), TMath::Power(par[1]/x, par[2]))+par[3]*x;
1289}
1290
393afcc4 1291Double_t AliAnalysisTaskThreeJets::Gauss(Double_t x,const Double_t * const par) const
95392a57 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
393afcc4 1296Double_t AliAnalysisTaskThreeJets::Total(Double_t x,const Double_t * const par) const
95392a57 1297{
1298 return Exponent(x, par)+Gauss(x, &par[4]);
1299}
1300
1301