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