]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/RESONANCES/AliRsnMiniAnalysisTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
st^Cic functions for p-value: cleanup and add all code that was used for approval
[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / AliRsnMiniAnalysisTask.cxx
CommitLineData
03d23846 1//
2// Analysis task for 'mini' sub-package
3// Contains all definitions needed for running an analysis:
4// -- global event cut
5// -- a list of track cuts (any number)
6// -- definitions of output histograms
7// -- values to be computed.
8// Each one must be defined using the "CREATE" methods, which
9// add directly a new element in the task collections, and don't
10// need an external object to be passed to the task itself.
11//
12
13#include <Riostream.h>
14
03d23846 15#include <TH1.h>
6aaeb33c 16#include <TList.h>
17#include <TTree.h>
52e6652d 18#include <TStopwatch.h>
03d23846 19
20#include "AliLog.h"
6aaeb33c 21#include "AliEventplane.h"
22#include "AliMultiplicity.h"
23#include "AliTriggerAnalysis.h"
03d23846 24#include "AliAnalysisManager.h"
6aaeb33c 25#include "AliInputEventHandler.h"
26
03d23846 27#include "AliESDtrackCuts.h"
28#include "AliESDUtils.h"
6aaeb33c 29
03d23846 30#include "AliAODEvent.h"
03d23846 31#include "AliAODMCParticle.h"
32
33#include "AliRsnCutSet.h"
6aaeb33c 34#include "AliRsnMiniPair.h"
03d23846 35#include "AliRsnMiniEvent.h"
36#include "AliRsnMiniParticle.h"
03d23846 37
38#include "AliRsnMiniAnalysisTask.h"
39
52e6652d 40
03d23846 41ClassImp(AliRsnMiniAnalysisTask)
42
43//__________________________________________________________________________________________________
44AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() :
45 AliAnalysisTaskSE(),
c80938f5 46 fUseMC(kFALSE),
03d23846 47 fEvNum(0),
4fb0dfa3 48 fTriggerMask(0),
03d23846 49 fUseCentrality(kFALSE),
50 fCentralityType("QUALITY"),
1cdc3671 51 fUseAOD049CentralityPatch(kFALSE),
caef2e3c 52 fContinuousMix(kTRUE),
03d23846 53 fNMix(0),
54 fMaxDiffMult(10),
55 fMaxDiffVz(1.0),
56 fMaxDiffAngle(1E20),
57 fOutput(0x0),
58 fHistograms("AliRsnMiniOutput", 0),
59 fValues("AliRsnMiniValue", 0),
60 fHEventStat(0x0),
955ca046 61 fHAEventsVsMulti(0x0),
5654276f 62 fHAEventVz(0x0),
63 fHAEventMultiCent(0x0),
64 fHAEventPlane(0x0),
03d23846 65 fEventCuts(0x0),
66 fTrackCuts(0),
67 fRsnEvent(),
c80938f5 68 fEvBuffer(0x0),
c80938f5 69 fTriggerAna(0x0),
6aaeb33c 70 fESDtrackCuts(0x0),
caef2e3c 71 fMiniEvent(0x0),
61f275d1 72 fBigOutput(kFALSE),
088ca370 73 fMixPrintRefresh(-1),
74 fMaxNDaughters(-1),
31dbef4e 75 fCheckP(kFALSE),
76 fCheckFeedDown(kFALSE),
77 fOriginDselection(kFALSE),
78 fKeepDfromB(kFALSE),
79 fKeepDfromBOnly(kFALSE),
80 fRejectIfNoQuark(kFALSE)
03d23846 81{
82//
83// Dummy constructor ALWAYS needed for I/O.
84//
85}
86
87//__________________________________________________________________________________________________
c80938f5 88AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name, Bool_t useMC) :
03d23846 89 AliAnalysisTaskSE(name),
c80938f5 90 fUseMC(useMC),
03d23846 91 fEvNum(0),
4fb0dfa3 92 fTriggerMask(AliVEvent::kMB),
03d23846 93 fUseCentrality(kFALSE),
94 fCentralityType("QUALITY"),
1cdc3671 95 fUseAOD049CentralityPatch(kFALSE),
caef2e3c 96 fContinuousMix(kTRUE),
03d23846 97 fNMix(0),
98 fMaxDiffMult(10),
99 fMaxDiffVz(1.0),
100 fMaxDiffAngle(1E20),
101 fOutput(0x0),
102 fHistograms("AliRsnMiniOutput", 0),
103 fValues("AliRsnMiniValue", 0),
104 fHEventStat(0x0),
955ca046 105 fHAEventsVsMulti(0x0),
5654276f 106 fHAEventVz(0x0),
107 fHAEventMultiCent(0x0),
108 fHAEventPlane(0x0),
03d23846 109 fEventCuts(0x0),
110 fTrackCuts(0),
111 fRsnEvent(),
c80938f5 112 fEvBuffer(0x0),
c80938f5 113 fTriggerAna(0x0),
6aaeb33c 114 fESDtrackCuts(0x0),
caef2e3c 115 fMiniEvent(0x0),
61f275d1 116 fBigOutput(kFALSE),
088ca370 117 fMixPrintRefresh(-1),
118 fMaxNDaughters(-1),
31dbef4e 119 fCheckP(kFALSE),
120 fCheckFeedDown(kFALSE),
121 fOriginDselection(kFALSE),
122 fKeepDfromB(kFALSE),
123 fKeepDfromBOnly(kFALSE),
124 fRejectIfNoQuark(kFALSE)
03d23846 125{
126//
127// Default constructor.
128// Define input and output slots here (never in the dummy constructor)
129// Input slot #0 works with a TChain - it is connected to the default input container
130// Output slot #1 writes into a TH1 container
131//
132
133 DefineOutput(1, TList::Class());
134}
135
136//__________________________________________________________________________________________________
52e6652d 137AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask &copy) :
03d23846 138 AliAnalysisTaskSE(copy),
c80938f5 139 fUseMC(copy.fUseMC),
03d23846 140 fEvNum(0),
4fb0dfa3 141 fTriggerMask(copy.fTriggerMask),
03d23846 142 fUseCentrality(copy.fUseCentrality),
143 fCentralityType(copy.fCentralityType),
1cdc3671 144 fUseAOD049CentralityPatch(copy.fUseAOD049CentralityPatch),
caef2e3c 145 fContinuousMix(copy.fContinuousMix),
03d23846 146 fNMix(copy.fNMix),
147 fMaxDiffMult(copy.fMaxDiffMult),
148 fMaxDiffVz(copy.fMaxDiffVz),
149 fMaxDiffAngle(copy.fMaxDiffAngle),
150 fOutput(0x0),
151 fHistograms(copy.fHistograms),
152 fValues(copy.fValues),
153 fHEventStat(0x0),
955ca046 154 fHAEventsVsMulti(0x0),
5654276f 155 fHAEventVz(0x0),
156 fHAEventMultiCent(0x0),
157 fHAEventPlane(0x0),
03d23846 158 fEventCuts(copy.fEventCuts),
159 fTrackCuts(copy.fTrackCuts),
160 fRsnEvent(),
c80938f5 161 fEvBuffer(0x0),
c80938f5 162 fTriggerAna(copy.fTriggerAna),
6aaeb33c 163 fESDtrackCuts(copy.fESDtrackCuts),
caef2e3c 164 fMiniEvent(0x0),
61f275d1 165 fBigOutput(copy.fBigOutput),
088ca370 166 fMixPrintRefresh(copy.fMixPrintRefresh),
167 fMaxNDaughters(copy.fMaxNDaughters),
31dbef4e 168 fCheckP(copy.fCheckP),
169 fCheckFeedDown(copy.fCheckFeedDown),
170 fOriginDselection(copy.fOriginDselection),
171 fKeepDfromB(copy.fOriginDselection),
172 fKeepDfromBOnly(copy.fKeepDfromBOnly),
173 fRejectIfNoQuark(copy.fRejectIfNoQuark)
03d23846 174{
175//
176// Copy constructor.
177// Implemented as requested by C++ standards.
178// Can be used in PROOF and by plugins.
179//
180}
181
182//__________________________________________________________________________________________________
52e6652d 183AliRsnMiniAnalysisTask &AliRsnMiniAnalysisTask::operator=(const AliRsnMiniAnalysisTask &copy)
03d23846 184{
185//
186// Assignment operator.
187// Implemented as requested by C++ standards.
188// Can be used in PROOF and by plugins.
189//
190
191 AliAnalysisTaskSE::operator=(copy);
e6f3a909 192 if (this == &copy)
61f275d1 193 return *this;
c80938f5 194 fUseMC = copy.fUseMC;
4fb0dfa3 195 fEvNum = copy.fEvNum;
196 fTriggerMask = copy.fTriggerMask;
03d23846 197 fUseCentrality = copy.fUseCentrality;
198 fCentralityType = copy.fCentralityType;
1cdc3671 199 fUseAOD049CentralityPatch = copy.fUseAOD049CentralityPatch;
caef2e3c 200 fContinuousMix = copy.fContinuousMix;
03d23846 201 fNMix = copy.fNMix;
202 fMaxDiffMult = copy.fMaxDiffMult;
203 fMaxDiffVz = copy.fMaxDiffVz;
204 fMaxDiffAngle = copy.fMaxDiffAngle;
205 fHistograms = copy.fHistograms;
206 fValues = copy.fValues;
955ca046 207 fHEventStat = copy.fHEventStat;
208 fHAEventsVsMulti = copy.fHAEventsVsMulti;
5654276f 209 fHAEventVz = copy.fHAEventVz;
210 fHAEventMultiCent = copy.fHAEventMultiCent;
211 fHAEventPlane = copy.fHAEventPlane;
03d23846 212 fEventCuts = copy.fEventCuts;
213 fTrackCuts = copy.fTrackCuts;
214 fTriggerAna = copy.fTriggerAna;
c80938f5 215 fESDtrackCuts = copy.fESDtrackCuts;
caef2e3c 216 fBigOutput = copy.fBigOutput;
61f275d1 217 fMixPrintRefresh = copy.fMixPrintRefresh;
088ca370 218 fMaxNDaughters = copy.fMaxNDaughters;
219 fCheckP = copy.fCheckP;
31dbef4e 220 fCheckFeedDown = copy.fCheckFeedDown;
221 fOriginDselection = copy.fOriginDselection;
222 fKeepDfromB = copy.fOriginDselection;
223 fKeepDfromBOnly = copy.fKeepDfromBOnly;
224 fRejectIfNoQuark = copy.fRejectIfNoQuark;
03d23846 225 return (*this);
226}
227
228//__________________________________________________________________________________________________
229AliRsnMiniAnalysisTask::~AliRsnMiniAnalysisTask()
230{
231//
52e6652d 232// Destructor.
03d23846 233// Clean-up the output list, but not the histograms that are put inside
234// (the list is owner and will clean-up these histograms). Protect in PROOF case.
235//
236
d573d2fb 237
03d23846 238 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
239 delete fOutput;
d573d2fb 240 delete fEvBuffer;
03d23846 241 }
242}
243
244//__________________________________________________________________________________________________
245Int_t AliRsnMiniAnalysisTask::AddTrackCuts(AliRsnCutSet *cuts)
246{
247//
248// Add a new cut set for a new criterion for track selection.
249// A user can add as many as he wants, and each one corresponds
250// to one of the available bits in the AliRsnMiniParticle mask.
251// The only check is the following: if a cut set with the same name
252// as the argument is there, this is not added.
253// Return value is the array position of this set.
254//
255
256 TObject *obj = fTrackCuts.FindObject(cuts->GetName());
52e6652d 257
03d23846 258 if (obj) {
259 AliInfo(Form("A cut set named '%s' already exists", cuts->GetName()));
260 return fTrackCuts.IndexOf(obj);
261 } else {
262 fTrackCuts.AddLast(cuts);
263 return fTrackCuts.IndexOf(cuts);
264 }
265}
266
267//__________________________________________________________________________________________________
268void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
269{
270//
271// Initialization of outputs.
272// This is called once per worker node.
273//
274
275 // reset counter
9e7b94f5 276 fEvNum = -1;
52e6652d 277
d573d2fb 278 // message
279 AliInfo(Form("Selected event characterization: %s (%s)", (fUseCentrality ? "centrality" : "multiplicity"), fCentralityType.Data()));
03d23846 280
281 // initialize trigger analysis
c80938f5 282 if (fTriggerAna) delete fTriggerAna;
03d23846 283 fTriggerAna = new AliTriggerAnalysis;
52e6652d 284
c80938f5 285 // initialize ESD quality cuts
286 if (fESDtrackCuts) delete fESDtrackCuts;
287 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
03d23846 288
289 // create list and set it as owner of its content (MANDATORY)
caef2e3c 290 if (fBigOutput) OpenFile(1);
03d23846 291 fOutput = new TList();
292 fOutput->SetOwner();
52e6652d 293
03d23846 294 // initialize event statistics counter
295 fHEventStat = new TH1F("hEventStat", "Event statistics", 4, 0.0, 4.0);
296 fHEventStat->GetXaxis()->SetBinLabel(1, "CINT1B");
297 fHEventStat->GetXaxis()->SetBinLabel(2, "V0AND");
298 fHEventStat->GetXaxis()->SetBinLabel(3, "Candle");
299 fHEventStat->GetXaxis()->SetBinLabel(4, "Accepted");
300 fOutput->Add(fHEventStat);
52e6652d 301
955ca046 302 if (fUseCentrality)
303 fHAEventsVsMulti = new TH1F("hAEventsVsMulti", "Accepted events vs Centrality", 100, 0, 100.0);
304 else
e6952ec7 305 fHAEventsVsMulti = new TH1F("hAEventsVsMulti", "Accepted events vs Multiplicity",1000, 0, 1000.0);
955ca046 306 fOutput->Add(fHAEventsVsMulti);
e6952ec7 307
5654276f 308 if(fHAEventVz) fOutput->Add(fHAEventVz);
309 if(fHAEventMultiCent) fOutput->Add(fHAEventMultiCent);
310 if(fHAEventPlane) fOutput->Add(fHAEventPlane);
311
61f275d1 312 TIter next(&fTrackCuts);
313 AliRsnCutSet *cs;
314 while ((cs = (AliRsnCutSet *) next())) {
315 cs->Init(fOutput);
316 }
317
03d23846 318 // create temporary tree for filtered events
6aaeb33c 319 if (fMiniEvent) delete fMiniEvent;
c80938f5 320 fEvBuffer = new TTree("EventBuffer", "Temporary buffer for mini events");
6aaeb33c 321 fEvBuffer->Branch("events", "AliRsnMiniEvent", &fMiniEvent);
52e6652d 322
03d23846 323 // create one histogram per each stored definition (event histograms)
324 Int_t i, ndef = fHistograms.GetEntries();
325 AliRsnMiniOutput *def = 0x0;
326 for (i = 0; i < ndef; i++) {
52e6652d 327 def = (AliRsnMiniOutput *)fHistograms[i];
03d23846 328 if (!def) continue;
329 if (!def->Init(GetName(), fOutput)) {
330 AliError(Form("Def '%s': failed initialization", def->GetName()));
331 continue;
332 }
333 }
52e6652d 334
03d23846 335 // post data for ALL output slots >0 here, to get at least an empty histogram
336 PostData(1, fOutput);
337}
338
339//__________________________________________________________________________________________________
340void AliRsnMiniAnalysisTask::UserExec(Option_t *)
341{
342//
6aaeb33c 343// Computation loop.
344// In this case, it checks if the event is acceptable, and eventually
345// creates the corresponding mini-event and stores it in the buffer.
346// The real histogram filling is done at the end, in "FinishTaskOutput".
03d23846 347//
348
9e7b94f5 349 // increment event counter
03d23846 350 fEvNum++;
9e7b94f5 351
03d23846 352 // check current event
353 Char_t check = CheckCurrentEvent();
354 if (!check) return;
52e6652d 355
d573d2fb 356 // setup PID response
52e6652d 357 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
358 AliInputEventHandler *inputHandler = (AliInputEventHandler *)man->GetInputEventHandler();
359 fRsnEvent.SetPIDResponse(inputHandler->GetPIDResponse());
360
6aaeb33c 361 // fill a mini-event from current
362 // and skip this event if no tracks were accepted
45aa62b9 363 FillMiniEvent(check);
52e6652d 364
6aaeb33c 365 // fill MC based histograms on mothers,
366 // which do need the original event
c80938f5 367 if (fUseMC) {
368 if (fRsnEvent.IsESD() && fMCEvent)
6aaeb33c 369 FillTrueMotherESD(fMiniEvent);
c80938f5 370 else if (fRsnEvent.IsAOD() && fRsnEvent.GetAODList())
6aaeb33c 371 FillTrueMotherAOD(fMiniEvent);
03d23846 372 }
52e6652d 373
9e7b94f5 374 // if the event is not empty, store it
375 if (fMiniEvent->IsEmpty()) {
376 AliDebugClass(2, Form("Rejecting empty event #%d", fEvNum));
377 } else {
378 Int_t id = fEvBuffer->GetEntries();
379 AliDebugClass(2, Form("Adding event #%d with ID = %d", fEvNum, id));
380 fMiniEvent->ID() = id;
381 fEvBuffer->Fill();
382 }
52e6652d 383
6aaeb33c 384 // post data for computed stuff
03d23846 385 PostData(1, fOutput);
386}
387
388//__________________________________________________________________________________________________
389void AliRsnMiniAnalysisTask::FinishTaskOutput()
390{
391//
6aaeb33c 392// This function is called at the end of the loop on available events,
393// and then the buffer will be full with all the corresponding mini-events,
394// each one containing all tracks selected by each of the available track cuts.
395// Here a loop is done on each of these events, and both single-event and mixing are computed
03d23846 396//
397
6aaeb33c 398 // security code: reassign the buffer to the mini-event cursor
399 fEvBuffer->SetBranchAddress("events", &fMiniEvent);
52e6652d 400 TStopwatch timer;
0efb7042 401 // prepare variables
6aaeb33c 402 Int_t ievt, nEvents = (Int_t)fEvBuffer->GetEntries();
403 Int_t idef, nDefs = fHistograms.GetEntries();
45aa62b9 404 Int_t imix, iloop, ifill;
6aaeb33c 405 AliRsnMiniOutput *def = 0x0;
406 AliRsnMiniOutput::EComputation compType;
52e6652d 407
61f275d1 408 Int_t printNum = fMixPrintRefresh;
409 if (printNum < 0) {
410 if (nEvents>1e5) printNum=nEvents/100;
411 else if (nEvents>1e4) printNum=nEvents/10;
412 else printNum = 0;
413 }
414
6aaeb33c 415 // loop on events, and for each one fill all outputs
416 // using the appropriate procedure depending on its type
417 // only mother-related histograms are filled in UserExec,
418 // since they require direct access to MC event
52e6652d 419 timer.Start();
6aaeb33c 420 for (ievt = 0; ievt < nEvents; ievt++) {
421 // get next entry
422 fEvBuffer->GetEntry(ievt);
52e6652d 423 if (printNum&&(ievt%printNum==0)) {
424 AliInfo(Form("[%s] Std.Event %d/%d",GetName(), ievt,nEvents));
61f275d1 425 timer.Stop(); timer.Print(); fflush(stdout); timer.Start(kFALSE);
52e6652d 426 }
45aa62b9 427 // fill
6aaeb33c 428 for (idef = 0; idef < nDefs; idef++) {
52e6652d 429 def = (AliRsnMiniOutput *)fHistograms[idef];
6aaeb33c 430 if (!def) continue;
431 compType = def->GetComputation();
432 // execute computation in the appropriate way
433 switch (compType) {
434 case AliRsnMiniOutput::kEventOnly:
9e7b94f5 435 //AliDebugClass(1, Form("Event %d, def '%s': event-value histogram filling", ievt, def->GetName()));
45aa62b9 436 ifill = 1;
0efb7042 437 def->FillEvent(fMiniEvent, &fValues);
6aaeb33c 438 break;
439 case AliRsnMiniOutput::kTruePair:
9e7b94f5 440 //AliDebugClass(1, Form("Event %d, def '%s': true-pair histogram filling", ievt, def->GetName()));
0efb7042 441 ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
6aaeb33c 442 break;
443 case AliRsnMiniOutput::kTrackPair:
9e7b94f5 444 //AliDebugClass(1, Form("Event %d, def '%s': pair-value histogram filling", ievt, def->GetName()));
0efb7042 445 ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
6aaeb33c 446 break;
447 case AliRsnMiniOutput::kTrackPairRotated1:
9e7b94f5 448 //AliDebugClass(1, Form("Event %d, def '%s': rotated (1) background histogram filling", ievt, def->GetName()));
0efb7042 449 ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
6aaeb33c 450 break;
451 case AliRsnMiniOutput::kTrackPairRotated2:
9e7b94f5 452 //AliDebugClass(1, Form("Event %d, def '%s': rotated (2) background histogram filling", ievt, def->GetName()));
0efb7042 453 ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
6aaeb33c 454 break;
455 default:
456 // other kinds are processed elsewhere
45aa62b9 457 ifill = 0;
6aaeb33c 458 AliDebugClass(2, Form("Computation = %d", (Int_t)compType));
459 }
45aa62b9 460 // message
461 AliDebugClass(1, Form("Event %6d: def = '%15s' -- fills = %5d", ievt, def->GetName(), ifill));
03d23846 462 }
463 }
52e6652d 464
0efb7042 465 // if no mixing is required, stop here and post the output
466 if (fNMix < 1) {
467 AliDebugClass(2, "Stopping here, since no mixing is required");
468 PostData(1, fOutput);
469 return;
470 }
52e6652d 471
0efb7042 472 // initialize mixing counter
9e7b94f5 473 Int_t nmatched[nEvents];
474 TString *smatched = new TString[nEvents];
0efb7042 475 for (ievt = 0; ievt < nEvents; ievt++) {
9e7b94f5 476 smatched[ievt] = "|";
477 nmatched[ievt] = 0;
0efb7042 478 }
61f275d1 479
480
52e6652d 481 AliInfo(Form("[%s] Std.Event %d/%d",GetName(), nEvents,nEvents));
61f275d1 482 timer.Stop(); timer.Print(); timer.Start(); fflush(stdout);
52e6652d 483
0efb7042 484 // search for good matchings
485 for (ievt = 0; ievt < nEvents; ievt++) {
52e6652d 486 if (printNum&&(ievt%printNum==0)) {
487 AliInfo(Form("[%s] EventMixing searching %d/%d",GetName(),ievt,nEvents));
488 timer.Stop(); timer.Print(); timer.Start(kFALSE); fflush(stdout);
489 }
9e7b94f5 490 if (nmatched[ievt] >= fNMix) continue;
0efb7042 491 fEvBuffer->GetEntry(ievt);
492 AliRsnMiniEvent evMain(*fMiniEvent);
493 for (iloop = 1; iloop < nEvents; iloop++) {
494 imix = ievt + iloop;
495 if (imix >= nEvents) imix -= nEvents;
496 if (imix == ievt) continue;
497 // text next entry
498 fEvBuffer->GetEntry(imix);
499 // skip if events are not matched
caef2e3c 500 if (!EventsMatch(&evMain, fMiniEvent)) continue;
0efb7042 501 // check that the array of good matches for mixed does not already contain main event
9e7b94f5 502 if (smatched[imix].Contains(Form("|%d|", ievt))) continue;
503 // check that the found good events has not enough matches already
504 if (nmatched[imix] >= fNMix) continue;
0efb7042 505 // add new mixing candidate
9e7b94f5 506 smatched[ievt].Append(Form("%d|", imix));
507 nmatched[ievt]++;
508 nmatched[imix]++;
509 if (nmatched[ievt] >= fNMix) break;
0efb7042 510 }
9e7b94f5 511 AliDebugClass(1, Form("Matches for event %5d = %d [%s] (missing are declared above)", evMain.ID(), nmatched[ievt], smatched[ievt].Data()));
0efb7042 512 }
61f275d1 513
52e6652d 514 AliInfo(Form("[%s] EventMixing searching %d/%d",GetName(),nEvents,nEvents));
61f275d1 515 timer.Stop(); timer.Print(); fflush(stdout); timer.Start();
52e6652d 516
0efb7042 517 // perform mixing
518 TObjArray *list = 0x0;
519 TObjString *os = 0x0;
520 for (ievt = 0; ievt < nEvents; ievt++) {
52e6652d 521 if (printNum&&(ievt%printNum==0)) {
522 AliInfo(Form("[%s] EventMixing %d/%d",GetName(),ievt,nEvents));
523 timer.Stop(); timer.Print(); timer.Start(kFALSE); fflush(stdout);
524 }
0efb7042 525 ifill = 0;
526 fEvBuffer->GetEntry(ievt);
527 AliRsnMiniEvent evMain(*fMiniEvent);
9e7b94f5 528 list = smatched[ievt].Tokenize("|");
0efb7042 529 TObjArrayIter next(list);
52e6652d 530 while ( (os = (TObjString *)next()) ) {
0efb7042 531 imix = os->GetString().Atoi();
0efb7042 532 fEvBuffer->GetEntry(imix);
533 for (idef = 0; idef < nDefs; idef++) {
52e6652d 534 def = (AliRsnMiniOutput *)fHistograms[idef];
0efb7042 535 if (!def) continue;
536 if (!def->IsTrackPairMix()) continue;
d573d2fb 537 ifill += def->FillPair(&evMain, fMiniEvent, &fValues, kTRUE);
9e7b94f5 538 if (!def->IsSymmetric()) {
539 AliDebugClass(2, "Reflecting non symmetric pair");
540 ifill += def->FillPair(fMiniEvent, &evMain, &fValues, kFALSE);
541 }
0efb7042 542 }
543 }
544 delete list;
545 }
52e6652d 546
9e7b94f5 547 delete [] smatched;
52e6652d 548
549 AliInfo(Form("[%s] EventMixing %d/%d",GetName(),nEvents,nEvents));
61f275d1 550 timer.Stop(); timer.Print(); fflush(stdout);
52e6652d 551
0efb7042 552 /*
52e6652d 553 OLD
0efb7042 554 ifill = 0;
555 for (iloop = 1; iloop < nEvents; iloop++) {
556 imix = ievt + iloop;
557 // restart from beginning if reached last event
558 if (imix >= nEvents) imix -= nEvents;
559 // avoid to mix an event with itself
560 if (imix == ievt) continue;
561 // skip all events already mixed enough times
562 if (fNMixed[ievt] >= fNMix) break;
563 if (fNMixed[imix] >= fNMix) continue;
564 fEvBuffer->GetEntry(imix);
565 // skip if events are not matched
566 if (TMath::Abs(evMain.Vz() - fMiniEvent->Vz() ) > fMaxDiffVz ) continue;
567 if (TMath::Abs(evMain.Mult() - fMiniEvent->Mult() ) > fMaxDiffMult ) continue;
568 if (TMath::Abs(evMain.Angle() - fMiniEvent->Angle()) > fMaxDiffAngle) continue;
569 // found a match: increment counter for both events
570 AliDebugClass(1, Form("Event %d, def '%s': event mixing (%d with %d)", ievt, def->GetName(), ievt, imix));
571 fNMixed[ievt]++;
572 fNMixed[imix]++;
573 // process mixing
574 ifill += def->FillPair(&evMain, fMiniEvent, &fValues);
575 // stop if mixed enough times
576 if (fNMixed[ievt] >= fNMix) break;
577 }
52e6652d 578 break;
6aaeb33c 579 // print number of mixings done with each event
45aa62b9 580 for (ievt = 0; ievt < nEvents; ievt++) {
581 AliDebugClass(2, Form("Event %6d: mixed %2d times", ievt, fNMixed[ievt]));
582 }
0efb7042 583 */
52e6652d 584
6aaeb33c 585 // post computed data
03d23846 586 PostData(1, fOutput);
587}
588
589//__________________________________________________________________________________________________
590void AliRsnMiniAnalysisTask::Terminate(Option_t *)
591{
592//
593// Draw result to screen, or perform fitting, normalizations
594// Called once at the end of the query
595//
596
52e6652d 597 fOutput = dynamic_cast<TList *>(GetOutputData(1));
598 if (!fOutput) {
599 AliError("Could not retrieve TList fOutput");
600 return;
03d23846 601 }
602}
603
604//__________________________________________________________________________________________________
605Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent()
606{
607//
c80938f5 608// This method checks if current event is OK for analysis.
609// In case it is, the pointers of the local AliRsnEvent data member
610// will point to it, in order to allow cut checking, otherwise the
611// function exits with a failure message.
612// ---
613// ESD events must pass the physics selection, AOD are supposed to do.
614// ---
615// While checking the event, a histogram is filled to count the number
616// of CINT1B, V0AND and CANDLE events, which are needed for normalization
617// ---
618// Return values can be:
03d23846 619// -- 'E' if the event is accepted and is ESD
620// -- 'A' if the event is accepted and is AOD
621// -- 0 if the event is not accepted
622//
623
624 // string to sum messages
625 TString msg("");
52e6652d 626
03d23846 627 // check input type
c80938f5 628 // exit points are provided in all cases an event is bad
629 // if this block is passed, an event can be rejected only
630 // if it does not pass the set of event cuts defined in the task
03d23846 631 Char_t output = 0;
c80938f5 632 Bool_t isSelected;
633 if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
634 // type ESD
03d23846 635 output = 'E';
c80938f5 636 // ESD specific check: Physics Selection
637 // --> if this is failed, the event is rejected
4fb0dfa3 638 isSelected = (((AliInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
639
c80938f5 640 if (!isSelected) {
45aa62b9 641 AliDebugClass(2, "Event does not pass physics selections");
c80938f5 642 fRsnEvent.SetRef(0x0);
643 fRsnEvent.SetRefMC(0x0);
644 return 0;
645 }
646 // set reference to input
647 fRsnEvent.SetRef(fInputEvent);
648 // add MC if requested and available
649 if (fUseMC) {
52e6652d 650 if (fMCEvent)
c80938f5 651 fRsnEvent.SetRefMC(fMCEvent);
652 else {
653 AliWarning("MC event requested but not available");
654 fRsnEvent.SetRefMC(0x0);
655 }
656 }
657 } else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
658 // type AOD
03d23846 659 output = 'A';
c80938f5 660 // set reference to input
661 fRsnEvent.SetRef(fInputEvent);
662 // add MC if requested and available (it is in the same object)
663 if (fUseMC) {
664 fRsnEvent.SetRefMC(fInputEvent);
665 if (!fRsnEvent.GetAODList()) {
666 AliWarning("MC event requested but not available");
667 fRsnEvent.SetRefMC(0x0);
668 }
669 }
670 } else {
03d23846 671 AliError(Form("Bad input event class: %s", fInputEvent->ClassName()));
c80938f5 672 // reset pointers in local AliRsnEvent object
673 fRsnEvent.SetRef(0x0);
674 fRsnEvent.SetRefMC(0x0);
03d23846 675 return 0;
676 }
52e6652d 677
c80938f5 678 // fill counter of accepted events
679 fHEventStat->Fill(0.1);
52e6652d 680
03d23846 681 // check if it is V0AND
c80938f5 682 // --> uses a cast to AliESDEvent even if the input is an AliAODEvent
52e6652d 683 Bool_t v0A = fTriggerAna->IsOfflineTriggerFired((AliESDEvent *)fInputEvent, AliTriggerAnalysis::kV0A);
684 Bool_t v0C = fTriggerAna->IsOfflineTriggerFired((AliESDEvent *)fInputEvent, AliTriggerAnalysis::kV0C);
03d23846 685 if (v0A && v0C) {
686 msg += " -- VOAND = YES";
687 fHEventStat->Fill(1.1);
688 } else {
689 msg += " -- VOAND = NO ";
690 }
52e6652d 691
03d23846 692 // check candle
c80938f5 693 // --> requires at least one good quality track with Pt > 0.5 and |eta| <= 0.8
694 Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
03d23846 695 Bool_t candle = kFALSE;
52e6652d 696 for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
697 AliVTrack *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
698 AliESDtrack *esdt = dynamic_cast<AliESDtrack *>(track);
699 AliAODTrack *aodt = dynamic_cast<AliAODTrack *>(track);
03d23846 700 if (track->Pt() < 0.5) continue;
701 if(TMath::Abs(track->Eta()) > 0.8) continue;
c80938f5 702 if (esdt) if (!fESDtrackCuts->AcceptTrack(esdt)) continue;
703 if (aodt) if (!aodt->TestFilterBit(5)) continue;
03d23846 704 candle = kTRUE;
705 break;
706 }
707 if (candle) {
52e6652d 708 msg += " -- CANDLE = YES";
03d23846 709 fHEventStat->Fill(2.1);
710 } else {
52e6652d 711 msg += " -- CANDLE = NO ";
03d23846 712 }
52e6652d 713
c80938f5 714 // if event cuts are defined, they are checked here
715 // final decision on the event depends on this
716 isSelected = kTRUE;
03d23846 717 if (fEventCuts) {
718 if (!fEventCuts->IsSelected(&fRsnEvent)) {
719 msg += " -- Local cuts = REJECTED";
c80938f5 720 isSelected = kFALSE;
03d23846 721 } else {
722 msg += " -- Local cuts = ACCEPTED";
c80938f5 723 isSelected = kTRUE;
03d23846 724 }
725 } else {
726 msg += " -- Local cuts = NONE";
c80938f5 727 isSelected = kTRUE;
03d23846 728 }
52e6652d 729
03d23846 730 // if the above exit point is not taken, the event is accepted
9e7b94f5 731 AliDebugClass(2, Form("Stats: %s", msg.Data()));
c80938f5 732 if (isSelected) {
733 fHEventStat->Fill(3.1);
955ca046 734 Double_t multi = ComputeCentrality((output == 'E'));
735 fHAEventsVsMulti->Fill(multi);
5654276f 736 if(fHAEventVz) fHAEventVz->Fill(multi,fInputEvent->GetPrimaryVertex()->GetZ());
737 if(fHAEventMultiCent) fHAEventMultiCent->Fill(multi,ComputeMultiplicity(output == 'E',fHAEventMultiCent->GetYaxis()->GetTitle()));
738 if(fHAEventPlane) fHAEventPlane->Fill(multi,ComputeAngle());
c80938f5 739 return output;
740 } else {
741 return 0;
742 }
03d23846 743}
744
6aaeb33c 745//__________________________________________________________________________________________________
45aa62b9 746void AliRsnMiniAnalysisTask::FillMiniEvent(Char_t evType)
6aaeb33c 747{
748//
749// Refresh cursor mini-event data member to fill with current event.
750// Returns the total number of tracks selected.
751//
752
753 // assign event-related values
45aa62b9 754 if (fMiniEvent) delete fMiniEvent;
6aaeb33c 755 fMiniEvent = new AliRsnMiniEvent;
756 fMiniEvent->Vz() = fInputEvent->GetPrimaryVertex()->GetZ();
757 fMiniEvent->Angle() = ComputeAngle();
758 fMiniEvent->Mult() = ComputeCentrality((evType == 'E'));
45aa62b9 759 AliDebugClass(2, Form("Event %d: type = %c -- vz = %f -- mult = %f -- angle = %f", fEvNum, evType, fMiniEvent->Vz(), fMiniEvent->Mult(), fMiniEvent->Angle()));
52e6652d 760
6aaeb33c 761 // loop on daughters and assign track-related values
762 Int_t ic, ncuts = fTrackCuts.GetEntries();
763 Int_t ip, npart = fRsnEvent.GetAbsoluteSum();
45aa62b9 764 Int_t npos = 0, nneg = 0, nneu = 0;
6aaeb33c 765 AliRsnDaughter cursor;
766 AliRsnMiniParticle miniParticle;
767 for (ip = 0; ip < npart; ip++) {
768 // point cursor to next particle
769 fRsnEvent.SetDaughter(cursor, ip);
770 // copy momentum and MC info if present
771 miniParticle.CopyDaughter(&cursor);
45aa62b9 772 miniParticle.Index() = ip;
6aaeb33c 773 // switch on the bits corresponding to passed cuts
774 for (ic = 0; ic < ncuts; ic++) {
52e6652d 775 AliRsnCutSet *cuts = (AliRsnCutSet *)fTrackCuts[ic];
6aaeb33c 776 if (cuts->IsSelected(&cursor)) miniParticle.SetCutBit(ic);
777 }
778 // if a track passes at least one track cut, it is added to the pool
45aa62b9 779 if (miniParticle.CutBits()) {
780 fMiniEvent->AddParticle(miniParticle);
781 if (miniParticle.Charge() == '+') npos++;
782 else if (miniParticle.Charge() == '-') nneg++;
783 else nneu++;
784 }
6aaeb33c 785 }
52e6652d 786
6aaeb33c 787 // get number of accepted tracks
45aa62b9 788 AliDebugClass(1, Form("Event %6d: total = %5d, accepted = %4d (pos %4d, neg %4d, neu %4d)", fEvNum, npart, (Int_t)fMiniEvent->Particles().GetEntriesFast(), npos, nneg, nneu));
6aaeb33c 789}
790
3dced4f8 791//__________________________________________________________________________________________________
792Double_t AliRsnMiniAnalysisTask::ComputeAngle()
793{
794//
795// Get the plane angle
796//
797
cd273d20 798 AliEventplane *plane = 0x0;
52e6652d 799
cd273d20 800 if (fInputEvent->InheritsFrom(AliESDEvent::Class()))
801 plane = fInputEvent->GetEventplane();
802 else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
52e6652d 803 AliAODEvent *aodEvent = (AliAODEvent *)fInputEvent;
cd273d20 804 plane = aodEvent->GetHeader()->GetEventplaneP();
805 }
806
52e6652d 807 if (plane)
3dced4f8 808 return plane->GetEventplane("Q");
809 else {
810 AliWarning("No event plane defined");
811 return 1E20;
812 }
813}
814
03d23846 815//__________________________________________________________________________________________________
816Double_t AliRsnMiniAnalysisTask::ComputeCentrality(Bool_t isESD)
817{
818//
819// Computes event centrality/multiplicity according to the criterion defined
820// by two elements: (1) choice between multiplicity and centrality and
821// (2) the string defining what criterion must be used for specific computation.
822//
823
824 if (fUseCentrality) {
3da8cef7 825
826 if ((!fUseMC) && (!isESD) && (fUseAOD049CentralityPatch)) {
827 return ApplyCentralityPatchAOD049();
828 } else {
829 AliCentrality *centrality = fInputEvent->GetCentrality();
830 if (!centrality) {
831 AliError("Cannot compute centrality!");
832 return -1.0;
833 }
834 return centrality->GetCentralityPercentile(fCentralityType.Data());
835 }
03d23846 836 } else {
03d23846 837 if (!fCentralityType.CompareTo("TRACKS"))
838 return fInputEvent->GetNumberOfTracks();
839 else if (!fCentralityType.CompareTo("QUALITY"))
c80938f5 840 if (isESD)
52e6652d 841 return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent *)fInputEvent, kTRUE);
c80938f5 842 else {
843 Double_t count = 0.;
844 Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
52e6652d 845 for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
846 AliVTrack *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
847 AliAODTrack *aodt = dynamic_cast<AliAODTrack *>(track);
c80938f5 848 if (!aodt) continue;
849 if (!aodt->TestFilterBit(5)) continue;
850 count++;
851 }
852 return count;
853 }
03d23846 854 else if (!fCentralityType.CompareTo("TRACKLETS")) {
c80938f5 855 if (isESD) {
52e6652d 856 const AliMultiplicity *mult = ((AliESDEvent *)fInputEvent)->GetMultiplicity();
c80938f5 857 Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
858 for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
859 return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ());
860 } else {
861 AliWarning("Cannot compute multiplicity with SPD tracklets from AOD");
862 return 1E20;
863 }
03d23846 864 } else {
865 AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", fCentralityType.Data()));
866 return -1.0;
867 }
868 }
869}
870
5654276f 871//__________________________________________________________________________________________________
872Double_t AliRsnMiniAnalysisTask::ComputeMultiplicity(Bool_t isESD,TString type)
873{
874//
875// Computes event multiplicity according to the string defining
876// what criterion must be used for specific computation.
877//
878
e6952ec7 879 type.ToUpper();
880
881 if (!type.CompareTo("TRACKS"))
882 return fInputEvent->GetNumberOfTracks();
883 else if (!type.CompareTo("QUALITY"))
884 if (isESD)
885 return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent *)fInputEvent, kTRUE);
886 else {
887 Double_t count = 0.;
888 Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
889 for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
890 AliVTrack *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
891 AliAODTrack *aodt = dynamic_cast<AliAODTrack *>(track);
892 if (!aodt) continue;
893 if (!aodt->TestFilterBit(5)) continue;
894 count++;
895 }
896 return count;
5654276f 897 }
e6952ec7 898 else if (!type.CompareTo("TRACKLETS")) {
899 if (isESD) {
900 const AliMultiplicity *mult = ((AliESDEvent *)fInputEvent)->GetMultiplicity();
901 Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
902 for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
903 return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ());
904 } else {
905 AliWarning("Cannot compute multiplicity with SPD tracklets from AOD");
906 return 1E20;
907 }
908 } else {
909 AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", type.Data()));
910 return -1.0;
911 }
5654276f 912
e6952ec7 913 return 1E20;
5654276f 914}
915
03d23846 916//__________________________________________________________________________________________________
917void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
918{
919//
920// Fills the histograms with true mother (ESD version)
921//
922
923 Bool_t okMatch;
924 Int_t id, ndef = fHistograms.GetEntries();
925 Int_t ip, label1, label2, npart = fMCEvent->GetNumberOfTracks();
926 static AliRsnMiniPair miniPair;
927 AliMCParticle *daughter1, *daughter2;
928 TLorentzVector p1, p2;
929 AliRsnMiniOutput *def = 0x0;
52e6652d 930
03d23846 931 for (id = 0; id < ndef; id++) {
52e6652d 932 def = (AliRsnMiniOutput *)fHistograms[id];
03d23846 933 if (!def) continue;
934 if (!def->IsMother()) continue;
935 for (ip = 0; ip < npart; ip++) {
52e6652d 936 AliMCParticle *part = (AliMCParticle *)fMCEvent->GetTrack(ip);
e6952ec7 937 //get mother pdg code
64129b35 938 if (part->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
03d23846 939 // check that daughters match expected species
088ca370 940 if (part->Particle()->GetNDaughters() < 2) continue;
941 if (fMaxNDaughters > 0 && part->Particle()->GetNDaughters() > fMaxNDaughters) continue;
03d23846 942 label1 = part->Particle()->GetDaughter(0);
943 label2 = part->Particle()->GetDaughter(1);
52e6652d 944 daughter1 = (AliMCParticle *)fMCEvent->GetTrack(label1);
945 daughter2 = (AliMCParticle *)fMCEvent->GetTrack(label2);
03d23846 946 okMatch = kFALSE;
947 if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(1)) {
948 okMatch = kTRUE;
949 p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
950 p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
951 } else if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(0)) {
952 okMatch = kTRUE;
953 p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
954 p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
955 }
956 if (!okMatch) continue;
088ca370 957 if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
958 (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
959 (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
31dbef4e 960 if(fCheckFeedDown){
961 Int_t pdgGranma = 0;
962 Int_t mother = 0;
963 mother = part->GetMother();
964 Int_t istep = 0;
965 Int_t abspdgGranma =0;
966 Bool_t isFromB=kFALSE;
967 Bool_t isQuarkFound=kFALSE;
968 while (mother >=0 ){
969 istep++;
970 AliDebug(2,Form("mother at step %d = %d", istep, mother));
971 AliMCParticle* mcGranma = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(mother));
972 if (mcGranma){
973 pdgGranma = mcGranma->PdgCode();
974 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
975 abspdgGranma = TMath::Abs(pdgGranma);
976 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
977 isFromB=kTRUE;
978 }
979 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
980 mother = mcGranma->GetMother();
981 }else{
982 AliError("Failed casting the mother particle!");
983 break;
984 }
985 }
986 if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
987 if(isFromB){
988 if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
989 }
990 else{
991 if (fKeepDfromBOnly) pdgGranma = -999;
992
993 if (pdgGranma == -99999){
994 AliDebug(2,"This particle does not have a quark in his genealogy\n");
995 continue;
996 }
997 if (pdgGranma == -9999){
998 AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");
999 continue;
1000 }
1001
1002 if (pdgGranma == -999){
1003 AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");
1004 continue;
1005 }
1006
1007 }
1008 }
03d23846 1009 // assign momenta to computation object
1010 miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
1011 miniPair.FillRef(def->GetMotherMass());
1012 // do computations
45aa62b9 1013 def->FillMother(&miniPair, miniEvent, &fValues);
03d23846 1014 }
1015 }
1016}
1017
1018//__________________________________________________________________________________________________
1019void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
1020{
1021//
1022// Fills the histograms with true mother (AOD version)
1023//
1024
1025 Bool_t okMatch;
1026 TClonesArray *list = fRsnEvent.GetAODList();
1027 Int_t id, ndef = fHistograms.GetEntries();
1028 Int_t ip, label1, label2, npart = list->GetEntries();
1029 static AliRsnMiniPair miniPair;
1030 AliAODMCParticle *daughter1, *daughter2;
1031 TLorentzVector p1, p2;
1032 AliRsnMiniOutput *def = 0x0;
52e6652d 1033
03d23846 1034 for (id = 0; id < ndef; id++) {
52e6652d 1035 def = (AliRsnMiniOutput *)fHistograms[id];
03d23846 1036 if (!def) continue;
1037 if (!def->IsMother()) continue;
1038 for (ip = 0; ip < npart; ip++) {
52e6652d 1039 AliAODMCParticle *part = (AliAODMCParticle *)list->At(ip);
64129b35 1040 if (part->GetPdgCode() != def->GetMotherPDG()) continue;
03d23846 1041 // check that daughters match expected species
1042 if (part->GetNDaughters() < 2) continue;
088ca370 1043 if (fMaxNDaughters > 0 && part->GetNDaughters() > fMaxNDaughters) continue;
03d23846 1044 label1 = part->GetDaughter(0);
1045 label2 = part->GetDaughter(1);
52e6652d 1046 daughter1 = (AliAODMCParticle *)list->At(label1);
1047 daughter2 = (AliAODMCParticle *)list->At(label2);
03d23846 1048 okMatch = kFALSE;
1049 if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(1)) {
1050 okMatch = kTRUE;
1051 p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
1052 p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
1053 } else if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(0)) {
1054 okMatch = kTRUE;
1055 p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
1056 p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
1057 }
1058 if (!okMatch) continue;
088ca370 1059 if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
1060 (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
31dbef4e 1061 (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
1062 if(fCheckFeedDown){
1063 Int_t pdgGranma = 0;
1064 Int_t mother = 0;
1065 mother = part->GetMother();
1066 Int_t istep = 0;
1067 Int_t abspdgGranma =0;
1068 Bool_t isFromB=kFALSE;
1069 Bool_t isQuarkFound=kFALSE;
1070 while (mother >=0 ){
1071 istep++;
1072 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1073 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(list->At(mother));
1074 if (mcGranma){
1075 pdgGranma = mcGranma->GetPdgCode();
1076 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1077 abspdgGranma = TMath::Abs(pdgGranma);
1078 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1079 isFromB=kTRUE;
1080 }
1081 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1082 mother = mcGranma->GetMother();
1083 }else{
1084 AliError("Failed casting the mother particle!");
1085 break;
1086 }
1087 }
1088 if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
1089 if(isFromB){
1090 if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
1091 }
1092 else{
1093 if (fKeepDfromBOnly) pdgGranma = -999;
1094 }
1095
1096 if (pdgGranma == -99999){
1097 AliDebug(2,"This particle does not have a quark in his genealogy\n");
1098 continue;
1099 }
1100 if (pdgGranma == -9999){
1101 AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");
1102 continue;
1103 }
1104
1105 if (pdgGranma == -999){
1106 AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");
1107 continue;
1108 }
1109 }
088ca370 1110 // assign momenta to computation object
03d23846 1111 miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
1112 miniPair.FillRef(def->GetMotherMass());
1113 // do computations
45aa62b9 1114 def->FillMother(&miniPair, miniEvent, &fValues);
03d23846 1115 }
1116 }
1117}
31dbef4e 1118//___________________________________________________________
1119void AliRsnMiniAnalysisTask::SetDselection(UShort_t originDselection)
1120{
1121 // setting the way the D0 will be selected
1122 // 0 --> only from c quarks
1123 // 1 --> only from b quarks
1124 // 2 --> from both c quarks and b quarks
1125
1126 fOriginDselection = originDselection;
1127
1128 if (fOriginDselection == 0) {
1129 fKeepDfromB = kFALSE;
1130 fKeepDfromBOnly = kFALSE;
1131 }
1132
1133 if (fOriginDselection == 1) {
1134 fKeepDfromB = kTRUE;
1135 fKeepDfromBOnly = kTRUE;
1136 }
1137
1138 if (fOriginDselection == 2) {
1139 fKeepDfromB = kTRUE;
1140 fKeepDfromBOnly = kFALSE;
1141 }
1142
1143 return;
1144}
caef2e3c 1145//__________________________________________________________________________________________________
1146Bool_t AliRsnMiniAnalysisTask::EventsMatch(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2)
1147{
1148//
1149// Check if two events are compatible.
1150// If the mixing is continuous, this is true if differences in vz, mult and angle are smaller than
52e6652d 1151// the specified values.
caef2e3c 1152// If the mixing is binned, this is true if the events are in the same bin.
1153//
1154
1155 if (!event1 || !event2) return kFALSE;
1156 Int_t ivz1, ivz2, imult1, imult2, iangle1, iangle2;
9e7b94f5 1157 Double_t dv, dm, da;
52e6652d 1158
caef2e3c 1159 if (fContinuousMix) {
9e7b94f5 1160 dv = TMath::Abs(event1->Vz() - event2->Vz() );
1161 dm = TMath::Abs(event1->Mult() - event2->Mult() );
1162 da = TMath::Abs(event1->Angle() - event2->Angle());
1163 if (dv > fMaxDiffVz) {
1164 //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Vz = %f", event1->ID(), event2->ID(), dv));
1165 return kFALSE;
1166 }
1167 if (dm > fMaxDiffMult ) {
1168 //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Mult = %f", event1->ID(), event2->ID(), dm));
1169 return kFALSE;
1170 }
1171 if (da > fMaxDiffAngle) {
1172 //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Angle = %f", event1->ID(), event2->ID(), da));
1173 return kFALSE;
1174 }
caef2e3c 1175 return kTRUE;
1176 } else {
1177 ivz1 = (Int_t)(event1->Vz() / fMaxDiffVz);
1178 ivz2 = (Int_t)(event2->Vz() / fMaxDiffVz);
1179 imult1 = (Int_t)(event1->Mult() / fMaxDiffMult);
1180 imult2 = (Int_t)(event2->Mult() / fMaxDiffMult);
1181 iangle1 = (Int_t)(event1->Angle() / fMaxDiffAngle);
1182 iangle2 = (Int_t)(event2->Angle() / fMaxDiffAngle);
1183 if (ivz1 != ivz2) return kFALSE;
1184 if (imult1 != imult2) return kFALSE;
1185 if (iangle1 != iangle2) return kFALSE;
1186 return kTRUE;
1187 }
1188}
1189
52e6652d 1190
1cdc3671 1191//---------------------------------------------------------------------
1192Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchAOD049()
1193{
3da8cef7 1194 //
1195 //Apply centrality patch for AOD049 outliers
1196 //
1197 if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
1198 AliWarning("Requested patch for AOD049 for ESD. ");
1cdc3671 1199 return -999.0;
3da8cef7 1200 }
1201
1202 if (fCentralityType!="V0M") {
1203 AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
1cdc3671 1204 return -999.0;
3da8cef7 1205 }
1206
1207 AliCentrality *centrality = fInputEvent->GetCentrality();
1208 if (!centrality) {
1209 AliWarning("Cannot get centrality from AOD event.");
1210 return -999.0;
1211 }
1212
1213 Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
1214 /*
1215 Bool_t isSelRun = kFALSE;
1216 Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
1217 if(cent<0){
1218 Int_t quality = centrality->GetQuality();
1219 if(quality<=1){
1220 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1221 } else {
1222 Int_t runnum=aodEvent->GetRunNumber();
1223 for(Int_t ir=0;ir<5;ir++){
1224 if(runnum==selRun[ir]){
1225 isSelRun=kTRUE;
1226 break;
1227 }
1228 }
1229 if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1230 }
1231 }
1232 */
1233
1234 if(cent>=0.0) {
1235 Float_t v0 = 0.0;
1236 AliAODEvent *aodEvent = (AliAODEvent *)fInputEvent;
1237 AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
1238 v0+=aodV0->GetMTotV0A();
1239 v0+=aodV0->GetMTotV0C();
1240 if ( (cent==0) && (v0<19500) ) {
1241 AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
1242 return -999.0;
1243 }
1244 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1245 Float_t val = 1.30552 + 0.147931 * v0;
1246
1247 Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
1248 120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
1249 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
1250 68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
1251 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
1252 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
1253 26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
1254 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
1255 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
1256 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
1257 };
1258
1259 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] ) {
1260 AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
1261 return -999.0;
1262 }
1263 } else {
1264 //force it to be -999. whatever the negative value was
1265 cent = -999.;
1266 }
1267 return cent;
1cdc3671 1268}
5654276f 1269
1270//----------------------------------------------------------------------------------
e6952ec7 1271void AliRsnMiniAnalysisTask::SetEventQAHist(TString type,TH2F *histo)
5654276f 1272{
e6952ec7 1273 if(!histo) {
1274 AliWarning(Form("event QA histogram pointer not defined for slot %s",type.Data()));
1275 return;
1276 }
1277
1278 type.ToLower();
1279
1280 if(!type.CompareTo("vz")) fHAEventVz = histo;
1281 else if(!type.CompareTo("multicent")) {
1282 TString mtype(histo->GetYaxis()->GetTitle());
1283 mtype.ToUpper();
1284 if(mtype.CompareTo("QUALITY") && mtype.CompareTo("TRACKS") && mtype.CompareTo("TRACKLETS")) {
1285 AliWarning(Form("multiplicity vs. centrality histogram y-axis %s unknown, setting to TRACKS",mtype.Data()));
1286 histo->GetYaxis()->SetTitle("TRACKS");
1287 }
1288 fHAEventMultiCent = histo;
1289 }
1290 else if(!type.CompareTo("eventplane")) fHAEventPlane = histo;
1291 else AliWarning(Form("event QA histogram slot %s undefined",type.Data()));
1292
1293 return;
5654276f 1294}
4fb0dfa3 1295
1296//----------------------------------------------------------------------------------
1297Int_t AliRsnMiniAnalysisTask::CreateValue(AliRsnMiniValue::EType type, Bool_t useMC)
1298{
1299//
1300// Create a new value in the task,
1301// and returns its ID, which is needed for setting up histograms.
1302// If that value was already initialized, returns its ID and does not recreate it.
1303//
1304
1305 Int_t valID = ValueID(type, useMC);
1306 if (valID >= 0 && valID < fValues.GetEntries()) {
1307 AliInfo(Form("Value '%s' is already created in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
1308 } else {
1309 valID = fValues.GetEntries();
1310 AliInfo(Form("Creating value '%s' in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
1311 new (fValues[valID]) AliRsnMiniValue(type, useMC);
1312 }
1313
1314 return valID;
1315}
1316
1317//----------------------------------------------------------------------------------
1318Int_t AliRsnMiniAnalysisTask::ValueID(AliRsnMiniValue::EType type, Bool_t useMC)
1319{
1320//
1321// Searches if a value computation is initialized
1322//
1323
1324 const char *name = AliRsnMiniValue::ValueName(type, useMC);
1325 TObject *obj = fValues.FindObject(name);
1326 if (obj)
1327 return fValues.IndexOf(obj);
1328 else
1329 return -1;
1330}
1331
1332//----------------------------------------------------------------------------------
1333AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, AliRsnMiniOutput::EOutputType type, AliRsnMiniOutput::EComputation src)
1334{
1335//
1336// Create a new histogram definition in the task,
1337// which is then returned to the user for its configuration
1338//
1339
1340 Int_t n = fHistograms.GetEntries();
1341 AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, type, src);
1342
1343 return newDef;
1344}
1345
1346//----------------------------------------------------------------------------------
1347AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, const char *outType, const char *compType)
1348{
1349//
1350// Create a new histogram definition in the task,
1351// which is then returned to the user for its configuration
1352//
1353
1354 Int_t n = fHistograms.GetEntries();
1355 AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, outType, compType);
1356
1357 return newDef;
1358}