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
 921 //__________________________________________________________________________________________________
 922 void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
 923 {
 924 //
 925 // Fills the histograms with true mother (ESD version)
 926 //
 927
 928    Bool_t okMatch;
 929    Int_t id, ndef = fHistograms.GetEntries();
 930    Int_t ip, label1, label2, npart = fMCEvent->GetNumberOfTracks();
 931    static AliRsnMiniPair miniPair;
 932    AliMCParticle *daughter1, *daughter2;
 933    TLorentzVector p1, p2;
 934    AliRsnMiniOutput *def = 0x0;
 935
 936    for (id = 0; id < ndef; id++) {
 937       def = (AliRsnMiniOutput *)fHistograms[id];
 938       if (!def) continue;
 939       if (!def->IsMother()) continue;
 940       for (ip = 0; ip < npart; ip++) {
 941          AliMCParticle *part = (AliMCParticle *)fMCEvent->GetTrack(ip);
 942          //get mother pdg code
 943          if (part->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
 944          // check that daughters match expected species
 945          if (part->Particle()->GetNDaughters() < 2) continue;
 946          if (fMaxNDaughters > 0 && part->Particle()->GetNDaughters() > fMaxNDaughters) continue;
 947          label1 = part->Particle()->GetDaughter(0);
 948          label2 = part->Particle()->GetDaughter(1);
 949          daughter1 = (AliMCParticle *)fMCEvent->GetTrack(label1);
 950          daughter2 = (AliMCParticle *)fMCEvent->GetTrack(label2);
 951          okMatch = kFALSE;
 952          if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(1)) {
 953             okMatch = kTRUE;
 954             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
 955             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
 956          } else if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(0)) {
 957             okMatch = kTRUE;
 958             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
 959             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
 960          }
 961          if (!okMatch) continue;
 962          if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
 963                                 (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
 964                                 (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
 965          if(fCheckFeedDown){
 966                 Int_t pdgGranma = 0;
 967                 Int_t mother = 0;
 968                 mother = part->GetMother();
 969                 Int_t istep = 0;
 970                 Int_t abspdgGranma =0;
 971                 Bool_t isFromB=kFALSE;
 972                 Bool_t isQuarkFound=kFALSE;
 973                 while (mother >=0 ){
 974                         istep++;
 975                         AliDebug(2,Form("mother at step %d = %d", istep, mother));
 976                         AliMCParticle* mcGranma = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(mother));
 977                         if (mcGranma){
 978                                 pdgGranma = mcGranma->PdgCode();
 979                                 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
 980                                 abspdgGranma = TMath::Abs(pdgGranma);
 981                                 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 982                                   isFromB=kTRUE;
 983                                 }
 984                                 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 985                                 mother = mcGranma->GetMother();
 986                         }else{
 987                                 AliError("Failed casting the mother particle!");
 988                                 break;
 989                         }
 990                 }
 991                 if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
 992                 if(isFromB){
 993                   if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
 994                 }
 995                 else{
 996                   if (fKeepDfromBOnly) pdgGranma = -999;
 997
 998                 if (pdgGranma == -99999){
 999                         AliDebug(2,"This particle does not have a quark in his genealogy\n");
1000                         continue;
1001                 }
1002                 if (pdgGranma == -9999){
1003                         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");
1004                         continue;
1005                 }
1006
1007                 if (pdgGranma == -999){
1008                         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");
1009                         continue;
1010                 }
1011
1012               }
1013          }
1014          // assign momenta to computation object
1015          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
1016          miniPair.FillRef(def->GetMotherMass());
1017          // do computations
1018          def->FillMother(&miniPair, miniEvent, &fValues);
1019       }
1020    }
1021 }
1022
1023 //__________________________________________________________________________________________________
1024 void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
1025 {
1026 //
1027 // Fills the histograms with true mother (AOD version)
1028 //
1029
1030    Bool_t okMatch;
1031    TClonesArray *list = fRsnEvent.GetAODList();
1032    Int_t id, ndef = fHistograms.GetEntries();
1033    Int_t ip, label1, label2, npart = list->GetEntries();
1034    static AliRsnMiniPair miniPair;
1035    AliAODMCParticle *daughter1, *daughter2;
1036    TLorentzVector p1, p2;
1037    AliRsnMiniOutput *def = 0x0;
1038
1039    for (id = 0; id < ndef; id++) {
1040       def = (AliRsnMiniOutput *)fHistograms[id];
1041       if (!def) continue;
1042       if (!def->IsMother()) continue;
1043       for (ip = 0; ip < npart; ip++) {
1044          AliAODMCParticle *part = (AliAODMCParticle *)list->At(ip);
1045          if (part->GetPdgCode() != def->GetMotherPDG()) continue;
1046          // check that daughters match expected species
1047          if (part->GetNDaughters() < 2) continue;
1048          if (fMaxNDaughters > 0 && part->GetNDaughters() > fMaxNDaughters) continue;
1049          label1 = part->GetDaughter(0);
1050          label2 = part->GetDaughter(1);
1051          daughter1 = (AliAODMCParticle *)list->At(label1);
1052          daughter2 = (AliAODMCParticle *)list->At(label2);
1053          okMatch = kFALSE;
1054          if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(1)) {
1055             okMatch = kTRUE;
1056             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0));
1057             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1));
1058          } else if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(0)) {
1059             okMatch = kTRUE;
1060             p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1));
1061             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
1062          }
1063          if (!okMatch) continue;
1064          if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
1065                                 (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
1066                                 (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
1067          if(fCheckFeedDown){
1068                 Int_t pdgGranma = 0;
1069                 Int_t mother = 0;
1070                 mother = part->GetMother();
1071                 Int_t istep = 0;
1072                 Int_t abspdgGranma =0;
1073                 Bool_t isFromB=kFALSE;
1074                 Bool_t isQuarkFound=kFALSE;
1075                 while (mother >=0 ){
1076                         istep++;
1077                         AliDebug(2,Form("mother at step %d = %d", istep, mother));
1078                         AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(list->At(mother));
1079                         if (mcGranma){
1080                                 pdgGranma = mcGranma->GetPdgCode();
1081                                 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1082                                 abspdgGranma = TMath::Abs(pdgGranma);
1083                                 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1084                                   isFromB=kTRUE;
1085                                 }
1086                                 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1087                                 mother = mcGranma->GetMother();
1088                         }else{
1089                                 AliError("Failed casting the mother particle!");
1090                                 break;
1091                         }
1092                 }
1093                 if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
1094                 if(isFromB){
1095                   if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
1096                 }
1097                 else{
1098                   if (fKeepDfromBOnly) pdgGranma = -999;
1099                   }
1100
1101                 if (pdgGranma == -99999){
1102                         AliDebug(2,"This particle does not have a quark in his genealogy\n");
1103                         continue;
1104                 }
1105                 if (pdgGranma == -9999){
1106                         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");
1107                         continue;
1108                 }
1109
1110                 if (pdgGranma == -999){
1111                         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");
1112                         continue;
1113                 }
1114          }
1115          // assign momenta to computation object
1116          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
1117          miniPair.FillRef(def->GetMotherMass());
1118          // do computations
1119          def->FillMother(&miniPair, miniEvent, &fValues);
1120       }
1121    }
1122 }
1123 //___________________________________________________________
1124 void AliRsnMiniAnalysisTask::SetDselection(UShort_t originDselection)
1125 {
1126         // setting the way the D0 will be selected
1127         // 0 --> only from c quarks
1128         // 1 --> only from b quarks
1129         // 2 --> from both c quarks and b quarks
1130
1131         fOriginDselection = originDselection;
1132
1133         if (fOriginDselection == 0) {
1134                 fKeepDfromB = kFALSE;
1135                 fKeepDfromBOnly = kFALSE;
1136         }
1137
1138         if (fOriginDselection == 1) {
1139                 fKeepDfromB = kTRUE;
1140                 fKeepDfromBOnly = kTRUE;
1141         }
1142
1143         if (fOriginDselection == 2) {
1144                 fKeepDfromB = kTRUE;
1145                 fKeepDfromBOnly = kFALSE;
1146         }
1147
1148         return;
1149 }
1150 //__________________________________________________________________________________________________
1151 Bool_t AliRsnMiniAnalysisTask::EventsMatch(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2)
1152 {
1153 //
1154 // Check if two events are compatible.
1155 // If the mixing is continuous, this is true if differences in vz, mult and angle are smaller than
1156 // the specified values.
1157 // If the mixing is binned, this is true if the events are in the same bin.
1158 //
1159
1160    if (!event1 || !event2) return kFALSE;
1161    Int_t ivz1, ivz2, imult1, imult2, iangle1, iangle2;
1162    Double_t dv, dm, da;
1163
1164    if (fContinuousMix) {
1165       dv = TMath::Abs(event1->Vz()    - event2->Vz()   );
1166       dm = TMath::Abs(event1->Mult()  - event2->Mult() );
1167       da = TMath::Abs(event1->Angle() - event2->Angle());
1168       if (dv > fMaxDiffVz) {
1169          //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Vz = %f", event1->ID(), event2->ID(), dv));
1170          return kFALSE;
1171       }
1172       if (dm > fMaxDiffMult ) {
1173          //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Mult = %f", event1->ID(), event2->ID(), dm));
1174          return kFALSE;
1175       }
1176       if (da > fMaxDiffAngle) {
1177          //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Angle = %f", event1->ID(), event2->ID(), da));
1178          return kFALSE;
1179       }
1180       return kTRUE;
1181    } else {
1182       ivz1 = (Int_t)(event1->Vz() / fMaxDiffVz);
1183       ivz2 = (Int_t)(event2->Vz() / fMaxDiffVz);
1184       imult1 = (Int_t)(event1->Mult() / fMaxDiffMult);
1185       imult2 = (Int_t)(event2->Mult() / fMaxDiffMult);
1186       iangle1 = (Int_t)(event1->Angle() / fMaxDiffAngle);
1187       iangle2 = (Int_t)(event2->Angle() / fMaxDiffAngle);
1188       if (ivz1 != ivz2) return kFALSE;
1189       if (imult1 != imult2) return kFALSE;
1190       if (iangle1 != iangle2) return kFALSE;
1191       return kTRUE;
1192    }
1193 }
1194
1195 //---------------------------------------------------------------------
1196 Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchPbPb2011(){
1197   //This part rejects randomly events such that the centrality gets flat for LHC11h Pb-Pb data
1198   //for 0-5% and 10-20% centrality bin
1199
1200   if (fCentralityType!="V0M") {
1201     AliWarning("Wrong value (not centrality from V0).");
1202     return -999.0;
1203   }
1204
1205   AliCentrality *centrality = fInputEvent->GetCentrality();
1206   if (!centrality) {
1207     AliWarning("Cannot get centrality from AOD event.");
1208     return -999.0;
1209   }
1210
1211   Double_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
1212   Double_t rnd_hc = -1., testf = 0.0, ff = 0, N1 = -1., N2 = -1.;
1213
1214   if(fUseCentralityPatchPbPb2011==510){
1215     N1 = 1.9404e+06;
1216     N2 = 1.56435e+06; //N2 is the reference
1217     ff = 5.04167e+06 - 1.49885e+06*cent + 2.35998e+05*cent*cent -1.22873e+04*cent*cent*cent;
1218   } else {
1219     if(fUseCentralityPatchPbPb2011==1020){
1220       N2 = 2.0e+05; //N2 is the reference
1221       N1 = 3.7e+05;
1222       ff = -1.73979e+06 - 3.05316e+06*cent + 1.05517e+06*cent*cent - 133205*cent*cent*cent + 8187.45*cent*cent*cent*cent - 247.875*cent*cent*cent*cent*cent + 2.9676*cent*cent*cent*cent*cent*cent;
1223     } else {
1224       AliError(Form("Patch for the requested centrality (%i) is not available", fUseCentralityPatchPbPb2011));
1225       return -999.0;
1226     }
1227   }
1228   testf = ( N2 + (N1-ff) ) / N1;
1229   rnd_hc = gRandom->Rndm();
1230
1231   //AliDebugClass(1, Form("Flat Centrality %d", fUseCentralityPatchPbPb2011));
1232
1233   if (rnd_hc < 0 || rnd_hc > 1 )
1234     {
1235       AliWarning("Wrong Random number generated");
1236       return -999.0;
1237     }
1238
1239   if (rnd_hc < testf)
1240     return cent;
1241   else
1242     return -999.0;
1243 }
1244 //---------------------------------------------------------------------
1245 Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchAOD049()
1246 {
1247    //
1248    //Apply centrality patch for AOD049 outliers
1249    //
1250    if (fInputEvent->InheritsFrom(AliESDEvent::Class())) {
1251       AliWarning("Requested patch for AOD049 for ESD. ");
1252       return -999.0;
1253    }
1254
1255    if (fCentralityType!="V0M") {
1256       AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
1257       return -999.0;
1258    }
1259
1260    AliCentrality *centrality = fInputEvent->GetCentrality();
1261    if (!centrality) {
1262       AliWarning("Cannot get centrality from AOD event.");
1263       return -999.0;
1264    }
1265
1266    Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
1267    /*
1268    Bool_t isSelRun = kFALSE;
1269    Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
1270    if(cent<0){
1271      Int_t quality = centrality->GetQuality();
1272      if(quality<=1){
1273        cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1274      } else {
1275        Int_t runnum=aodEvent->GetRunNumber();
1276        for(Int_t ir=0;ir<5;ir++){
1277    if(runnum==selRun[ir]){
1278      isSelRun=kTRUE;
1279      break;
1280    }
1281        }
1282        if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1283      }
1284    }
1285    */
1286
1287    if(cent>=0.0) {
1288       Float_t v0 = 0.0;
1289       AliAODEvent *aodEvent = (AliAODEvent *)fInputEvent;
1290       AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
1291       v0+=aodV0->GetMTotV0A();
1292       v0+=aodV0->GetMTotV0C();
1293       if ( (cent==0) && (v0<19500) ) {
1294          AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
1295          return -999.0;
1296       }
1297       Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1298       Float_t val = 1.30552 +  0.147931 * v0;
1299
1300       Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
1301                                120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
1302                                92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
1303                                68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
1304                                51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
1305                                37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
1306                                26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
1307                                19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
1308                                12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
1309                                13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
1310                               };
1311
1312       if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )  {
1313          AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
1314          return -999.0;
1315       }
1316    } else {
1317       //force it to be -999. whatever the negative value was
1318       cent = -999.;
1319    }
1320    return cent;
1321 }
1322
1323 //----------------------------------------------------------------------------------
1324 void AliRsnMiniAnalysisTask::SetEventQAHist(TString type,TH2F *histo)
1325 {
1326    if(!histo) {
1327       AliWarning(Form("event QA histogram pointer not defined for slot %s",type.Data()));
1328       return;
1329    }
1330
1331    type.ToLower();
1332
1333    if(!type.CompareTo("vz")) fHAEventVz = histo;
1334    else if(!type.CompareTo("multicent")) {
1335       TString mtype(histo->GetYaxis()->GetTitle());
1336       mtype.ToUpper();
1337       if(mtype.CompareTo("QUALITY") && mtype.CompareTo("TRACKS") && mtype.CompareTo("TRACKLETS")) {
1338          AliWarning(Form("multiplicity vs. centrality histogram y-axis %s unknown, setting to TRACKS",mtype.Data()));
1339          histo->GetYaxis()->SetTitle("TRACKS");
1340       }
1341       fHAEventMultiCent = histo;
1342    }
1343    else if(!type.CompareTo("eventplane")) fHAEventPlane = histo;
1344    else AliWarning(Form("event QA histogram slot %s undefined",type.Data()));
1345
1346    return;
1347 }
1348
1349 //----------------------------------------------------------------------------------
1350 Int_t AliRsnMiniAnalysisTask::CreateValue(AliRsnMiniValue::EType type, Bool_t useMC)
1351 {
1352 //
1353 // Create a new value in the task,
1354 // and returns its ID, which is needed for setting up histograms.
1355 // If that value was already initialized, returns its ID and does not recreate it.
1356 //
1357
1358    Int_t valID = ValueID(type, useMC);
1359    if (valID >= 0 && valID < fValues.GetEntries()) {
1360       AliInfo(Form("Value '%s' is already created in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
1361    } else {
1362       valID = fValues.GetEntries();
1363       AliInfo(Form("Creating value '%s' in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
1364       new (fValues[valID]) AliRsnMiniValue(type, useMC);
1365    }
1366
1367    return valID;
1368 }
1369
1370 //----------------------------------------------------------------------------------
1371 Int_t AliRsnMiniAnalysisTask::ValueID(AliRsnMiniValue::EType type, Bool_t useMC)
1372 {
1373 //
1374 // Searches if a value computation is initialized
1375 //
1376
1377    const char *name = AliRsnMiniValue::ValueName(type, useMC);
1378    TObject *obj = fValues.FindObject(name);
1379    if (obj)
1380       return fValues.IndexOf(obj);
1381    else
1382       return -1;
1383 }
1384
1385 //----------------------------------------------------------------------------------
1386 AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, AliRsnMiniOutput::EOutputType type, AliRsnMiniOutput::EComputation src)
1387 {
1388 //
1389 // Create a new histogram definition in the task,
1390 // which is then returned to the user for its configuration
1391 //
1392
1393    Int_t n = fHistograms.GetEntries();
1394    AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, type, src);
1395
1396    return newDef;
1397 }
1398
1399 //----------------------------------------------------------------------------------
1400 AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, const char *outType, const char *compType)
1401 {
1402 //
1403 // Create a new histogram definition in the task,
1404 // which is then returned to the user for its configuration
1405 //
1406
1407    Int_t n = fHistograms.GetEntries();
1408    AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, outType, compType);
1409
1410    return newDef;
1411 }