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