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 #include "TRandom.h"
  20
  21 #include "AliLog.h"
  22 #include "AliEventplane.h"
  23 #include "AliMultiplicity.h"
  24 #include "AliTriggerAnalysis.h"
  25 #include "AliAnalysisManager.h"
  26 #include "AliInputEventHandler.h"
  27
  28 #include "AliESDtrackCuts.h"
  29 #include "AliESDUtils.h"
  30
  31 #include "AliAODEvent.h"
  32 #include "AliAODMCParticle.h"
  33
  34 #include "AliRsnCutSet.h"
  35 #include "AliRsnMiniPair.h"
  36 #include "AliRsnMiniEvent.h"
  37 #include "AliRsnMiniParticle.h"
  38
  39 #include "AliRsnMiniAnalysisTask.h"
  40
  41
  42 ClassImp(AliRsnMiniAnalysisTask)
  43
  44 //__________________________________________________________________________________________________
  45 AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() :
  46    AliAnalysisTaskSE(),
  47    fUseMC(kFALSE),
  48    fEvNum(0),
  49    fTriggerMask(0),
  50    fUseCentrality(kFALSE),
  51    fCentralityType("QUALITY"),
  52    fUseAOD049CentralityPatch(kFALSE),
  53    fUseCentralityPatchPbPb2011(0),
  54    fContinuousMix(kTRUE),
  55    fNMix(0),
  56    fMaxDiffMult(10),
  57    fMaxDiffVz(1.0),
  58    fMaxDiffAngle(1E20),
  59    fOutput(0x0),
  60    fHistograms("AliRsnMiniOutput", 0),
  61    fValues("AliRsnMiniValue", 0),
  62    fHEventStat(0x0),
  63    fHAEventsVsMulti(0x0),
  64    fHAEventVz(0x0),
  65    fHAEventMultiCent(0x0),
  66    fHAEventPlane(0x0),
  67    fEventCuts(0x0),
  68    fTrackCuts(0),
  69    fRsnEvent(),
  70    fEvBuffer(0x0),
  71    fTriggerAna(0x0),
  72    fESDtrackCuts(0x0),
  73    fMiniEvent(0x0),
  74    fBigOutput(kFALSE),
  75    fMixPrintRefresh(-1),
  76    fMaxNDaughters(-1),
  77    fCheckP(kFALSE),
  78    fCheckFeedDown(kFALSE),
  79    fOriginDselection(kFALSE),
  80    fKeepDfromB(kFALSE),
  81    fKeepDfromBOnly(kFALSE),
  82    fRejectIfNoQuark(kFALSE)
  83 {
  84 //
  85 // Dummy constructor ALWAYS needed for I/O.
  86 //
  87 }
  88
  89 //__________________________________________________________________________________________________
  90 AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name, Bool_t useMC) :
  91    AliAnalysisTaskSE(name),
  92    fUseMC(useMC),
  93    fEvNum(0),
  94    fTriggerMask(AliVEvent::kMB),
  95    fUseCentrality(kFALSE),
  96    fCentralityType("QUALITY"),
  97    fUseAOD049CentralityPatch(kFALSE),
  98   fUseCentralityPatchPbPb2011(0),
  99    fContinuousMix(kTRUE),
 100    fNMix(0),
 101    fMaxDiffMult(10),
 102    fMaxDiffVz(1.0),
 103    fMaxDiffAngle(1E20),
 104    fOutput(0x0),
 105    fHistograms("AliRsnMiniOutput", 0),
 106    fValues("AliRsnMiniValue", 0),
 107    fHEventStat(0x0),
 108    fHAEventsVsMulti(0x0),
 109    fHAEventVz(0x0),
 110    fHAEventMultiCent(0x0),
 111    fHAEventPlane(0x0),
 112    fEventCuts(0x0),
 113    fTrackCuts(0),
 114    fRsnEvent(),
 115    fEvBuffer(0x0),
 116    fTriggerAna(0x0),
 117    fESDtrackCuts(0x0),
 118    fMiniEvent(0x0),
 119    fBigOutput(kFALSE),
 120    fMixPrintRefresh(-1),
 121    fMaxNDaughters(-1),
 122    fCheckP(kFALSE),
 123    fCheckFeedDown(kFALSE),
 124    fOriginDselection(kFALSE),
 125    fKeepDfromB(kFALSE),
 126    fKeepDfromBOnly(kFALSE),
 127    fRejectIfNoQuark(kFALSE)
 128 {
 129 //
 130 // Default constructor.
 131 // Define input and output slots here (never in the dummy constructor)
 132 // Input slot #0 works with a TChain - it is connected to the default input container
 133 // Output slot #1 writes into a TH1 container
 134 //
 135
 136    DefineOutput(1, TList::Class());
 137 }
 138
 139 //__________________________________________________________________________________________________
 140 AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask &copy) :
 141    AliAnalysisTaskSE(copy),
 142    fUseMC(copy.fUseMC),
 143    fEvNum(0),
 144    fTriggerMask(copy.fTriggerMask),
 145    fUseCentrality(copy.fUseCentrality),
 146    fCentralityType(copy.fCentralityType),
 147    fUseAOD049CentralityPatch(copy.fUseAOD049CentralityPatch),
 148    fUseCentralityPatchPbPb2011(copy.fUseCentralityPatchPbPb2011),
 149    fContinuousMix(copy.fContinuousMix),
 150    fNMix(copy.fNMix),
 151    fMaxDiffMult(copy.fMaxDiffMult),
 152    fMaxDiffVz(copy.fMaxDiffVz),
 153    fMaxDiffAngle(copy.fMaxDiffAngle),
 154    fOutput(0x0),
 155    fHistograms(copy.fHistograms),
 156    fValues(copy.fValues),
 157    fHEventStat(0x0),
 158    fHAEventsVsMulti(0x0),
 159    fHAEventVz(0x0),
 160    fHAEventMultiCent(0x0),
 161    fHAEventPlane(0x0),
 162    fEventCuts(copy.fEventCuts),
 163    fTrackCuts(copy.fTrackCuts),
 164    fRsnEvent(),
 165    fEvBuffer(0x0),
 166    fTriggerAna(copy.fTriggerAna),
 167    fESDtrackCuts(copy.fESDtrackCuts),
 168    fMiniEvent(0x0),
 169    fBigOutput(copy.fBigOutput),
 170    fMixPrintRefresh(copy.fMixPrintRefresh),
 171    fMaxNDaughters(copy.fMaxNDaughters),
 172    fCheckP(copy.fCheckP),
 173    fCheckFeedDown(copy.fCheckFeedDown),
 174    fOriginDselection(copy.fOriginDselection),
 175    fKeepDfromB(copy.fOriginDselection),
 176    fKeepDfromBOnly(copy.fKeepDfromBOnly),
 177    fRejectIfNoQuark(copy.fRejectIfNoQuark)
 178 {
 179 //
 180 // Copy constructor.
 181 // Implemented as requested by C++ standards.
 182 // Can be used in PROOF and by plugins.
 183 //
 184 }
 185
 186 //__________________________________________________________________________________________________
 187 AliRsnMiniAnalysisTask &AliRsnMiniAnalysisTask::operator=(const AliRsnMiniAnalysisTask &copy)
 188 {
 189 //
 190 // Assignment operator.
 191 // Implemented as requested by C++ standards.
 192 // Can be used in PROOF and by plugins.
 193 //
 194
 195    AliAnalysisTaskSE::operator=(copy);
 196    if (this == &copy)
 197       return *this;
 198    fUseMC = copy.fUseMC;
 199    fEvNum = copy.fEvNum;
 200    fTriggerMask = copy.fTriggerMask;
 201    fUseCentrality = copy.fUseCentrality;
 202    fCentralityType = copy.fCentralityType;
 203    fUseAOD049CentralityPatch = copy.fUseAOD049CentralityPatch;
 204    fUseCentralityPatchPbPb2011 = copy.fUseCentralityPatchPbPb2011;
 205    fContinuousMix = copy.fContinuousMix;
 206    fNMix = copy.fNMix;
 207    fMaxDiffMult = copy.fMaxDiffMult;
 208    fMaxDiffVz = copy.fMaxDiffVz;
 209    fMaxDiffAngle = copy.fMaxDiffAngle;
 210    fHistograms = copy.fHistograms;
 211    fValues = copy.fValues;
 212    fHEventStat = copy.fHEventStat;
 213    fHAEventsVsMulti = copy.fHAEventsVsMulti;
 214    fHAEventVz = copy.fHAEventVz;
 215    fHAEventMultiCent = copy.fHAEventMultiCent;
 216    fHAEventPlane = copy.fHAEventPlane;
 217    fEventCuts = copy.fEventCuts;
 218    fTrackCuts = copy.fTrackCuts;
 219    fTriggerAna = copy.fTriggerAna;
 220    fESDtrackCuts = copy.fESDtrackCuts;
 221    fBigOutput = copy.fBigOutput;
 222    fMixPrintRefresh = copy.fMixPrintRefresh;
 223    fMaxNDaughters = copy.fMaxNDaughters;
 224    fCheckP = copy.fCheckP;
 225    fCheckFeedDown = copy.fCheckFeedDown;
 226    fOriginDselection = copy.fOriginDselection;
 227    fKeepDfromB = copy.fOriginDselection;
 228    fKeepDfromBOnly = copy.fKeepDfromBOnly;
 229    fRejectIfNoQuark = copy.fRejectIfNoQuark;
 230    return (*this);
 231 }
 232
 233 //__________________________________________________________________________________________________
 234 AliRsnMiniAnalysisTask::~AliRsnMiniAnalysisTask()
 235 {
 236 //
 237 // Destructor.
 238 // Clean-up the output list, but not the histograms that are put inside
 239 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 240 //
 241
 242
 243    if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 244       delete fOutput;
 245       delete fEvBuffer;
 246    }
 247 }
 248
 249 //__________________________________________________________________________________________________
 250 Int_t AliRsnMiniAnalysisTask::AddTrackCuts(AliRsnCutSet *cuts)
 251 {
 252 //
 253 // Add a new cut set for a new criterion for track selection.
 254 // A user can add as many as he wants, and each one corresponds
 255 // to one of the available bits in the AliRsnMiniParticle mask.
 256 // The only check is the following: if a cut set with the same name
 257 // as the argument is there, this is not added.
 258 // Return value is the array position of this set.
 259 //
 260
 261    TObject *obj = fTrackCuts.FindObject(cuts->GetName());
 262
 263    if (obj) {
 264       AliInfo(Form("A cut set named '%s' already exists", cuts->GetName()));
 265       return fTrackCuts.IndexOf(obj);
 266    } else {
 267       fTrackCuts.AddLast(cuts);
 268       return fTrackCuts.IndexOf(cuts);
 269    }
 270 }
 271
 272 //__________________________________________________________________________________________________
 273 void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
 274 {
 275 //
 276 // Initialization of outputs.
 277 // This is called once per worker node.
 278 //
 279
 280    // reset counter
 281    fEvNum = -1;
 282
 283    // message
 284    AliInfo(Form("Selected event characterization: %s (%s)", (fUseCentrality ? "centrality" : "multiplicity"), fCentralityType.Data()));
 285
 286    // initialize trigger analysis
 287    if (fTriggerAna) delete fTriggerAna;
 288    fTriggerAna = new AliTriggerAnalysis;
 289
 290    // initialize ESD quality cuts
 291    if (fESDtrackCuts) delete fESDtrackCuts;
 292    fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
 293
 294    // create list and set it as owner of its content (MANDATORY)
 295    if (fBigOutput) OpenFile(1);
 296    fOutput = new TList();
 297    fOutput->SetOwner();
 298
 299    // initialize event statistics counter
 300    fHEventStat = new TH1F("hEventStat", "Event statistics", 4, 0.0, 4.0);
 301    fHEventStat->GetXaxis()->SetBinLabel(1, "CINT1B");
 302    fHEventStat->GetXaxis()->SetBinLabel(2, "V0AND");
 303    fHEventStat->GetXaxis()->SetBinLabel(3, "Candle");
 304    fHEventStat->GetXaxis()->SetBinLabel(4, "Accepted");
 305    fOutput->Add(fHEventStat);
 306
 307    if (fUseCentrality)
 308       fHAEventsVsMulti = new TH1F("hAEventsVsMulti", "Accepted events vs Centrality", 100, 0, 100.0);
 309    else
 310       fHAEventsVsMulti = new TH1F("hAEventsVsMulti", "Accepted events vs Multiplicity",1000, 0, 1000.0);
 311    fOutput->Add(fHAEventsVsMulti);
 312
 313    if(fHAEventVz) fOutput->Add(fHAEventVz);
 314    if(fHAEventMultiCent) fOutput->Add(fHAEventMultiCent);
 315    if(fHAEventPlane) fOutput->Add(fHAEventPlane);
 316
 317    TIter next(&fTrackCuts);
 318    AliRsnCutSet *cs;
 319    while ((cs = (AliRsnCutSet *) next())) {
 320       cs->Init(fOutput);
 321    }
 322
 323    // create temporary tree for filtered events
 324    if (fMiniEvent) delete fMiniEvent;
 325    fEvBuffer = new TTree("EventBuffer", "Temporary buffer for mini events");
 326    fEvBuffer->Branch("events", "AliRsnMiniEvent", &fMiniEvent);
 327
 328    // create one histogram per each stored definition (event histograms)
 329    Int_t i, ndef = fHistograms.GetEntries();
 330    AliRsnMiniOutput *def = 0x0;
 331    for (i = 0; i < ndef; i++) {
 332       def = (AliRsnMiniOutput *)fHistograms[i];
 333       if (!def) continue;
 334       if (!def->Init(GetName(), fOutput)) {
 335          AliError(Form("Def '%s': failed initialization", def->GetName()));
 336          continue;
 337       }
 338    }
 339
 340    // post data for ALL output slots >0 here, to get at least an empty histogram
 341    PostData(1, fOutput);
 342 }
 343
 344 //__________________________________________________________________________________________________
 345 void AliRsnMiniAnalysisTask::UserExec(Option_t *)
 346 {
 347 //
 348 // Computation loop.
 349 // In this case, it checks if the event is acceptable, and eventually
 350 // creates the corresponding mini-event and stores it in the buffer.
 351 // The real histogram filling is done at the end, in "FinishTaskOutput".
 352 //
 353
 354    // increment event counter
 355    fEvNum++;
 356
 357    // check current event
 358    Char_t check = CheckCurrentEvent();
 359    if (!check) return;
 360
 361    // setup PID response
 362    AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 363    AliInputEventHandler *inputHandler = (AliInputEventHandler *)man->GetInputEventHandler();
 364    fRsnEvent.SetPIDResponse(inputHandler->GetPIDResponse());
 365
 366    // fill a mini-event from current
 367    // and skip this event if no tracks were accepted
 368    FillMiniEvent(check);
 369
 370    // fill MC based histograms on mothers,
 371    // which do need the original event
 372    if (fUseMC) {
 373       if (fRsnEvent.IsESD() && fMCEvent)
 374          FillTrueMotherESD(fMiniEvent);
 375       else if (fRsnEvent.IsAOD() && fRsnEvent.GetAODList())
 376          FillTrueMotherAOD(fMiniEvent);
 377    }
 378
 379    // if the event is not empty, store it
 380    if (fMiniEvent->IsEmpty()) {
 381       AliDebugClass(2, Form("Rejecting empty event #%d", fEvNum));
 382    } else {
 383       Int_t id = fEvBuffer->GetEntries();
 384       AliDebugClass(2, Form("Adding event #%d with ID = %d", fEvNum, id));
 385       fMiniEvent->ID() = id;
 386       fEvBuffer->Fill();
 387    }
 388
 389    // post data for computed stuff
 390    PostData(1, fOutput);
 391 }
 392
 393 //__________________________________________________________________________________________________
 394 void AliRsnMiniAnalysisTask::FinishTaskOutput()
 395 {
 396 //
 397 // This function is called at the end of the loop on available events,
 398 // and then the buffer will be full with all the corresponding mini-events,
 399 // each one containing all tracks selected by each of the available track cuts.
 400 // Here a loop is done on each of these events, and both single-event and mixing are computed
 401 //
 402
 403    // security code: reassign the buffer to the mini-event cursor
 404    fEvBuffer->SetBranchAddress("events", &fMiniEvent);
 405    TStopwatch timer;
 406    // prepare variables
 407    Int_t ievt, nEvents = (Int_t)fEvBuffer->GetEntries();
 408    Int_t idef, nDefs   = fHistograms.GetEntries();
 409    Int_t imix, iloop, ifill;
 410    AliRsnMiniOutput *def = 0x0;
 411    AliRsnMiniOutput::EComputation compType;
 412
 413    Int_t printNum = fMixPrintRefresh;
 414    if (printNum < 0) {
 415       if (nEvents>1e5) printNum=nEvents/100;
 416       else if (nEvents>1e4) printNum=nEvents/10;
 417       else printNum = 0;
 418    }
 419
 420    // loop on events, and for each one fill all outputs
 421    // using the appropriate procedure depending on its type
 422    // only mother-related histograms are filled in UserExec,
 423    // since they require direct access to MC event
 424    timer.Start();
 425    for (ievt = 0; ievt < nEvents; ievt++) {
 426       // get next entry
 427       fEvBuffer->GetEntry(ievt);
 428       if (printNum&&(ievt%printNum==0)) {
 429          AliInfo(Form("[%s] Std.Event %d/%d",GetName(), ievt,nEvents));
 430          timer.Stop(); timer.Print(); fflush(stdout); timer.Start(kFALSE);
 431       }
 432       // fill
 433       for (idef = 0; idef < nDefs; idef++) {
 434          def = (AliRsnMiniOutput *)fHistograms[idef];
 435          if (!def) continue;
 436          compType = def->GetComputation();
 437          // execute computation in the appropriate way
 438          switch (compType) {
 439             case AliRsnMiniOutput::kEventOnly:
 440                //AliDebugClass(1, Form("Event %d, def '%s': event-value histogram filling", ievt, def->GetName()));
 441                ifill = 1;
 442                def->FillEvent(fMiniEvent, &fValues);
 443                break;
 444             case AliRsnMiniOutput::kTruePair:
 445                //AliDebugClass(1, Form("Event %d, def '%s': true-pair histogram filling", ievt, def->GetName()));
 446                ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
 447                break;
 448             case AliRsnMiniOutput::kTrackPair:
 449                //AliDebugClass(1, Form("Event %d, def '%s': pair-value histogram filling", ievt, def->GetName()));
 450                ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
 451                break;
 452             case AliRsnMiniOutput::kTrackPairRotated1:
 453                //AliDebugClass(1, Form("Event %d, def '%s': rotated (1) background histogram filling", ievt, def->GetName()));
 454                ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
 455                break;
 456             case AliRsnMiniOutput::kTrackPairRotated2:
 457                //AliDebugClass(1, Form("Event %d, def '%s': rotated (2) background histogram filling", ievt, def->GetName()));
 458                ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
 459                break;
 460             default:
 461                // other kinds are processed elsewhere
 462                ifill = 0;
 463                AliDebugClass(2, Form("Computation = %d", (Int_t)compType));
 464          }
 465          // message
 466          AliDebugClass(1, Form("Event %6d: def = '%15s' -- fills = %5d", ievt, def->GetName(), ifill));
 467       }
 468    }
 469
 470    // if no mixing is required, stop here and post the output
 471    if (fNMix < 1) {
 472       AliDebugClass(2, "Stopping here, since no mixing is required");
 473       PostData(1, fOutput);
 474       return;
 475    }
 476
 477    // initialize mixing counter
 478    Int_t    nmatched[nEvents];
 479    TString *smatched = new TString[nEvents];
 480    for (ievt = 0; ievt < nEvents; ievt++) {
 481       smatched[ievt] = "|";
 482       nmatched[ievt] = 0;
 483    }
 484
 485
 486    AliInfo(Form("[%s] Std.Event %d/%d",GetName(), nEvents,nEvents));
 487    timer.Stop(); timer.Print(); timer.Start(); fflush(stdout);
 488
 489    // search for good matchings
 490    for (ievt = 0; ievt < nEvents; ievt++) {
 491       if (printNum&&(ievt%printNum==0)) {
 492          AliInfo(Form("[%s] EventMixing searching %d/%d",GetName(),ievt,nEvents));
 493          timer.Stop(); timer.Print(); timer.Start(kFALSE); fflush(stdout);
 494       }
 495       if (nmatched[ievt] >= fNMix) continue;
 496       fEvBuffer->GetEntry(ievt);
 497       AliRsnMiniEvent evMain(*fMiniEvent);
 498       for (iloop = 1; iloop < nEvents; iloop++) {
 499          imix = ievt + iloop;
 500          if (imix >= nEvents) imix -= nEvents;
 501          if (imix == ievt) continue;
 502          // text next entry
 503          fEvBuffer->GetEntry(imix);
 504          // skip if events are not matched
 505          if (!EventsMatch(&evMain, fMiniEvent)) continue;
 506          // check that the array of good matches for mixed does not already contain main event
 507          if (smatched[imix].Contains(Form("|%d|", ievt))) continue;
 508          // check that the found good events has not enough matches already
 509          if (nmatched[imix] >= fNMix) continue;
 510          // add new mixing candidate
 511          smatched[ievt].Append(Form("%d|", imix));
 512          nmatched[ievt]++;
 513          nmatched[imix]++;
 514          if (nmatched[ievt] >= fNMix) break;
 515       }
 516       AliDebugClass(1, Form("Matches for event %5d = %d [%s] (missing are declared above)", evMain.ID(), nmatched[ievt], smatched[ievt].Data()));
 517    }
 518
 519    AliInfo(Form("[%s] EventMixing searching %d/%d",GetName(),nEvents,nEvents));
 520    timer.Stop(); timer.Print(); fflush(stdout); timer.Start();
 521
 522    // perform mixing
 523    TObjArray *list = 0x0;
 524    TObjString *os = 0x0;
 525    for (ievt = 0; ievt < nEvents; ievt++) {
 526       if (printNum&&(ievt%printNum==0)) {
 527          AliInfo(Form("[%s] EventMixing %d/%d",GetName(),ievt,nEvents));
 528          timer.Stop(); timer.Print(); timer.Start(kFALSE); fflush(stdout);
 529       }
 530       ifill = 0;
 531       fEvBuffer->GetEntry(ievt);
 532       AliRsnMiniEvent evMain(*fMiniEvent);
 533       list = smatched[ievt].Tokenize("|");
 534       TObjArrayIter next(list);
 535       while ( (os = (TObjString *)next()) ) {
 536          imix = os->GetString().Atoi();
 537          fEvBuffer->GetEntry(imix);
 538          for (idef = 0; idef < nDefs; idef++) {
 539             def = (AliRsnMiniOutput *)fHistograms[idef];
 540             if (!def) continue;
 541             if (!def->IsTrackPairMix()) continue;
 542             ifill += def->FillPair(&evMain, fMiniEvent, &fValues, kTRUE);
 543             if (!def->IsSymmetric()) {
 544                AliDebugClass(2, "Reflecting non symmetric pair");
 545                ifill += def->FillPair(fMiniEvent, &evMain, &fValues, kFALSE);
 546             }
 547          }
 548       }
 549       delete list;
 550    }
 551
 552    delete [] smatched;
 553
 554    AliInfo(Form("[%s] EventMixing %d/%d",GetName(),nEvents,nEvents));
 555    timer.Stop(); timer.Print(); fflush(stdout);
 556
 557    /*
 558    OLD
 559    ifill = 0;
 560    for (iloop = 1; iloop < nEvents; iloop++) {
 561       imix = ievt + iloop;
 562       // restart from beginning if reached last event
 563       if (imix >= nEvents) imix -= nEvents;
 564       // avoid to mix an event with itself
 565       if (imix == ievt) continue;
 566       // skip all events already mixed enough times
 567       if (fNMixed[ievt] >= fNMix) break;
 568       if (fNMixed[imix] >= fNMix) continue;
 569       fEvBuffer->GetEntry(imix);
 570       // skip if events are not matched
 571       if (TMath::Abs(evMain.Vz()    - fMiniEvent->Vz()   ) > fMaxDiffVz   ) continue;
 572       if (TMath::Abs(evMain.Mult()  - fMiniEvent->Mult() ) > fMaxDiffMult ) continue;
 573       if (TMath::Abs(evMain.Angle() - fMiniEvent->Angle()) > fMaxDiffAngle) continue;
 574       // found a match: increment counter for both events
 575       AliDebugClass(1, Form("Event %d, def '%s': event mixing (%d with %d)", ievt, def->GetName(), ievt, imix));
 576       fNMixed[ievt]++;
 577       fNMixed[imix]++;
 578       // process mixing
 579       ifill += def->FillPair(&evMain, fMiniEvent, &fValues);
 580       // stop if mixed enough times
 581       if (fNMixed[ievt] >= fNMix) break;
 582    }
 583    break;
 584    // print number of mixings done with each event
 585    for (ievt = 0; ievt < nEvents; ievt++) {
 586       AliDebugClass(2, Form("Event %6d: mixed %2d times", ievt, fNMixed[ievt]));
 587    }
 588    */
 589
 590    // post computed data
 591    PostData(1, fOutput);
 592 }
 593
 594 //__________________________________________________________________________________________________
 595 void AliRsnMiniAnalysisTask::Terminate(Option_t *)
 596 {
 597 //
 598 // Draw result to screen, or perform fitting, normalizations
 599 // Called once at the end of the query
 600 //
 601
 602    fOutput = dynamic_cast<TList *>(GetOutputData(1));
 603    if (!fOutput) {
 604       AliError("Could not retrieve TList fOutput");
 605       return;
 606    }
 607 }
 608
 609 //__________________________________________________________________________________________________
 610 Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent()
 611 {
 612 //
 613 // This method checks if current event is OK for analysis.
 614 // In case it is, the pointers of the local AliRsnEvent data member
 615 // will point to it, in order to allow cut checking, otherwise the
 616 // function exits with a failure message.
 617 // ---
 618 // ESD events must pass the physics selection, AOD are supposed to do.
 619 // ---
 620 // While checking the event, a histogram is filled to count the number
 621 // of CINT1B, V0AND and CANDLE events, which are needed for normalization
 622 // ---
 623 // Return values can be:
 624 //    -- 'E' if the event is accepted and is ESD
 625 //    -- 'A' if the event is accepted and is AOD
 626 //    --  0  if the event is not accepted
 627 //
 628
 629    // string to sum messages
 630    TString msg("");
 631
 632    // check input type
 633    // exit points are provided in all cases an event is bad
 634    // if this block is passed, an event can be rejected only
 635    // if it does not pass the set of event cuts defined in the task
 636    Char_t output = 0;
 637    Bool_t isSelected;
 638    if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
 639       // type ESD
 640       output = 'E';
 641       // ESD specific check: Physics Selection
 642       // --> if this is failed, the event is rejected
 643       isSelected = (((AliInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
 644
 645       if (!isSelected) {
 646          AliDebugClass(2, "Event does not pass physics selections");
 647          fRsnEvent.SetRef(0x0);
 648          fRsnEvent.SetRefMC(0x0);
 649          return 0;
 650       }
 651       // set reference to input
 652       fRsnEvent.SetRef(fInputEvent);
 653       // add MC if requested and available
 654       if (fUseMC) {
 655          if (fMCEvent)
 656             fRsnEvent.SetRefMC(fMCEvent);
 657          else {
 658             AliWarning("MC event requested but not available");
 659             fRsnEvent.SetRefMC(0x0);
 660          }
 661       }
 662    } else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
 663       // type AOD
 664       output = 'A';
 665       // set reference to input
 666       fRsnEvent.SetRef(fInputEvent);
 667       // add MC if requested and available (it is in the same object)
 668       if (fUseMC) {
 669          fRsnEvent.SetRefMC(fInputEvent);
 670          if (!fRsnEvent.GetAODList()) {
 671             AliWarning("MC event requested but not available");
 672             fRsnEvent.SetRefMC(0x0);
 673          }
 674       }
 675    } else {
 676       AliError(Form("Bad input event class: %s", fInputEvent->ClassName()));
 677       // reset pointers in local AliRsnEvent object
 678       fRsnEvent.SetRef(0x0);
 679       fRsnEvent.SetRefMC(0x0);
 680       return 0;
 681    }
 682
 683    // fill counter of accepted events
 684    fHEventStat->Fill(0.1);
 685
 686    // check if it is V0AND
 687    // --> uses a cast to AliESDEvent even if the input is an AliAODEvent
 688    Bool_t v0A = fTriggerAna->IsOfflineTriggerFired((AliESDEvent *)fInputEvent, AliTriggerAnalysis::kV0A);
 689    Bool_t v0C = fTriggerAna->IsOfflineTriggerFired((AliESDEvent *)fInputEvent, AliTriggerAnalysis::kV0C);
 690    if (v0A && v0C) {
 691       msg += " -- VOAND = YES";
 692       fHEventStat->Fill(1.1);
 693    } else {
 694       msg += " -- VOAND = NO ";
 695    }
 696
 697    // check candle
 698    // --> requires at least one good quality track with Pt > 0.5 and |eta| <= 0.8
 699    Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
 700    Bool_t candle = kFALSE;
 701    for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
 702       AliVTrack   *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
 703       AliESDtrack *esdt  = dynamic_cast<AliESDtrack *>(track);
 704       AliAODTrack *aodt  = dynamic_cast<AliAODTrack *>(track);
 705       if (track->Pt() < 0.5) continue;
 706       if(TMath::Abs(track->Eta()) > 0.8) continue;
 707       if (esdt) if (!fESDtrackCuts->AcceptTrack(esdt)) continue;
 708       if (aodt) if (!aodt->TestFilterBit(5)) continue;
 709       candle = kTRUE;
 710       break;
 711    }
 712    if (candle) {
 713       msg += " -- CANDLE = YES";
 714       fHEventStat->Fill(2.1);
 715    } else {
 716       msg += " -- CANDLE = NO ";
 717    }
 718
 719    // if event cuts are defined, they are checked here
 720    // final decision on the event depends on this
 721    isSelected = kTRUE;
 722    if (fEventCuts) {
 723       if (!fEventCuts->IsSelected(&fRsnEvent)) {
 724          msg += " -- Local cuts = REJECTED";
 725          isSelected = kFALSE;
 726       } else {
 727          msg += " -- Local cuts = ACCEPTED";
 728          isSelected = kTRUE;
 729       }
 730    } else {
 731       msg += " -- Local cuts = NONE";
 732       isSelected = kTRUE;
 733    }
 734
 735    // if the above exit point is not taken, the event is accepted
 736    AliDebugClass(2, Form("Stats: %s", msg.Data()));
 737    if (isSelected) {
 738       fHEventStat->Fill(3.1);
 739       Double_t multi = ComputeCentrality((output == 'E'));
 740       fHAEventsVsMulti->Fill(multi);
 741       if(fHAEventVz) fHAEventVz->Fill(multi,fInputEvent->GetPrimaryVertex()->GetZ());
 742       if(fHAEventMultiCent) fHAEventMultiCent->Fill(multi,ComputeMultiplicity(output == 'E',fHAEventMultiCent->GetYaxis()->GetTitle()));
 743       if(fHAEventPlane) fHAEventPlane->Fill(multi,ComputeAngle());
 744       return output;
 745    } else {
 746       return 0;
 747    }
 748 }
 749
 750 //__________________________________________________________________________________________________
 751 void AliRsnMiniAnalysisTask::FillMiniEvent(Char_t evType)
 752 {
 753 //
 754 // Refresh cursor mini-event data member to fill with current event.
 755 // Returns the total number of tracks selected.
 756 //
 757
 758    // assign event-related values
 759    if (fMiniEvent) delete fMiniEvent;
 760    fMiniEvent = new AliRsnMiniEvent;
 761    fMiniEvent->Vz()    = fInputEvent->GetPrimaryVertex()->GetZ();
 762    fMiniEvent->Angle() = ComputeAngle();
 763    fMiniEvent->Mult()  = ComputeCentrality((evType == 'E'));
 764    AliDebugClass(2, Form("Event %d: type = %c -- vz = %f -- mult = %f -- angle = %f", fEvNum, evType, fMiniEvent->Vz(), fMiniEvent->Mult(), fMiniEvent->Angle()));
 765
 766    // loop on daughters and assign track-related values
 767    Int_t ic, ncuts = fTrackCuts.GetEntries();
 768    Int_t ip, npart = fRsnEvent.GetAbsoluteSum();
 769    Int_t npos = 0, nneg = 0, nneu = 0;
 770    AliRsnDaughter cursor;
 771    AliRsnMiniParticle miniParticle;
 772    for (ip = 0; ip < npart; ip++) {
 773       // point cursor to next particle
 774       fRsnEvent.SetDaughter(cursor, ip);
 775       // copy momentum and MC info if present
 776       miniParticle.CopyDaughter(&cursor);
 777       miniParticle.Index() = ip;
 778       // switch on the bits corresponding to passed cuts
 779       for (ic = 0; ic < ncuts; ic++) {
 780          AliRsnCutSet *cuts = (AliRsnCutSet *)fTrackCuts[ic];
 781          if (cuts->IsSelected(&cursor)) miniParticle.SetCutBit(ic);
 782       }
 783       // if a track passes at least one track cut, it is added to the pool
 784       if (miniParticle.CutBits()) {
 785          fMiniEvent->AddParticle(miniParticle);
 786          if (miniParticle.Charge() == '+') npos++;
 787          else if (miniParticle.Charge() == '-') nneg++;
 788          else nneu++;
 789       }
 790    }
 791
 792    // get number of accepted tracks
 793    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));
 794 }
 795
 796 //__________________________________________________________________________________________________
 797 Double_t AliRsnMiniAnalysisTask::ComputeAngle()
 798 {
 799 //
 800 // Get the plane angle
 801 //
 802
 803    AliEventplane *plane = 0x0;
 804
 805    if (fInputEvent->InheritsFrom(AliESDEvent::Class()))
 806       plane = fInputEvent->GetEventplane();
 807    else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
 808       AliAODEvent *aodEvent = (AliAODEvent *)fInputEvent;
 809       plane = aodEvent->GetHeader()->GetEventplaneP();
 810    }
 811
 812    if (plane)
 813       return plane->GetEventplane("Q");
 814    else {
 815       AliWarning("No event plane defined");
 816       return 1E20;
 817    }
 818 }
 819
 820 //__________________________________________________________________________________________________
 821 Double_t AliRsnMiniAnalysisTask::ComputeCentrality(Bool_t isESD)
 822 {
 823 //
 824 // Computes event centrality/multiplicity according to the criterion defined
 825 // by two elements: (1) choice between multiplicity and centrality and
 826 // (2) the string defining what criterion must be used for specific computation.
 827 //
 828
 829    if (fUseCentrality) {
 830      if ((!fUseMC) && (fUseCentralityPatchPbPb2011)) {
 831        return ApplyCentralityPatchPbPb2011();//
 832     }
 833      if ((!fUseMC) && (!isESD) && (fUseAOD049CentralityPatch)) {
 834        return ApplyCentralityPatchAOD049();
 835      } else {
 836        AliCentrality *centrality = fInputEvent->GetCentrality();
 837          if (!centrality) {
 838            AliError("Cannot compute centrality!");
 839            return -1.0;
 840          }
 841          return centrality->GetCentralityPercentile(fCentralityType.Data());
 842      }
 843    } else {
 844       if (!fCentralityType.CompareTo("TRACKS"))
 845          return fInputEvent->GetNumberOfTracks();
 846       else if (!fCentralityType.CompareTo("QUALITY"))
 847          if (isESD)
 848             return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent *)fInputEvent, kTRUE);
 849          else {
 850             Double_t count = 0.;
 851             Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
 852             for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
 853                AliVTrack   *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
 854                AliAODTrack *aodt  = dynamic_cast<AliAODTrack *>(track);
 855                if (!aodt) continue;
 856                if (!aodt->TestFilterBit(5)) continue;
 857                count++;
 858             }
 859             return count;
 860          }
 861       else if (!fCentralityType.CompareTo("TRACKLETS")) {
 862          if (isESD) {
 863             const AliMultiplicity *mult = ((AliESDEvent *)fInputEvent)->GetMultiplicity();
 864             Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
 865             for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
 866             return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ());
 867          } else {
 868             AliWarning("Cannot compute multiplicity with SPD tracklets from AOD");
 869             return 1E20;
 870          }
 871       } else {
 872          AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", fCentralityType.Data()));
 873          return -1.0;
 874       }
 875    }
 876 }
 877
 878 //__________________________________________________________________________________________________
 879 Double_t AliRsnMiniAnalysisTask::ComputeMultiplicity(Bool_t isESD,TString type)
 880 {
 881 //
 882 // Computes event multiplicity according to the string defining
 883 // what criterion must be used for specific computation.
 884 //
 885
 886    type.ToUpper();
 887
 888    if (!type.CompareTo("TRACKS"))
 889       return fInputEvent->GetNumberOfTracks();
 890    else if (!type.CompareTo("QUALITY"))
 891       if (isESD)
 892          return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent *)fInputEvent, kTRUE);
 893       else {
 894          Double_t count = 0.;
 895          Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
 896          for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
 897             AliVTrack   *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
 898             AliAODTrack *aodt  = dynamic_cast<AliAODTrack *>(track);
 899             if (!aodt) continue;
 900             if (!aodt->TestFilterBit(5)) continue;
 901             count++;
 902          }
 903          return count;
 904       }
 905    else if (!type.CompareTo("TRACKLETS")) {
 906       if (isESD) {
 907          const AliMultiplicity *mult = ((AliESDEvent *)fInputEvent)->GetMultiplicity();
 908          Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
 909          for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
 910          return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ());
 911       } else {
 912          AliWarning("Cannot compute multiplicity with SPD tracklets from AOD");
 913          return 1E20;
 914       }
 915    } else {
 916       AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", type.Data()));
 917       return -1.0;
 918    }
 919
 920    return 1E20;
 921 }
 922
 923 //__________________________________________________________________________________________________
 924 void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
 925 {
 926 //
 927 // Fills the histograms with true mother (ESD version)
 928 //
 929
 930    Bool_t okMatch;
 931    Int_t id, ndef = fHistograms.GetEntries();
 932    Int_t ip, label1, label2, npart = fMCEvent->GetNumberOfTracks();
 933    static AliRsnMiniPair miniPair;
 934    AliMCParticle *daughter1, *daughter2;
 935    TLorentzVector p1, p2;
 936    AliRsnMiniOutput *def = 0x0;
 937
 938    for (id = 0; id < ndef; id++) {
 939       def = (AliRsnMiniOutput *)fHistograms[id];
 940       if (!def) continue;
 941       if (!def->IsMother()) continue;
 942       for (ip = 0; ip < npart; ip++) {
 943          AliMCParticle *part = (AliMCParticle *)fMCEvent->GetTrack(ip);
 944          //get mother pdg code
 945          if (part->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
 946          // check that daughters match expected species
 947          if (part->Particle()->GetNDaughters() < 2) continue;
 948          if (fMaxNDaughters > 0 && part->Particle()->GetNDaughters() > fMaxNDaughters) continue;
 949          label1 = part->Particle()->GetDaughter(0);
 950          label2 = part->Particle()->GetDaughter(1);
 951          daughter1 = (AliMCParticle *)fMCEvent->GetTrack(label1);
 952          daughter2 = (AliMCParticle *)fMCEvent->GetTrack(label2);
 953          okMatch = kFALSE;
 954          if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(1)) {
 955             okMatch = kTRUE;
 956             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
 957             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
 958          } else if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(0)) {
 959             okMatch = kTRUE;
 960             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
 961             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
 962          }
 963          if (!okMatch) continue;
 964          if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
 965                                 (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
 966                                 (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
 967          if(fCheckFeedDown){
 968                 Int_t pdgGranma = 0;
 969                 Int_t mother = 0;
 970                 mother = part->GetMother();
 971                 Int_t istep = 0;
 972                 Int_t abspdgGranma =0;
 973                 Bool_t isFromB=kFALSE;
 974                 Bool_t isQuarkFound=kFALSE;
 975                 while (mother >=0 ){
 976                         istep++;
 977                         AliDebug(2,Form("mother at step %d = %d", istep, mother));
 978                         AliMCParticle* mcGranma = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(mother));
 979                         if (mcGranma){
 980                                 pdgGranma = mcGranma->PdgCode();
 981                                 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
 982                                 abspdgGranma = TMath::Abs(pdgGranma);
 983                                 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 984                                   isFromB=kTRUE;
 985                                 }
 986                                 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 987                                 mother = mcGranma->GetMother();
 988                         }else{
 989                                 AliError("Failed casting the mother particle!");
 990                                 break;
 991                         }
 992                 }
 993                 if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
 994                 if(isFromB){
 995                   if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
 996                 }
 997                 else{
 998                   if (fKeepDfromBOnly) pdgGranma = -999;
 999
1000                 if (pdgGranma == -99999){
1001                         AliDebug(2,"This particle does not have a quark in his genealogy\n");
1002                         continue;
1003                 }
1004                 if (pdgGranma == -9999){
1005                         AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");
1006                         continue;
1007                 }
1008
1009                 if (pdgGranma == -999){
1010                         AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");
1011                         continue;
1012                 }
1013
1014               }
1015          }
1016          // assign momenta to computation object
1017          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
1018          miniPair.FillRef(def->GetMotherMass());
1019          // do computations
1020          def->FillMother(&miniPair, miniEvent, &fValues);
1021       }
1022    }
1023 }
1024
1025 //__________________________________________________________________________________________________
1026 void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
1027 {
1028 //
1029 // Fills the histograms with true mother (AOD version)
1030 //
1031
1032    Bool_t okMatch;
1033    TClonesArray *list = fRsnEvent.GetAODList();
1034    Int_t id, ndef = fHistograms.GetEntries();
1035    Int_t ip, label1, label2, npart = list->GetEntries();
1036    static AliRsnMiniPair miniPair;
1037    AliAODMCParticle *daughter1, *daughter2;
1038    TLorentzVector p1, p2;
1039    AliRsnMiniOutput *def = 0x0;
1040
1041    for (id = 0; id < ndef; id++) {
1042       def = (AliRsnMiniOutput *)fHistograms[id];
1043       if (!def) continue;
1044       if (!def->IsMother()) continue;
1045       for (ip = 0; ip < npart; ip++) {
1046          AliAODMCParticle *part = (AliAODMCParticle *)list->At(ip);
1047          if (part->GetPdgCode() != def->GetMotherPDG()) continue;
1048          // check that daughters match expected species
1049          if (part->GetNDaughters() < 2) continue;
1050          if (fMaxNDaughters > 0 && part->GetNDaughters() > fMaxNDaughters) continue;
1051          label1 = part->GetDaughter(0);
1052          label2 = part->GetDaughter(1);
1053          daughter1 = (AliAODMCParticle *)list->At(label1);
1054          daughter2 = (AliAODMCParticle *)list->At(label2);
1055          okMatch = kFALSE;
1056          if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(1)) {
1057             okMatch = kTRUE;
1058             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
1059             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
1060          } else if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(0)) {
1061             okMatch = kTRUE;
1062             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
1063             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
1064          }
1065          if (!okMatch) continue;
1066          if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
1067                                 (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
1068                                 (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
1069          if(fCheckFeedDown){
1070                 Int_t pdgGranma = 0;
1071                 Int_t mother = 0;
1072                 mother = part->GetMother();
1073                 Int_t istep = 0;
1074                 Int_t abspdgGranma =0;
1075                 Bool_t isFromB=kFALSE;
1076                 Bool_t isQuarkFound=kFALSE;
1077                 while (mother >=0 ){
1078                         istep++;
1079                         AliDebug(2,Form("mother at step %d = %d", istep, mother));
1080                         AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(list->At(mother));
1081                         if (mcGranma){
1082                                 pdgGranma = mcGranma->GetPdgCode();
1083                                 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1084                                 abspdgGranma = TMath::Abs(pdgGranma);
1085                                 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1086                                   isFromB=kTRUE;
1087                                 }
1088                                 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1089                                 mother = mcGranma->GetMother();
1090                         }else{
1091                                 AliError("Failed casting the mother particle!");
1092                                 break;
1093                         }
1094                 }
1095                 if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
1096                 if(isFromB){
1097                   if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
1098                 }
1099                 else{
1100                   if (fKeepDfromBOnly) pdgGranma = -999;
1101                   }
1102
1103                 if (pdgGranma == -99999){
1104                         AliDebug(2,"This particle does not have a quark in his genealogy\n");
1105                         continue;
1106                 }
1107                 if (pdgGranma == -9999){
1108                         AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");
1109                         continue;
1110                 }
1111
1112                 if (pdgGranma == -999){
1113                         AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");
1114                         continue;
1115                 }
1116          }
1117          // assign momenta to computation object
1118          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
1119          miniPair.FillRef(def->GetMotherMass());
1120          // do computations
1121          def->FillMother(&miniPair, miniEvent, &fValues);
1122       }
1123    }
1124 }
1125 //___________________________________________________________
1126 void AliRsnMiniAnalysisTask::SetDselection(UShort_t originDselection)
1127 {
1128         // setting the way the D0 will be selected
1129         // 0 --> only from c quarks
1130         // 1 --> only from b quarks
1131         // 2 --> from both c quarks and b quarks
1132
1133         fOriginDselection = originDselection;
1134
1135         if (fOriginDselection == 0) {
1136                 fKeepDfromB = kFALSE;
1137                 fKeepDfromBOnly = kFALSE;
1138         }
1139
1140         if (fOriginDselection == 1) {
1141                 fKeepDfromB = kTRUE;
1142                 fKeepDfromBOnly = kTRUE;
1143         }
1144
1145         if (fOriginDselection == 2) {
1146                 fKeepDfromB = kTRUE;
1147                 fKeepDfromBOnly = kFALSE;
1148         }
1149
1150         return;
1151 }
1152 //__________________________________________________________________________________________________
1153 Bool_t AliRsnMiniAnalysisTask::EventsMatch(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2)
1154 {
1155 //
1156 // Check if two events are compatible.
1157 // If the mixing is continuous, this is true if differences in vz, mult and angle are smaller than
1158 // the specified values.
1159 // If the mixing is binned, this is true if the events are in the same bin.
1160 //
1161
1162    if (!event1 || !event2) return kFALSE;
1163    Int_t ivz1, ivz2, imult1, imult2, iangle1, iangle2;
1164    Double_t dv, dm, da;
1165
1166    if (fContinuousMix) {
1167       dv = TMath::Abs(event1->Vz()    - event2->Vz()   );
1168       dm = TMath::Abs(event1->Mult()  - event2->Mult() );
1169       da = TMath::Abs(event1->Angle() - event2->Angle());
1170       if (dv > fMaxDiffVz) {
1171          //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Vz = %f", event1->ID(), event2->ID(), dv));
1172          return kFALSE;
1173       }
1174       if (dm > fMaxDiffMult ) {
1175          //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Mult = %f", event1->ID(), event2->ID(), dm));
1176          return kFALSE;
1177       }
1178       if (da > fMaxDiffAngle) {
1179          //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Angle = %f", event1->ID(), event2->ID(), da));
1180          return kFALSE;
1181       }
1182       return kTRUE;
1183    } else {
1184       ivz1 = (Int_t)(event1->Vz() / fMaxDiffVz);
1185       ivz2 = (Int_t)(event2->Vz() / fMaxDiffVz);
1186       imult1 = (Int_t)(event1->Mult() / fMaxDiffMult);
1187       imult2 = (Int_t)(event2->Mult() / fMaxDiffMult);
1188       iangle1 = (Int_t)(event1->Angle() / fMaxDiffAngle);
1189       iangle2 = (Int_t)(event2->Angle() / fMaxDiffAngle);
1190       if (ivz1 != ivz2) return kFALSE;
1191       if (imult1 != imult2) return kFALSE;
1192       if (iangle1 != iangle2) return kFALSE;
1193       return kTRUE;
1194    }
1195 }
1196
1197 //---------------------------------------------------------------------
1198 Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchPbPb2011(){
1199   //This part rejects randomly events such that the centrality gets flat for LHC11h Pb-Pb data
1200   //for 0-5% and 10-20% centrality bin
1201
1202   if (fCentralityType!="V0M") {
1203     AliWarning("Wrong value (not centrality from V0).");
1204     return -999.0;
1205   }
1206
1207   AliCentrality *centrality = fInputEvent->GetCentrality();
1208   if (!centrality) {
1209     AliWarning("Cannot get centrality from AOD event.");
1210     return -999.0;
1211   }
1212
1213
1214   Double_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
1215   Double_t rnd_hc, testf, ff, N1, N2;
1216
1217   if(fUseCentralityPatchPbPb2011==5){
1218     N1 = 1.9404e+06;
1219     N2 = 1.56435e+06;
1220     ff = 5.04167e+06 - 1.49885e+06*cent + 2.35998e+05*cent*cent -1.22873e+04*cent*cent*cent;
1221   }
1222
1223   if(fUseCentralityPatchPbPb2011==1020){
1224     N1 = 1.56435e+06;
1225     N2 = 4.20e+05;
1226     ff = 1.68062e+08 - 5.19673e+07*cent + 6.4068e+06*cent*cent + 6.4068e+06*cent*cent*cent - 392687*cent*cent*cent*cent - 145.07*cent*cent*cent*cent*cent;
1227   }
1228
1229   testf = ( N2 + (N1-ff) ) / N1;
1230   rnd_hc = gRandom->Rndm();
1231
1232   //AliDebugClass(1, Form("Flat Centrality %d", fUseCentralityPatchPbPb2011));
1233
1234   if (rnd_hc < 0 || rnd_hc > 1 )
1235     {
1236       AliWarning("Wrong Random number generated");
1237       return -999.0;
1238     }
1239
1240   if (rnd_hc < testf)
1241     return cent;
1242   else
1243     return -999.0;
1244 }
1245 //---------------------------------------------------------------------
1246 Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchAOD049()
1247 {
1248    //
1249    //Apply centrality patch for AOD049 outliers
1250    //
1251    if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
1252       AliWarning("Requested patch for AOD049 for ESD. ");
1253       return -999.0;
1254    }
1255
1256    if (fCentralityType!="V0M") {
1257       AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
1258       return -999.0;
1259    }
1260
1261    AliCentrality *centrality = fInputEvent->GetCentrality();
1262    if (!centrality) {
1263       AliWarning("Cannot get centrality from AOD event.");
1264       return -999.0;
1265    }
1266
1267    Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
1268    /*
1269    Bool_t isSelRun = kFALSE;
1270    Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
1271    if(cent<0){
1272      Int_t quality = centrality->GetQuality();
1273      if(quality<=1){
1274        cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1275      } else {
1276        Int_t runnum=aodEvent->GetRunNumber();
1277        for(Int_t ir=0;ir<5;ir++){
1278    if(runnum==selRun[ir]){
1279      isSelRun=kTRUE;
1280      break;
1281    }
1282        }
1283        if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1284      }
1285    }
1286    */
1287
1288    if(cent>=0.0) {
1289       Float_t v0 = 0.0;
1290       AliAODEvent *aodEvent = (AliAODEvent *)fInputEvent;
1291       AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
1292       v0+=aodV0->GetMTotV0A();
1293       v0+=aodV0->GetMTotV0C();
1294       if ( (cent==0) && (v0<19500) ) {
1295          AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
1296          return -999.0;
1297       }
1298       Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1299       Float_t val = 1.30552 +  0.147931 * v0;
1300
1301       Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
1302                                120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
1303                                92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
1304                                68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
1305                                51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
1306                                37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
1307                                26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
1308                                19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
1309                                12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
1310                                13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
1311                               };
1312
1313       if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )  {
1314          AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
1315          return -999.0;
1316       }
1317    } else {
1318       //force it to be -999. whatever the negative value was
1319       cent = -999.;
1320    }
1321    return cent;
1322 }
1323
1324 //----------------------------------------------------------------------------------
1325 void AliRsnMiniAnalysisTask::SetEventQAHist(TString type,TH2F *histo)
1326 {
1327    if(!histo) {
1328       AliWarning(Form("event QA histogram pointer not defined for slot %s",type.Data()));
1329       return;
1330    }
1331
1332    type.ToLower();
1333
1334    if(!type.CompareTo("vz")) fHAEventVz = histo;
1335    else if(!type.CompareTo("multicent")) {
1336       TString mtype(histo->GetYaxis()->GetTitle());
1337       mtype.ToUpper();
1338       if(mtype.CompareTo("QUALITY") && mtype.CompareTo("TRACKS") && mtype.CompareTo("TRACKLETS")) {
1339          AliWarning(Form("multiplicity vs. centrality histogram y-axis %s unknown, setting to TRACKS",mtype.Data()));
1340          histo->GetYaxis()->SetTitle("TRACKS");
1341       }
1342       fHAEventMultiCent = histo;
1343    }
1344    else if(!type.CompareTo("eventplane")) fHAEventPlane = histo;
1345    else AliWarning(Form("event QA histogram slot %s undefined",type.Data()));
1346
1347    return;
1348 }
1349
1350 //----------------------------------------------------------------------------------
1351 Int_t AliRsnMiniAnalysisTask::CreateValue(AliRsnMiniValue::EType type, Bool_t useMC)
1352 {
1353 //
1354 // Create a new value in the task,
1355 // and returns its ID, which is needed for setting up histograms.
1356 // If that value was already initialized, returns its ID and does not recreate it.
1357 //
1358
1359    Int_t valID = ValueID(type, useMC);
1360    if (valID >= 0 && valID < fValues.GetEntries()) {
1361       AliInfo(Form("Value '%s' is already created in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
1362    } else {
1363       valID = fValues.GetEntries();
1364       AliInfo(Form("Creating value '%s' in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
1365       new (fValues[valID]) AliRsnMiniValue(type, useMC);
1366    }
1367
1368    return valID;
1369 }
1370
1371 //----------------------------------------------------------------------------------
1372 Int_t AliRsnMiniAnalysisTask::ValueID(AliRsnMiniValue::EType type, Bool_t useMC)
1373 {
1374 //
1375 // Searches if a value computation is initialized
1376 //
1377
1378    const char *name = AliRsnMiniValue::ValueName(type, useMC);
1379    TObject *obj = fValues.FindObject(name);
1380    if (obj)
1381       return fValues.IndexOf(obj);
1382    else
1383       return -1;
1384 }
1385
1386 //----------------------------------------------------------------------------------
1387 AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, AliRsnMiniOutput::EOutputType type, AliRsnMiniOutput::EComputation src)
1388 {
1389 //
1390 // Create a new histogram definition in the task,
1391 // which is then returned to the user for its configuration
1392 //
1393
1394    Int_t n = fHistograms.GetEntries();
1395    AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, type, src);
1396
1397    return newDef;
1398 }
1399
1400 //----------------------------------------------------------------------------------
1401 AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, const char *outType, const char *compType)
1402 {
1403 //
1404 // Create a new histogram definition in the task,
1405 // which is then returned to the user for its configuration
1406 //
1407
1408    Int_t n = fHistograms.GetEntries();
1409    AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, outType, compType);
1410
1411    return newDef;
1412 }