PWGLF/RESONANCES/AliRsnMiniAnalysisTask.cxx

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