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          //get mother pdg code
 904          if (part->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
 905          // check that daughters match expected species
 906          if (part->Particle()->GetNDaughters() < 2) continue;
 907          label1 = part->Particle()->GetDaughter(0);
 908          label2 = part->Particle()->GetDaughter(1);
 909          daughter1 = (AliMCParticle *)fMCEvent->GetTrack(label1);
 910          daughter2 = (AliMCParticle *)fMCEvent->GetTrack(label2);
 911          okMatch = kFALSE;
 912          if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(1)) {
 913             okMatch = kTRUE;
 914             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
 915             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
 916          } else if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(0)) {
 917             okMatch = kTRUE;
 918             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
 919             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
 920          }
 921          if (!okMatch) continue;
 922          // assign momenta to computation object
 923          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
 924          miniPair.FillRef(def->GetMotherMass());
 925          // do computations
 926          def->FillMother(&miniPair, miniEvent, &fValues);
 927       }
 928    }
 929 }
 930
 931 //__________________________________________________________________________________________________
 932 void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
 933 {
 934 //
 935 // Fills the histograms with true mother (AOD version)
 936 //
 937
 938    Bool_t okMatch;
 939    TClonesArray *list = fRsnEvent.GetAODList();
 940    Int_t id, ndef = fHistograms.GetEntries();
 941    Int_t ip, label1, label2, npart = list->GetEntries();
 942    static AliRsnMiniPair miniPair;
 943    AliAODMCParticle *daughter1, *daughter2;
 944    TLorentzVector p1, p2;
 945    AliRsnMiniOutput *def = 0x0;
 946
 947    for (id = 0; id < ndef; id++) {
 948       def = (AliRsnMiniOutput *)fHistograms[id];
 949       if (!def) continue;
 950       if (!def->IsMother()) continue;
 951       for (ip = 0; ip < npart; ip++) {
 952          AliAODMCParticle *part = (AliAODMCParticle *)list->At(ip);
 953          if (part->GetPdgCode() != def->GetMotherPDG()) continue;
 954          // check that daughters match expected species
 955          if (part->GetNDaughters() < 2) continue;
 956          label1 = part->GetDaughter(0);
 957          label2 = part->GetDaughter(1);
 958          daughter1 = (AliAODMCParticle *)list->At(label1);
 959          daughter2 = (AliAODMCParticle *)list->At(label2);
 960          okMatch = kFALSE;
 961          if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(1)) {
 962             okMatch = kTRUE;
 963             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
 964             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
 965          } else if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(0)) {
 966             okMatch = kTRUE;
 967             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
 968             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
 969          }
 970          if (!okMatch) continue;
 971          // assign momenta to computation object
 972          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
 973          miniPair.FillRef(def->GetMotherMass());
 974          // do computations
 975          def->FillMother(&miniPair, miniEvent, &fValues);
 976       }
 977    }
 978 }
 979
 980 //__________________________________________________________________________________________________
 981 Bool_t AliRsnMiniAnalysisTask::EventsMatch(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2)
 982 {
 983 //
 984 // Check if two events are compatible.
 985 // If the mixing is continuous, this is true if differences in vz, mult and angle are smaller than
 986 // the specified values.
 987 // If the mixing is binned, this is true if the events are in the same bin.
 988 //
 989
 990    if (!event1 || !event2) return kFALSE;
 991    Int_t ivz1, ivz2, imult1, imult2, iangle1, iangle2;
 992    Double_t dv, dm, da;
 993
 994    if (fContinuousMix) {
 995       dv = TMath::Abs(event1->Vz()    - event2->Vz()   );
 996       dm = TMath::Abs(event1->Mult()  - event2->Mult() );
 997       da = TMath::Abs(event1->Angle() - event2->Angle());
 998       if (dv > fMaxDiffVz) {
 999          //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Vz = %f", event1->ID(), event2->ID(), dv));
1000          return kFALSE;
1001       }
1002       if (dm > fMaxDiffMult ) {
1003          //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Mult = %f", event1->ID(), event2->ID(), dm));
1004          return kFALSE;
1005       }
1006       if (da > fMaxDiffAngle) {
1007          //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Angle = %f", event1->ID(), event2->ID(), da));
1008          return kFALSE;
1009       }
1010       return kTRUE;
1011    } else {
1012       ivz1 = (Int_t)(event1->Vz() / fMaxDiffVz);
1013       ivz2 = (Int_t)(event2->Vz() / fMaxDiffVz);
1014       imult1 = (Int_t)(event1->Mult() / fMaxDiffMult);
1015       imult2 = (Int_t)(event2->Mult() / fMaxDiffMult);
1016       iangle1 = (Int_t)(event1->Angle() / fMaxDiffAngle);
1017       iangle2 = (Int_t)(event2->Angle() / fMaxDiffAngle);
1018       if (ivz1 != ivz2) return kFALSE;
1019       if (imult1 != imult2) return kFALSE;
1020       if (iangle1 != iangle2) return kFALSE;
1021       return kTRUE;
1022    }
1023 }
1024
1025
1026 //---------------------------------------------------------------------
1027 Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchAOD049()
1028 {
1029    //
1030    //Apply centrality patch for AOD049 outliers
1031    //
1032    if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
1033       AliWarning("Requested patch for AOD049 for ESD. ");
1034       return -999.0;
1035    }
1036
1037    if (fCentralityType!="V0M") {
1038       AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
1039       return -999.0;
1040    }
1041
1042    AliCentrality *centrality = fInputEvent->GetCentrality();
1043    if (!centrality) {
1044       AliWarning("Cannot get centrality from AOD event.");
1045       return -999.0;
1046    }
1047
1048    Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
1049    /*
1050    Bool_t isSelRun = kFALSE;
1051    Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
1052    if(cent<0){
1053      Int_t quality = centrality->GetQuality();
1054      if(quality<=1){
1055        cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1056      } else {
1057        Int_t runnum=aodEvent->GetRunNumber();
1058        for(Int_t ir=0;ir<5;ir++){
1059    if(runnum==selRun[ir]){
1060      isSelRun=kTRUE;
1061      break;
1062    }
1063        }
1064        if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1065      }
1066    }
1067    */
1068
1069    if(cent>=0.0) {
1070       Float_t v0 = 0.0;
1071       AliAODEvent *aodEvent = (AliAODEvent *)fInputEvent;
1072       AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
1073       v0+=aodV0->GetMTotV0A();
1074       v0+=aodV0->GetMTotV0C();
1075       if ( (cent==0) && (v0<19500) ) {
1076          AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
1077          return -999.0;
1078       }
1079       Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1080       Float_t val = 1.30552 +  0.147931 * v0;
1081
1082       Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
1083                                120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
1084                                92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
1085                                68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
1086                                51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
1087                                37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
1088                                26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
1089                                19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
1090                                12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
1091                                13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
1092                               };
1093
1094       if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )  {
1095          AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
1096          return -999.0;
1097       }
1098    } else {
1099       //force it to be -999. whatever the negative value was
1100       cent = -999.;
1101    }
1102    return cent;
1103 }
1104
1105 //----------------------------------------------------------------------------------
1106 void AliRsnMiniAnalysisTask::SetEventQAHist(TString type,TH2F *histo)
1107 {
1108    if(!histo) {
1109       AliWarning(Form("event QA histogram pointer not defined for slot %s",type.Data()));
1110       return;
1111    }
1112
1113    type.ToLower();
1114
1115    if(!type.CompareTo("vz")) fHAEventVz = histo;
1116    else if(!type.CompareTo("multicent")) {
1117       TString mtype(histo->GetYaxis()->GetTitle());
1118       mtype.ToUpper();
1119       if(mtype.CompareTo("QUALITY") && mtype.CompareTo("TRACKS") && mtype.CompareTo("TRACKLETS")) {
1120          AliWarning(Form("multiplicity vs. centrality histogram y-axis %s unknown, setting to TRACKS",mtype.Data()));
1121          histo->GetYaxis()->SetTitle("TRACKS");
1122       }
1123       fHAEventMultiCent = histo;
1124    }
1125    else if(!type.CompareTo("eventplane")) fHAEventPlane = histo;
1126    else AliWarning(Form("event QA histogram slot %s undefined",type.Data()));
1127
1128    return;
1129 }