]>
Commit | Line | Data |
---|---|---|
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 | ||
56 | ClassImp(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 | ||
137 | AliAnalysisTaskThreeJets::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 | ||
218 | Bool_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 | //___________________________________________________________________________________________________________________________________ | |
251 | void 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 | //__________________________________________________________________________________________________________________________________________ | |
453 | void AliAnalysisTaskThreeJets::Init() | |
454 | { | |
455 | printf("AliAnalysisJetCut::Init() \n"); | |
456 | } | |
457 | ||
458 | //____________________________________________________________________________________________________________________________________________ | |
459 | void 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 | //__________________________________________________________________________________________________________________________________________________ | |
1193 | void AliAnalysisTaskThreeJets::Terminate(Option_t *) | |
1194 | { | |
ec00a06a | 1195 | if(fDebug)printf(" AliAnalysisTaskThreeJets::Terminate()"); |
95392a57 | 1196 | |
1197 | } | |
1198 | ||
1199 | //_______________________________________User defined functions_____________________________________________________________________________________ | |
1200 | void 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 | ||
1213 | Bool_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 | 1280 | Double_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 | 1285 | Double_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 | 1290 | Double_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 | 1295 | Double_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 |