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