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