PWGLF/RESONANCES/AliRsnMiniAnalysisTask.cxx

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