2 // Analysis task for 'mini' sub-package
3 // Contains all definitions needed for running an analysis:
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.
13 #include <Riostream.h>
18 #include <TStopwatch.h>
22 #include "AliEventplane.h"
23 #include "AliMultiplicity.h"
24 #include "AliTriggerAnalysis.h"
25 #include "AliAnalysisManager.h"
26 #include "AliInputEventHandler.h"
28 #include "AliESDtrackCuts.h"
29 #include "AliESDUtils.h"
31 #include "AliAODEvent.h"
32 #include "AliAODMCParticle.h"
34 #include "AliRsnCutSet.h"
35 #include "AliRsnMiniPair.h"
36 #include "AliRsnMiniEvent.h"
37 #include "AliRsnMiniParticle.h"
39 #include "AliRsnMiniAnalysisTask.h"
42 ClassImp(AliRsnMiniAnalysisTask)
44 //__________________________________________________________________________________________________
45 AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() :
50 fUseCentrality(kFALSE),
51 fCentralityType("QUALITY"),
52 fUseAOD049CentralityPatch(kFALSE),
53 fUseCentralityPatchPbPb2011(0),
54 fContinuousMix(kTRUE),
60 fHistograms("AliRsnMiniOutput", 0),
61 fValues("AliRsnMiniValue", 0),
63 fHAEventsVsMulti(0x0),
64 fHAEventsVsTracklets(0x0),
66 fHAEventMultiCent(0x0),
79 fCheckFeedDown(kFALSE),
80 fOriginDselection(kFALSE),
82 fKeepDfromBOnly(kFALSE),
83 fRejectIfNoQuark(kFALSE),
84 fMotherAcceptanceCutMinPt(0.0),
85 fMotherAcceptanceCutMaxEta(0.9),
86 fKeepMotherInAcceptance(kFALSE)
89 // Dummy constructor ALWAYS needed for I/O.
93 //__________________________________________________________________________________________________
94 AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name, Bool_t useMC) :
95 AliAnalysisTaskSE(name),
98 fTriggerMask(AliVEvent::kMB),
99 fUseCentrality(kFALSE),
100 fCentralityType("QUALITY"),
101 fUseAOD049CentralityPatch(kFALSE),
102 fUseCentralityPatchPbPb2011(0),
103 fContinuousMix(kTRUE),
109 fHistograms("AliRsnMiniOutput", 0),
110 fValues("AliRsnMiniValue", 0),
112 fHAEventsVsMulti(0x0),
113 fHAEventsVsTracklets(0x0),
115 fHAEventMultiCent(0x0),
125 fMixPrintRefresh(-1),
128 fCheckFeedDown(kFALSE),
129 fOriginDselection(kFALSE),
131 fKeepDfromBOnly(kFALSE),
132 fRejectIfNoQuark(kFALSE),
133 fMotherAcceptanceCutMinPt(0.0),
134 fMotherAcceptanceCutMaxEta(0.9),
135 fKeepMotherInAcceptance(kFALSE)
138 // Default constructor.
139 // Define input and output slots here (never in the dummy constructor)
140 // Input slot #0 works with a TChain - it is connected to the default input container
141 // Output slot #1 writes into a TH1 container
144 DefineOutput(1, TList::Class());
147 //__________________________________________________________________________________________________
148 AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask ©) :
149 AliAnalysisTaskSE(copy),
152 fTriggerMask(copy.fTriggerMask),
153 fUseCentrality(copy.fUseCentrality),
154 fCentralityType(copy.fCentralityType),
155 fUseAOD049CentralityPatch(copy.fUseAOD049CentralityPatch),
156 fUseCentralityPatchPbPb2011(copy.fUseCentralityPatchPbPb2011),
157 fContinuousMix(copy.fContinuousMix),
159 fMaxDiffMult(copy.fMaxDiffMult),
160 fMaxDiffVz(copy.fMaxDiffVz),
161 fMaxDiffAngle(copy.fMaxDiffAngle),
163 fHistograms(copy.fHistograms),
164 fValues(copy.fValues),
166 fHAEventsVsMulti(0x0),
167 fHAEventsVsTracklets(0x0),
169 fHAEventMultiCent(0x0),
171 fEventCuts(copy.fEventCuts),
172 fTrackCuts(copy.fTrackCuts),
175 fTriggerAna(copy.fTriggerAna),
176 fESDtrackCuts(copy.fESDtrackCuts),
178 fBigOutput(copy.fBigOutput),
179 fMixPrintRefresh(copy.fMixPrintRefresh),
180 fMaxNDaughters(copy.fMaxNDaughters),
181 fCheckP(copy.fCheckP),
182 fCheckFeedDown(copy.fCheckFeedDown),
183 fOriginDselection(copy.fOriginDselection),
184 fKeepDfromB(copy.fOriginDselection),
185 fKeepDfromBOnly(copy.fKeepDfromBOnly),
186 fRejectIfNoQuark(copy.fRejectIfNoQuark),
187 fMotherAcceptanceCutMinPt(copy.fMotherAcceptanceCutMinPt),
188 fMotherAcceptanceCutMaxEta(copy.fMotherAcceptanceCutMaxEta),
189 fKeepMotherInAcceptance(copy.fKeepMotherInAcceptance)
193 // Implemented as requested by C++ standards.
194 // Can be used in PROOF and by plugins.
198 //__________________________________________________________________________________________________
199 AliRsnMiniAnalysisTask &AliRsnMiniAnalysisTask::operator=(const AliRsnMiniAnalysisTask ©)
202 // Assignment operator.
203 // Implemented as requested by C++ standards.
204 // Can be used in PROOF and by plugins.
207 AliAnalysisTaskSE::operator=(copy);
210 fUseMC = copy.fUseMC;
211 fEvNum = copy.fEvNum;
212 fTriggerMask = copy.fTriggerMask;
213 fUseCentrality = copy.fUseCentrality;
214 fCentralityType = copy.fCentralityType;
215 fUseAOD049CentralityPatch = copy.fUseAOD049CentralityPatch;
216 fUseCentralityPatchPbPb2011 = copy.fUseCentralityPatchPbPb2011;
217 fContinuousMix = copy.fContinuousMix;
219 fMaxDiffMult = copy.fMaxDiffMult;
220 fMaxDiffVz = copy.fMaxDiffVz;
221 fMaxDiffAngle = copy.fMaxDiffAngle;
222 fHistograms = copy.fHistograms;
223 fValues = copy.fValues;
224 fHEventStat = copy.fHEventStat;
225 fHAEventsVsMulti = copy.fHAEventsVsMulti;
226 fHAEventsVsTracklets = copy.fHAEventsVsTracklets;
227 fHAEventVz = copy.fHAEventVz;
228 fHAEventMultiCent = copy.fHAEventMultiCent;
229 fHAEventPlane = copy.fHAEventPlane;
230 fEventCuts = copy.fEventCuts;
231 fTrackCuts = copy.fTrackCuts;
232 fTriggerAna = copy.fTriggerAna;
233 fESDtrackCuts = copy.fESDtrackCuts;
234 fBigOutput = copy.fBigOutput;
235 fMixPrintRefresh = copy.fMixPrintRefresh;
236 fMaxNDaughters = copy.fMaxNDaughters;
237 fCheckP = copy.fCheckP;
238 fCheckFeedDown = copy.fCheckFeedDown;
239 fOriginDselection = copy.fOriginDselection;
240 fKeepDfromB = copy.fOriginDselection;
241 fKeepDfromBOnly = copy.fKeepDfromBOnly;
242 fRejectIfNoQuark = copy.fRejectIfNoQuark;
243 fMotherAcceptanceCutMinPt = copy.fMotherAcceptanceCutMinPt;
244 fMotherAcceptanceCutMaxEta = copy.fMotherAcceptanceCutMaxEta;
245 fKeepMotherInAcceptance = copy.fKeepMotherInAcceptance;
249 //__________________________________________________________________________________________________
250 AliRsnMiniAnalysisTask::~AliRsnMiniAnalysisTask()
254 // Clean-up the output list, but not the histograms that are put inside
255 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
259 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
265 //__________________________________________________________________________________________________
266 Int_t AliRsnMiniAnalysisTask::AddTrackCuts(AliRsnCutSet *cuts)
269 // Add a new cut set for a new criterion for track selection.
270 // A user can add as many as he wants, and each one corresponds
271 // to one of the available bits in the AliRsnMiniParticle mask.
272 // The only check is the following: if a cut set with the same name
273 // as the argument is there, this is not added.
274 // Return value is the array position of this set.
277 TObject *obj = fTrackCuts.FindObject(cuts->GetName());
280 AliInfo(Form("A cut set named '%s' already exists", cuts->GetName()));
281 return fTrackCuts.IndexOf(obj);
283 fTrackCuts.AddLast(cuts);
284 return fTrackCuts.IndexOf(cuts);
288 //__________________________________________________________________________________________________
289 void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
292 // Initialization of outputs.
293 // This is called once per worker node.
300 AliInfo(Form("Selected event characterization: %s (%s)", (fUseCentrality ? "centrality" : "multiplicity"), fCentralityType.Data()));
302 // initialize trigger analysis
303 if (fTriggerAna) delete fTriggerAna;
304 fTriggerAna = new AliTriggerAnalysis;
306 // initialize ESD quality cuts
307 if (fESDtrackCuts) delete fESDtrackCuts;
308 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
310 // create list and set it as owner of its content (MANDATORY)
311 if (fBigOutput) OpenFile(1);
312 fOutput = new TList();
315 // initialize event statistics counter
316 fHEventStat = new TH1F("hEventStat", "Event statistics", 8, 0.0, 8.0);
317 fHEventStat->GetXaxis()->SetBinLabel(1, "CINT1B");
318 fHEventStat->GetXaxis()->SetBinLabel(2, "V0AND");
319 fHEventStat->GetXaxis()->SetBinLabel(3, "Candle");
320 fHEventStat->GetXaxis()->SetBinLabel(4, "Accepted");
321 fHEventStat->GetXaxis()->SetBinLabel(5, "Not Accepted - Total");
322 fHEventStat->GetXaxis()->SetBinLabel(6, "Not Accepted - No Track Vertex");
323 fHEventStat->GetXaxis()->SetBinLabel(7, "Not Accepted - Not Enough Contributors");
324 fHEventStat->GetXaxis()->SetBinLabel(8, "Not Accepted - No Vertex inside |z| < 10 cm");
326 fOutput->Add(fHEventStat);
329 fHAEventsVsMulti = new TH1F("hAEventsVsMulti", "Accepted events vs Centrality", 100, 0, 100.0);
331 fHAEventsVsMulti = new TH1F("hAEventsVsMulti", "Accepted events vs Multiplicity",1000, 0, 1000.0);
332 fOutput->Add(fHAEventsVsMulti);
334 fHAEventsVsTracklets = new TH1F("hAEventsVsTracklets", "Accepted events vs Tracklet Number",1000, 0, 1000.0);
335 fOutput->Add(fHAEventsVsTracklets);
337 if(fHAEventVz) fOutput->Add(fHAEventVz);
338 if(fHAEventMultiCent) fOutput->Add(fHAEventMultiCent);
339 if(fHAEventPlane) fOutput->Add(fHAEventPlane);
341 TIter next(&fTrackCuts);
343 while ((cs = (AliRsnCutSet *) next())) {
347 // create temporary tree for filtered events
348 if (fMiniEvent) delete fMiniEvent;
349 fEvBuffer = new TTree("EventBuffer", "Temporary buffer for mini events");
350 fEvBuffer->Branch("events", "AliRsnMiniEvent", &fMiniEvent);
352 // create one histogram per each stored definition (event histograms)
353 Int_t i, ndef = fHistograms.GetEntries();
354 AliRsnMiniOutput *def = 0x0;
355 for (i = 0; i < ndef; i++) {
356 def = (AliRsnMiniOutput *)fHistograms[i];
358 if (!def->Init(GetName(), fOutput)) {
359 AliError(Form("Def '%s': failed initialization", def->GetName()));
364 // post data for ALL output slots >0 here, to get at least an empty histogram
365 PostData(1, fOutput);
368 //__________________________________________________________________________________________________
369 void AliRsnMiniAnalysisTask::UserExec(Option_t *)
373 // In this case, it checks if the event is acceptable, and eventually
374 // creates the corresponding mini-event and stores it in the buffer.
375 // The real histogram filling is done at the end, in "FinishTaskOutput".
378 // increment event counter
381 // check current event
382 Char_t check = CheckCurrentEvent();
385 // setup PID response
386 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
387 AliInputEventHandler *inputHandler = (AliInputEventHandler *)man->GetInputEventHandler();
388 fRsnEvent.SetPIDResponse(inputHandler->GetPIDResponse());
390 // fill a mini-event from current
391 // and skip this event if no tracks were accepted
392 FillMiniEvent(check);
394 // fill MC based histograms on mothers,
395 // which do need the original event
397 if (fRsnEvent.IsESD() && fMCEvent)
398 FillTrueMotherESD(fMiniEvent);
399 else if (fRsnEvent.IsAOD() && fRsnEvent.GetAODList())
400 FillTrueMotherAOD(fMiniEvent);
403 // if the event is not empty, store it
404 if (fMiniEvent->IsEmpty()) {
405 AliDebugClass(2, Form("Rejecting empty event #%d", fEvNum));
407 Int_t id = fEvBuffer->GetEntries();
408 AliDebugClass(2, Form("Adding event #%d with ID = %d", fEvNum, id));
409 fMiniEvent->ID() = id;
413 // post data for computed stuff
414 PostData(1, fOutput);
417 //__________________________________________________________________________________________________
418 void AliRsnMiniAnalysisTask::FinishTaskOutput()
421 // This function is called at the end of the loop on available events,
422 // and then the buffer will be full with all the corresponding mini-events,
423 // each one containing all tracks selected by each of the available track cuts.
424 // Here a loop is done on each of these events, and both single-event and mixing are computed
427 // security code: reassign the buffer to the mini-event cursor
428 fEvBuffer->SetBranchAddress("events", &fMiniEvent);
431 Int_t ievt, nEvents = (Int_t)fEvBuffer->GetEntries();
432 Int_t idef, nDefs = fHistograms.GetEntries();
433 Int_t imix, iloop, ifill;
434 AliRsnMiniOutput *def = 0x0;
435 AliRsnMiniOutput::EComputation compType;
437 Int_t printNum = fMixPrintRefresh;
439 if (nEvents>1e5) printNum=nEvents/100;
440 else if (nEvents>1e4) printNum=nEvents/10;
444 // loop on events, and for each one fill all outputs
445 // using the appropriate procedure depending on its type
446 // only mother-related histograms are filled in UserExec,
447 // since they require direct access to MC event
449 for (ievt = 0; ievt < nEvents; ievt++) {
451 fEvBuffer->GetEntry(ievt);
452 if (printNum&&(ievt%printNum==0)) {
453 AliInfo(Form("[%s] Std.Event %d/%d",GetName(), ievt,nEvents));
454 timer.Stop(); timer.Print(); fflush(stdout); timer.Start(kFALSE);
457 for (idef = 0; idef < nDefs; idef++) {
458 def = (AliRsnMiniOutput *)fHistograms[idef];
460 compType = def->GetComputation();
461 // execute computation in the appropriate way
463 case AliRsnMiniOutput::kEventOnly:
464 //AliDebugClass(1, Form("Event %d, def '%s': event-value histogram filling", ievt, def->GetName()));
466 def->FillEvent(fMiniEvent, &fValues);
468 case AliRsnMiniOutput::kTruePair:
469 //AliDebugClass(1, Form("Event %d, def '%s': true-pair histogram filling", ievt, def->GetName()));
470 ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
472 case AliRsnMiniOutput::kTrackPair:
473 //AliDebugClass(1, Form("Event %d, def '%s': pair-value histogram filling", ievt, def->GetName()));
474 ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
476 case AliRsnMiniOutput::kTrackPairRotated1:
477 //AliDebugClass(1, Form("Event %d, def '%s': rotated (1) background histogram filling", ievt, def->GetName()));
478 ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
480 case AliRsnMiniOutput::kTrackPairRotated2:
481 //AliDebugClass(1, Form("Event %d, def '%s': rotated (2) background histogram filling", ievt, def->GetName()));
482 ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
485 // other kinds are processed elsewhere
487 AliDebugClass(2, Form("Computation = %d", (Int_t)compType));
490 AliDebugClass(1, Form("Event %6d: def = '%15s' -- fills = %5d", ievt, def->GetName(), ifill));
494 // if no mixing is required, stop here and post the output
496 AliDebugClass(2, "Stopping here, since no mixing is required");
497 PostData(1, fOutput);
501 // initialize mixing counter
502 Int_t nmatched[nEvents];
503 TString *smatched = new TString[nEvents];
504 for (ievt = 0; ievt < nEvents; ievt++) {
505 smatched[ievt] = "|";
510 AliInfo(Form("[%s] Std.Event %d/%d",GetName(), nEvents,nEvents));
511 timer.Stop(); timer.Print(); timer.Start(); fflush(stdout);
513 // search for good matchings
514 for (ievt = 0; ievt < nEvents; ievt++) {
515 if (printNum&&(ievt%printNum==0)) {
516 AliInfo(Form("[%s] EventMixing searching %d/%d",GetName(),ievt,nEvents));
517 timer.Stop(); timer.Print(); timer.Start(kFALSE); fflush(stdout);
519 if (nmatched[ievt] >= fNMix) continue;
520 fEvBuffer->GetEntry(ievt);
521 AliRsnMiniEvent evMain(*fMiniEvent);
522 for (iloop = 1; iloop < nEvents; iloop++) {
524 if (imix >= nEvents) imix -= nEvents;
525 if (imix == ievt) continue;
527 fEvBuffer->GetEntry(imix);
528 // skip if events are not matched
529 if (!EventsMatch(&evMain, fMiniEvent)) continue;
530 // check that the array of good matches for mixed does not already contain main event
531 if (smatched[imix].Contains(Form("|%d|", ievt))) continue;
532 // check that the found good events has not enough matches already
533 if (nmatched[imix] >= fNMix) continue;
534 // add new mixing candidate
535 smatched[ievt].Append(Form("%d|", imix));
538 if (nmatched[ievt] >= fNMix) break;
540 AliDebugClass(1, Form("Matches for event %5d = %d [%s] (missing are declared above)", evMain.ID(), nmatched[ievt], smatched[ievt].Data()));
543 AliInfo(Form("[%s] EventMixing searching %d/%d",GetName(),nEvents,nEvents));
544 timer.Stop(); timer.Print(); fflush(stdout); timer.Start();
547 TObjArray *list = 0x0;
548 TObjString *os = 0x0;
549 for (ievt = 0; ievt < nEvents; ievt++) {
550 if (printNum&&(ievt%printNum==0)) {
551 AliInfo(Form("[%s] EventMixing %d/%d",GetName(),ievt,nEvents));
552 timer.Stop(); timer.Print(); timer.Start(kFALSE); fflush(stdout);
555 fEvBuffer->GetEntry(ievt);
556 AliRsnMiniEvent evMain(*fMiniEvent);
557 list = smatched[ievt].Tokenize("|");
558 TObjArrayIter next(list);
559 while ( (os = (TObjString *)next()) ) {
560 imix = os->GetString().Atoi();
561 fEvBuffer->GetEntry(imix);
562 for (idef = 0; idef < nDefs; idef++) {
563 def = (AliRsnMiniOutput *)fHistograms[idef];
565 if (!def->IsTrackPairMix()) continue;
566 ifill += def->FillPair(&evMain, fMiniEvent, &fValues, kTRUE);
567 if (!def->IsSymmetric()) {
568 AliDebugClass(2, "Reflecting non symmetric pair");
569 ifill += def->FillPair(fMiniEvent, &evMain, &fValues, kFALSE);
578 AliInfo(Form("[%s] EventMixing %d/%d",GetName(),nEvents,nEvents));
579 timer.Stop(); timer.Print(); fflush(stdout);
584 for (iloop = 1; iloop < nEvents; iloop++) {
586 // restart from beginning if reached last event
587 if (imix >= nEvents) imix -= nEvents;
588 // avoid to mix an event with itself
589 if (imix == ievt) continue;
590 // skip all events already mixed enough times
591 if (fNMixed[ievt] >= fNMix) break;
592 if (fNMixed[imix] >= fNMix) continue;
593 fEvBuffer->GetEntry(imix);
594 // skip if events are not matched
595 if (TMath::Abs(evMain.Vz() - fMiniEvent->Vz() ) > fMaxDiffVz ) continue;
596 if (TMath::Abs(evMain.Mult() - fMiniEvent->Mult() ) > fMaxDiffMult ) continue;
597 if (TMath::Abs(evMain.Angle() - fMiniEvent->Angle()) > fMaxDiffAngle) continue;
598 // found a match: increment counter for both events
599 AliDebugClass(1, Form("Event %d, def '%s': event mixing (%d with %d)", ievt, def->GetName(), ievt, imix));
603 ifill += def->FillPair(&evMain, fMiniEvent, &fValues);
604 // stop if mixed enough times
605 if (fNMixed[ievt] >= fNMix) break;
608 // print number of mixings done with each event
609 for (ievt = 0; ievt < nEvents; ievt++) {
610 AliDebugClass(2, Form("Event %6d: mixed %2d times", ievt, fNMixed[ievt]));
614 // post computed data
615 PostData(1, fOutput);
618 //__________________________________________________________________________________________________
619 void AliRsnMiniAnalysisTask::Terminate(Option_t *)
622 // Draw result to screen, or perform fitting, normalizations
623 // Called once at the end of the query
626 fOutput = dynamic_cast<TList *>(GetOutputData(1));
628 AliError("Could not retrieve TList fOutput");
633 //__________________________________________________________________________________________________
634 Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent()
637 // This method checks if current event is OK for analysis.
638 // In case it is, the pointers of the local AliRsnEvent data member
639 // will point to it, in order to allow cut checking, otherwise the
640 // function exits with a failure message.
642 // ESD events must pass the physics selection, AOD are supposed to do.
644 // While checking the event, a histogram is filled to count the number
645 // of CINT1B, V0AND and CANDLE events, which are needed for normalization
647 // Return values can be:
648 // -- 'E' if the event is accepted and is ESD
649 // -- 'A' if the event is accepted and is AOD
650 // -- 0 if the event is not accepted
653 // string to sum messages
657 // exit points are provided in all cases an event is bad
658 // if this block is passed, an event can be rejected only
659 // if it does not pass the set of event cuts defined in the task
662 if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
665 // ESD specific check: Physics Selection
666 // --> if this is failed, the event is rejected
667 isSelected = (((AliInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
670 AliDebugClass(2, "Event does not pass physics selections");
671 fRsnEvent.SetRef(0x0);
672 fRsnEvent.SetRefMC(0x0);
675 // set reference to input
676 fRsnEvent.SetRef(fInputEvent);
677 // add MC if requested and available
680 fRsnEvent.SetRefMC(fMCEvent);
682 AliWarning("MC event requested but not available");
683 fRsnEvent.SetRefMC(0x0);
686 } else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
689 // set reference to input
690 fRsnEvent.SetRef(fInputEvent);
691 // add MC if requested and available (it is in the same object)
693 fRsnEvent.SetRefMC(fInputEvent);
694 if (!fRsnEvent.GetAODList()) {
695 AliWarning("MC event requested but not available");
696 fRsnEvent.SetRefMC(0x0);
700 AliError(Form("Bad input event class: %s", fInputEvent->ClassName()));
701 // reset pointers in local AliRsnEvent object
702 fRsnEvent.SetRef(0x0);
703 fRsnEvent.SetRefMC(0x0);
707 // fill counter of accepted events
708 fHEventStat->Fill(0.1);
710 // check if it is V0AND
711 // --> uses a cast to AliESDEvent even if the input is an AliAODEvent
712 Bool_t v0A = fTriggerAna->IsOfflineTriggerFired((AliESDEvent *)fInputEvent, AliTriggerAnalysis::kV0A);
713 Bool_t v0C = fTriggerAna->IsOfflineTriggerFired((AliESDEvent *)fInputEvent, AliTriggerAnalysis::kV0C);
715 msg += " -- VOAND = YES";
716 fHEventStat->Fill(1.1);
718 msg += " -- VOAND = NO ";
722 // --> requires at least one good quality track with Pt > 0.5 and |eta| <= 0.8
723 Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
724 Bool_t candle = kFALSE;
725 for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
726 AliVTrack *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
727 AliESDtrack *esdt = dynamic_cast<AliESDtrack *>(track);
728 AliAODTrack *aodt = dynamic_cast<AliAODTrack *>(track);
729 if (track->Pt() < 0.5) continue;
730 if(TMath::Abs(track->Eta()) > 0.8) continue;
731 if (esdt) if (!fESDtrackCuts->AcceptTrack(esdt)) continue;
732 if (aodt) if (!aodt->TestFilterBit(5)) continue;
737 msg += " -- CANDLE = YES";
738 fHEventStat->Fill(2.1);
740 msg += " -- CANDLE = NO ";
743 // if event cuts are defined, they are checked here
744 // final decision on the event depends on this
747 if (!fEventCuts->IsSelected(&fRsnEvent)) {
748 msg += " -- Local cuts = REJECTED";
751 msg += " -- Local cuts = ACCEPTED";
755 msg += " -- Local cuts = NONE";
759 // if the above exit point is not taken, the event is accepted
760 AliDebugClass(2, Form("Stats: %s", msg.Data()));
762 fHEventStat->Fill(3.1);
763 Double_t multi = ComputeCentrality((output == 'E'));
764 Double_t tracklets = ComputeTracklets();
765 fHAEventsVsMulti->Fill(multi);
766 fHAEventsVsTracklets->Fill(tracklets);
767 if(fHAEventVz) fHAEventVz->Fill(multi,fInputEvent->GetPrimaryVertex()->GetZ());
768 if(fHAEventMultiCent) fHAEventMultiCent->Fill(multi,ComputeMultiplicity(output == 'E',fHAEventMultiCent->GetYaxis()->GetTitle()));
769 if(fHAEventPlane) fHAEventPlane->Fill(multi,ComputeAngle());
772 fHEventStat->Fill(4.1);
773 const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
774 if(!vertex) fHEventStat->Fill(5.1);
776 TString title=vertex->GetTitle();
777 if( (title.Contains("Z")) || (title.Contains("3D")) ) fHEventStat->Fill(5.1);
778 if(vertex->GetNContributors()<1.) fHEventStat->Fill(6.1);
779 if(TMath::Abs(vertex->GetZ())>10.) fHEventStat->Fill(7.1);
785 //__________________________________________________________________________________________________
786 void AliRsnMiniAnalysisTask::FillMiniEvent(Char_t evType)
789 // Refresh cursor mini-event data member to fill with current event.
790 // Returns the total number of tracks selected.
793 // assign event-related values
794 if (fMiniEvent) delete fMiniEvent;
795 fMiniEvent = new AliRsnMiniEvent;
796 fMiniEvent->Vz() = fInputEvent->GetPrimaryVertex()->GetZ();
797 fMiniEvent->Angle() = ComputeAngle();
798 fMiniEvent->Mult() = ComputeCentrality((evType == 'E'));
799 fMiniEvent->Tracklets() = ComputeTracklets();
800 AliDebugClass(2, Form("Event %d: type = %c -- vz = %f -- mult = %f -- angle = %f", fEvNum, evType, fMiniEvent->Vz(), fMiniEvent->Mult(), fMiniEvent->Angle()));
802 // loop on daughters and assign track-related values
803 Int_t ic, ncuts = fTrackCuts.GetEntries();
804 Int_t ip, npart = fRsnEvent.GetAbsoluteSum();
805 Int_t npos = 0, nneg = 0, nneu = 0;
806 AliRsnDaughter cursor;
807 AliRsnMiniParticle miniParticle;
808 for (ip = 0; ip < npart; ip++) {
809 // point cursor to next particle
810 fRsnEvent.SetDaughter(cursor, ip);
811 // copy momentum and MC info if present
812 miniParticle.CopyDaughter(&cursor);
813 miniParticle.Index() = ip;
814 // switch on the bits corresponding to passed cuts
815 for (ic = 0; ic < ncuts; ic++) {
816 AliRsnCutSet *cuts = (AliRsnCutSet *)fTrackCuts[ic];
817 if (cuts->IsSelected(&cursor)) miniParticle.SetCutBit(ic);
819 // if a track passes at least one track cut, it is added to the pool
820 if (miniParticle.CutBits()) {
821 fMiniEvent->AddParticle(miniParticle);
822 if (miniParticle.Charge() == '+') npos++;
823 else if (miniParticle.Charge() == '-') nneg++;
828 // get number of accepted tracks
829 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));
832 //__________________________________________________________________________________________________
833 Double_t AliRsnMiniAnalysisTask::ComputeAngle()
836 // Get the plane angle
839 AliEventplane *plane = 0x0;
841 if (fInputEvent->InheritsFrom(AliESDEvent::Class()))
842 plane = fInputEvent->GetEventplane();
843 else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
844 AliAODEvent *aodEvent = (AliAODEvent *)fInputEvent;
845 plane = ((AliVAODHeader*)aodEvent->GetHeader())->GetEventplaneP();
849 return plane->GetEventplane("Q");
851 AliWarning("No event plane defined");
856 //__________________________________________________________________________________________________
857 Double_t AliRsnMiniAnalysisTask::ComputeCentrality(Bool_t isESD)
860 // Computes event centrality/multiplicity according to the criterion defined
861 // by two elements: (1) choice between multiplicity and centrality and
862 // (2) the string defining what criterion must be used for specific computation.
865 if (fUseCentrality) {
866 if ((!fUseMC) && (fUseCentralityPatchPbPb2011)) {
867 return ApplyCentralityPatchPbPb2011();//
869 if ((!fUseMC) && (!isESD) && (fUseAOD049CentralityPatch)) {
870 return ApplyCentralityPatchAOD049();
872 AliCentrality *centrality = fInputEvent->GetCentrality();
874 AliError("Cannot compute centrality!");
877 return centrality->GetCentralityPercentile(fCentralityType.Data());
880 if (!fCentralityType.CompareTo("TRACKS"))
881 return fInputEvent->GetNumberOfTracks();
882 else if (!fCentralityType.CompareTo("QUALITY"))
884 return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent *)fInputEvent, kTRUE);
887 Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
888 for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
889 AliVTrack *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
890 AliAODTrack *aodt = dynamic_cast<AliAODTrack *>(track);
892 if (!aodt->TestFilterBit(5)) continue;
897 else if (!fCentralityType.CompareTo("TRACKLETS")) {
899 const AliMultiplicity *mult = ((AliESDEvent *)fInputEvent)->GetMultiplicity();
900 Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
901 for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
902 return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ());
904 AliWarning("Cannot compute multiplicity with SPD tracklets from AOD");
908 AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", fCentralityType.Data()));
914 //__________________________________________________________________________________________________
915 Double_t AliRsnMiniAnalysisTask::ComputeMultiplicity(Bool_t isESD,TString type)
918 // Computes event multiplicity according to the string defining
919 // what criterion must be used for specific computation.
924 if (!type.CompareTo("TRACKS"))
925 return fInputEvent->GetNumberOfTracks();
926 else if (!type.CompareTo("QUALITY"))
928 return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent *)fInputEvent, kTRUE);
931 Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
932 for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
933 AliVTrack *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
934 AliAODTrack *aodt = dynamic_cast<AliAODTrack *>(track);
936 if (!aodt->TestFilterBit(5)) continue;
941 else if (!type.CompareTo("TRACKLETS")) {
943 const AliMultiplicity *mult = ((AliESDEvent *)fInputEvent)->GetMultiplicity();
944 Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
945 for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
946 return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ());
948 AliWarning("Cannot compute multiplicity with SPD tracklets from AOD");
952 AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", type.Data()));
957 //__________________________________________________________________________________________________
958 Double_t AliRsnMiniAnalysisTask::ComputeTracklets()
961 // Get number of tracklets
964 Double_t count = 100;
966 if (fInputEvent->InheritsFrom(AliESDEvent::Class())){
967 AliESDEvent *esdEvent = (AliESDEvent *)fInputEvent;
968 const AliMultiplicity *spdmult = esdEvent->GetMultiplicity();
969 count = 1.0*spdmult->GetNumberOfTracklets();
971 else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
972 AliAODEvent *aodEvent = (AliAODEvent *)fInputEvent;
973 AliAODTracklets *spdmult = aodEvent->GetTracklets();
974 count = 1.0*spdmult->GetNumberOfTracklets();
980 //__________________________________________________________________________________________________
981 void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
984 // Fills the histograms with true mother (ESD version)
988 Int_t id, ndef = fHistograms.GetEntries();
989 Int_t ip, label1, label2, npart = fMCEvent->GetNumberOfTracks();
990 static AliRsnMiniPair miniPair;
991 AliMCParticle *daughter1, *daughter2;
992 TLorentzVector p1, p2;
993 AliRsnMiniOutput *def = 0x0;
995 for (id = 0; id < ndef; id++) {
996 def = (AliRsnMiniOutput *)fHistograms[id];
998 if (!def->IsMother() && !def->IsMotherInAcc()) continue;
999 for (ip = 0; ip < npart; ip++) {
1000 AliMCParticle *part = (AliMCParticle *)fMCEvent->GetTrack(ip);
1001 //get mother pdg code
1002 if (part->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
1003 // check that daughters match expected species
1004 if (part->Particle()->GetNDaughters() < 2) continue;
1005 if (fMaxNDaughters > 0 && part->Particle()->GetNDaughters() > fMaxNDaughters) continue;
1006 label1 = part->Particle()->GetDaughter(0);
1007 label2 = part->Particle()->GetDaughter(1);
1008 daughter1 = (AliMCParticle *)fMCEvent->GetTrack(label1);
1009 daughter2 = (AliMCParticle *)fMCEvent->GetTrack(label2);
1011 if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(1)) {
1013 p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
1014 p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
1015 } else if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(0)) {
1017 p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
1018 p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
1020 if (!okMatch) continue;
1021 if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
1022 (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
1023 (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
1025 Int_t pdgGranma = 0;
1027 mother = part->GetMother();
1029 Int_t abspdgGranma =0;
1030 Bool_t isFromB=kFALSE;
1031 Bool_t isQuarkFound=kFALSE;
1032 while (mother >=0 ){
1034 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1035 AliMCParticle* mcGranma = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(mother));
1037 pdgGranma = mcGranma->PdgCode();
1038 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1039 abspdgGranma = TMath::Abs(pdgGranma);
1040 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1043 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1044 mother = mcGranma->GetMother();
1046 AliError("Failed casting the mother particle!");
1050 if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
1052 if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
1055 if (fKeepDfromBOnly) pdgGranma = -999;
1057 if (pdgGranma == -99999){
1058 AliDebug(2,"This particle does not have a quark in his genealogy\n");
1061 if (pdgGranma == -9999){
1062 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");
1066 if (pdgGranma == -999){
1067 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");
1073 // assign momenta to computation object
1074 miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
1075 miniPair.FillRef(def->GetMotherMass());
1077 def->FillMother(&miniPair, miniEvent, &fValues);
1078 if(fKeepMotherInAcceptance){
1079 if(daughter1->Pt()<fMotherAcceptanceCutMinPt || daughter2->Pt()<fMotherAcceptanceCutMinPt || TMath::Abs(daughter1->Eta())>fMotherAcceptanceCutMaxEta || TMath::Abs(daughter2->Eta())>fMotherAcceptanceCutMaxEta) continue;
1080 def->FillMotherInAcceptance(&miniPair, miniEvent, &fValues);
1088 //__________________________________________________________________________________________________
1089 void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
1092 // Fills the histograms with true mother (AOD version)
1096 TClonesArray *list = fRsnEvent.GetAODList();
1097 Int_t id, ndef = fHistograms.GetEntries();
1098 Int_t ip, label1, label2, npart = list->GetEntries();
1099 static AliRsnMiniPair miniPair;
1100 AliAODMCParticle *daughter1, *daughter2;
1101 TLorentzVector p1, p2;
1102 AliRsnMiniOutput *def = 0x0;
1104 for (id = 0; id < ndef; id++) {
1105 def = (AliRsnMiniOutput *)fHistograms[id];
1107 if (!def->IsMother() && !def->IsMotherInAcc()) continue;
1108 for (ip = 0; ip < npart; ip++) {
1109 AliAODMCParticle *part = (AliAODMCParticle *)list->At(ip);
1110 if (part->GetPdgCode() != def->GetMotherPDG()) continue;
1111 // check that daughters match expected species
1112 if (part->GetNDaughters() < 2) continue;
1113 if (fMaxNDaughters > 0 && part->GetNDaughters() > fMaxNDaughters) continue;
1114 label1 = part->GetDaughter(0);
1115 label2 = part->GetDaughter(1);
1116 daughter1 = (AliAODMCParticle *)list->At(label1);
1117 daughter2 = (AliAODMCParticle *)list->At(label2);
1119 if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(1)) {
1121 p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
1122 p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
1123 } else if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(0)) {
1125 p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
1126 p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
1128 if (!okMatch) continue;
1129 if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
1130 (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
1131 (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
1133 Int_t pdgGranma = 0;
1135 mother = part->GetMother();
1137 Int_t abspdgGranma =0;
1138 Bool_t isFromB=kFALSE;
1139 Bool_t isQuarkFound=kFALSE;
1140 while (mother >=0 ){
1142 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1143 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(list->At(mother));
1145 pdgGranma = mcGranma->GetPdgCode();
1146 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1147 abspdgGranma = TMath::Abs(pdgGranma);
1148 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1151 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1152 mother = mcGranma->GetMother();
1154 AliError("Failed casting the mother particle!");
1158 if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
1160 if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
1163 if (fKeepDfromBOnly) pdgGranma = -999;
1166 if (pdgGranma == -99999){
1167 AliDebug(2,"This particle does not have a quark in his genealogy\n");
1170 if (pdgGranma == -9999){
1171 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");
1175 if (pdgGranma == -999){
1176 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");
1180 // assign momenta to computation object
1181 miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
1182 miniPair.FillRef(def->GetMotherMass());
1184 def->FillMother(&miniPair, miniEvent, &fValues);
1185 if(fKeepMotherInAcceptance){
1186 if(daughter1->Pt()<fMotherAcceptanceCutMinPt || daughter2->Pt()<fMotherAcceptanceCutMinPt || TMath::Abs(daughter1->Eta())>fMotherAcceptanceCutMaxEta || TMath::Abs(daughter2->Eta())>fMotherAcceptanceCutMaxEta) continue;
1187 def->FillMotherInAcceptance(&miniPair, miniEvent, &fValues);
1192 //___________________________________________________________
1193 void AliRsnMiniAnalysisTask::SetDselection(UShort_t originDselection)
1195 // setting the way the D0 will be selected
1196 // 0 --> only from c quarks
1197 // 1 --> only from b quarks
1198 // 2 --> from both c quarks and b quarks
1200 fOriginDselection = originDselection;
1202 if (fOriginDselection == 0) {
1203 fKeepDfromB = kFALSE;
1204 fKeepDfromBOnly = kFALSE;
1207 if (fOriginDselection == 1) {
1208 fKeepDfromB = kTRUE;
1209 fKeepDfromBOnly = kTRUE;
1212 if (fOriginDselection == 2) {
1213 fKeepDfromB = kTRUE;
1214 fKeepDfromBOnly = kFALSE;
1219 //__________________________________________________________________________________________________
1220 Bool_t AliRsnMiniAnalysisTask::EventsMatch(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2)
1223 // Check if two events are compatible.
1224 // If the mixing is continuous, this is true if differences in vz, mult and angle are smaller than
1225 // the specified values.
1226 // If the mixing is binned, this is true if the events are in the same bin.
1229 if (!event1 || !event2) return kFALSE;
1230 Int_t ivz1, ivz2, imult1, imult2, iangle1, iangle2;
1231 Double_t dv, dm, da;
1233 if (fContinuousMix) {
1234 dv = TMath::Abs(event1->Vz() - event2->Vz() );
1235 dm = TMath::Abs(event1->Mult() - event2->Mult() );
1236 da = TMath::Abs(event1->Angle() - event2->Angle());
1237 if (dv > fMaxDiffVz) {
1238 //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Vz = %f", event1->ID(), event2->ID(), dv));
1241 if (dm > fMaxDiffMult ) {
1242 //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Mult = %f", event1->ID(), event2->ID(), dm));
1245 if (da > fMaxDiffAngle) {
1246 //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Angle = %f", event1->ID(), event2->ID(), da));
1251 ivz1 = (Int_t)(event1->Vz() / fMaxDiffVz);
1252 ivz2 = (Int_t)(event2->Vz() / fMaxDiffVz);
1253 imult1 = (Int_t)(event1->Mult() / fMaxDiffMult);
1254 imult2 = (Int_t)(event2->Mult() / fMaxDiffMult);
1255 iangle1 = (Int_t)(event1->Angle() / fMaxDiffAngle);
1256 iangle2 = (Int_t)(event2->Angle() / fMaxDiffAngle);
1257 if (ivz1 != ivz2) return kFALSE;
1258 if (imult1 != imult2) return kFALSE;
1259 if (iangle1 != iangle2) return kFALSE;
1264 //---------------------------------------------------------------------
1265 Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchPbPb2011(){
1266 //This part rejects randomly events such that the centrality gets flat for LHC11h Pb-Pb data
1267 //for 0-5% and 10-20% centrality bin
1269 if (fCentralityType!="V0M") {
1270 AliWarning("Wrong value (not centrality from V0).");
1274 AliCentrality *centrality = fInputEvent->GetCentrality();
1276 AliWarning("Cannot get centrality from AOD event.");
1280 Double_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
1281 Double_t rnd_hc = -1., testf = 0.0, ff = 0, N1 = -1., N2 = -1.;
1283 if(fUseCentralityPatchPbPb2011==510){
1285 N2 = 1.56435e+06; //N2 is the reference
1286 ff = 5.04167e+06 - 1.49885e+06*cent + 2.35998e+05*cent*cent -1.22873e+04*cent*cent*cent;
1288 if(fUseCentralityPatchPbPb2011==1020){
1289 N2 = 2.0e+05; //N2 is the reference
1291 ff = -1.73979e+06 - 3.05316e+06*cent + 1.05517e+06*cent*cent - 133205*cent*cent*cent + 8187.45*cent*cent*cent*cent - 247.875*cent*cent*cent*cent*cent + 2.9676*cent*cent*cent*cent*cent*cent;
1293 AliError(Form("Patch for the requested centrality (%i) is not available", fUseCentralityPatchPbPb2011));
1297 testf = ( N2 + (N1-ff) ) / N1;
1298 rnd_hc = gRandom->Rndm();
1300 //AliDebugClass(1, Form("Flat Centrality %d", fUseCentralityPatchPbPb2011));
1302 if (rnd_hc < 0 || rnd_hc > 1 )
1304 AliWarning("Wrong Random number generated");
1313 //---------------------------------------------------------------------
1314 Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchAOD049()
1317 //Apply centrality patch for AOD049 outliers
1319 if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
1320 AliWarning("Requested patch for AOD049 for ESD. ");
1324 if (fCentralityType!="V0M") {
1325 AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
1329 AliCentrality *centrality = fInputEvent->GetCentrality();
1331 AliWarning("Cannot get centrality from AOD event.");
1335 Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
1337 Bool_t isSelRun = kFALSE;
1338 Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
1340 Int_t quality = centrality->GetQuality();
1342 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1344 Int_t runnum=aodEvent->GetRunNumber();
1345 for(Int_t ir=0;ir<5;ir++){
1346 if(runnum==selRun[ir]){
1351 if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1358 AliAODEvent *aodEvent = (AliAODEvent *)fInputEvent;
1359 AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
1360 v0+=aodV0->GetMTotV0A();
1361 v0+=aodV0->GetMTotV0C();
1362 if ( (cent==0) && (v0<19500) ) {
1363 AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
1366 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1367 Float_t val = 1.30552 + 0.147931 * v0;
1369 Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
1370 120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
1371 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
1372 68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
1373 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
1374 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
1375 26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
1376 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
1377 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
1378 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
1381 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] ) {
1382 AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
1386 //force it to be -999. whatever the negative value was
1392 //----------------------------------------------------------------------------------
1393 void AliRsnMiniAnalysisTask::SetEventQAHist(TString type,TH2F *histo)
1396 AliWarning(Form("event QA histogram pointer not defined for slot %s",type.Data()));
1402 if(!type.CompareTo("vz")) fHAEventVz = histo;
1403 else if(!type.CompareTo("multicent")) {
1404 TString mtype(histo->GetYaxis()->GetTitle());
1406 if(mtype.CompareTo("QUALITY") && mtype.CompareTo("TRACKS") && mtype.CompareTo("TRACKLETS")) {
1407 AliWarning(Form("multiplicity vs. centrality histogram y-axis %s unknown, setting to TRACKS",mtype.Data()));
1408 histo->GetYaxis()->SetTitle("TRACKS");
1410 fHAEventMultiCent = histo;
1412 else if(!type.CompareTo("eventplane")) fHAEventPlane = histo;
1413 else AliWarning(Form("event QA histogram slot %s undefined",type.Data()));
1418 //----------------------------------------------------------------------------------
1419 Int_t AliRsnMiniAnalysisTask::CreateValue(AliRsnMiniValue::EType type, Bool_t useMC)
1422 // Create a new value in the task,
1423 // and returns its ID, which is needed for setting up histograms.
1424 // If that value was already initialized, returns its ID and does not recreate it.
1427 Int_t valID = ValueID(type, useMC);
1428 if (valID >= 0 && valID < fValues.GetEntries()) {
1429 AliInfo(Form("Value '%s' is already created in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
1431 valID = fValues.GetEntries();
1432 AliInfo(Form("Creating value '%s' in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
1433 new (fValues[valID]) AliRsnMiniValue(type, useMC);
1439 //----------------------------------------------------------------------------------
1440 Int_t AliRsnMiniAnalysisTask::ValueID(AliRsnMiniValue::EType type, Bool_t useMC)
1443 // Searches if a value computation is initialized
1446 const char *name = AliRsnMiniValue::ValueName(type, useMC);
1447 TObject *obj = fValues.FindObject(name);
1449 return fValues.IndexOf(obj);
1454 //----------------------------------------------------------------------------------
1455 AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, AliRsnMiniOutput::EOutputType type, AliRsnMiniOutput::EComputation src)
1458 // Create a new histogram definition in the task,
1459 // which is then returned to the user for its configuration
1462 Int_t n = fHistograms.GetEntries();
1463 AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, type, src);
1468 //----------------------------------------------------------------------------------
1469 AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, const char *outType, const char *compType)
1472 // Create a new histogram definition in the task,
1473 // which is then returned to the user for its configuration
1476 Int_t n = fHistograms.GetEntries();
1477 AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, outType, compType);