PWG2/RESONANCES/AliRsnMiniAnalysisTask.cxx

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