PWGLF/RESONANCES/AliRsnMiniAnalysisTask.cxx

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