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