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