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