]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/RESONANCES/AliRsnMiniAnalysisTask.cxx
Bugfix in event mixing, addition of some features in output
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnMiniAnalysisTask.cxx
index 3db1437739c3d3e397ff6d0ca297d27c9e097997..033aae60bd50b579a022f735aa937e587e7ba9ab 100644 (file)
 
 #include <Riostream.h>
 
-#include <TList.h>
 #include <TH1.h>
-#include <TH2.h>
-#include <TH3.h>
-#include <THnSparse.h>
+#include <TList.h>
+#include <TTree.h>
 
 #include "AliLog.h"
+#include "AliEventplane.h"
+#include "AliMultiplicity.h"
+#include "AliTriggerAnalysis.h"
 #include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
+
 #include "AliESDtrackCuts.h"
 #include "AliESDUtils.h"
-#include "AliMultiplicity.h"
-#include "AliInputEventHandler.h"
-#include "AliTriggerAnalysis.h"
-#include "AliMCEvent.h"
+
 #include "AliAODEvent.h"
-#include "AliMCParticle.h"
 #include "AliAODMCParticle.h"
-#include "AliEventplane.h"
 
 #include "AliRsnCutSet.h"
+#include "AliRsnMiniPair.h"
 #include "AliRsnMiniEvent.h"
 #include "AliRsnMiniParticle.h"
-#include "AliRsnMiniPair.h"
 
 #include "AliRsnMiniAnalysisTask.h"
 
@@ -47,6 +45,7 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() :
    fEvNum(0),
    fUseCentrality(kFALSE),
    fCentralityType("QUALITY"),
+   fContinuousMix(kTRUE),
    fNMix(0),
    fMaxDiffMult(10),
    fMaxDiffVz(1.0),
@@ -59,9 +58,10 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() :
    fTrackCuts(0),
    fRsnEvent(),
    fEvBuffer(0x0),
-   fNMixed(0),
    fTriggerAna(0x0),
-   fESDtrackCuts(0x0)
+   fESDtrackCuts(0x0),
+   fMiniEvent(0x0),
+   fBigOutput(kFALSE)
 {
 //
 // Dummy constructor ALWAYS needed for I/O.
@@ -75,6 +75,7 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name, Bool_t useMC) :
    fEvNum(0),
    fUseCentrality(kFALSE),
    fCentralityType("QUALITY"),
+   fContinuousMix(kTRUE),
    fNMix(0),
    fMaxDiffMult(10),
    fMaxDiffVz(1.0),
@@ -87,9 +88,10 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name, Bool_t useMC) :
    fTrackCuts(0),
    fRsnEvent(),
    fEvBuffer(0x0),
-   fNMixed(0),
    fTriggerAna(0x0),
-   fESDtrackCuts(0x0)
+   fESDtrackCuts(0x0),
+   fMiniEvent(0x0),
+   fBigOutput(kFALSE)
 {
 //
 // Default constructor.
@@ -108,6 +110,7 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask& cop
    fEvNum(0),
    fUseCentrality(copy.fUseCentrality),
    fCentralityType(copy.fCentralityType),
+   fContinuousMix(copy.fContinuousMix),
    fNMix(copy.fNMix),
    fMaxDiffMult(copy.fMaxDiffMult),
    fMaxDiffVz(copy.fMaxDiffVz),
@@ -120,9 +123,10 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask& cop
    fTrackCuts(copy.fTrackCuts),
    fRsnEvent(),
    fEvBuffer(0x0),
-   fNMixed(0),
    fTriggerAna(copy.fTriggerAna),
-   fESDtrackCuts(copy.fESDtrackCuts)
+   fESDtrackCuts(copy.fESDtrackCuts),
+   fMiniEvent(0x0),
+   fBigOutput(copy.fBigOutput)
 {
 //
 // Copy constructor.
@@ -145,6 +149,7 @@ AliRsnMiniAnalysisTask& AliRsnMiniAnalysisTask::operator=(const AliRsnMiniAnalys
    fUseMC = copy.fUseMC;
    fUseCentrality = copy.fUseCentrality;
    fCentralityType = copy.fCentralityType;
+   fContinuousMix = copy.fContinuousMix;
    fNMix = copy.fNMix;
    fMaxDiffMult = copy.fMaxDiffMult;
    fMaxDiffVz = copy.fMaxDiffVz;
@@ -155,6 +160,7 @@ AliRsnMiniAnalysisTask& AliRsnMiniAnalysisTask::operator=(const AliRsnMiniAnalys
    fTrackCuts = copy.fTrackCuts;
    fTriggerAna = copy.fTriggerAna;
    fESDtrackCuts = copy.fESDtrackCuts;
+   fBigOutput = copy.fBigOutput;
    
    return (*this);
 }
@@ -168,8 +174,10 @@ AliRsnMiniAnalysisTask::~AliRsnMiniAnalysisTask()
 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 //
 
+
    if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
       delete fOutput;
+      delete fEvBuffer;
    }
 }
 
@@ -205,7 +213,10 @@ void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
 //
 
    // reset counter
-   fEvNum = 0;
+   fEvNum = -1;
+   
+   // message
+   AliInfo(Form("Selected event characterization: %s (%s)", (fUseCentrality ? "centrality" : "multiplicity"), fCentralityType.Data()));
 
    // initialize trigger analysis
    if (fTriggerAna) delete fTriggerAna;
@@ -216,6 +227,7 @@ void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
    fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
 
    // create list and set it as owner of its content (MANDATORY)
+   if (fBigOutput) OpenFile(1);
    fOutput = new TList();
    fOutput->SetOwner();
    
@@ -228,9 +240,9 @@ void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
    fOutput->Add(fHEventStat);
    
    // create temporary tree for filtered events
-   AliRsnMiniEvent *mini = 0x0;
+   if (fMiniEvent) delete fMiniEvent;
    fEvBuffer = new TTree("EventBuffer", "Temporary buffer for mini events");
-   fEvBuffer->Branch("events", "AliRsnMiniEvent", &mini);
+   fEvBuffer->Branch("events", "AliRsnMiniEvent", &fMiniEvent);
    
    // create one histogram per each stored definition (event histograms)
    Int_t i, ndef = fHistograms.GetEntries();
@@ -244,11 +256,6 @@ void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
       }
    }
    
-   // setup PID response
-   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager(); 
-       AliInputEventHandler *inputHandler = (AliInputEventHandler*)man->GetInputEventHandler(); 
-       fRsnEvent.SetPIDResponse(inputHandler->GetPIDResponse());
-   
    // post data for ALL output slots >0 here, to get at least an empty histogram
    PostData(1, fOutput);
 }
@@ -257,78 +264,48 @@ void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
 void AliRsnMiniAnalysisTask::UserExec(Option_t *)
 {
 //
-// Main computation loop.
-// To be precise, this loop only checks if events are OK and eventually
-// stores them in the temporary tree used for the analysis.
-// The computation of correlations is done only at the end,
-// in the automatic function which is called every time
+// Computation loop.
+// In this case, it checks if the event is acceptable, and eventually
+// creates the corresponding mini-event and stores it in the buffer.
+// The real histogram filling is done at the end, in "FinishTaskOutput".
 //
 
-   // event counter
+   // increment event counter
    fEvNum++;
-   
+
    // check current event
    Char_t check = CheckCurrentEvent();
    if (!check) return;
-      
-   // if the check is successful, the mini-event is created and stored
-   AliRsnMiniEvent *miniEvent = 0x0;
-   fEvBuffer->SetBranchAddress("events", &miniEvent);
    
-   // assign event-related values
-   miniEvent = new AliRsnMiniEvent;
-   miniEvent->Vz() = fInputEvent->GetPrimaryVertex()->GetZ();
-   miniEvent->Angle() = ComputeAngle();
-   miniEvent->Mult() = ComputeCentrality((check == 'E'));
-   AliDebugClass(1, Form("Event %d: vz = %f -- mult = %f -- angle = %f", fEvNum, miniEvent->Vz(), miniEvent->Mult(), miniEvent->Angle()));
+   // setup PID response
+   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager(); 
+       AliInputEventHandler *inputHandler = (AliInputEventHandler*)man->GetInputEventHandler(); 
+       fRsnEvent.SetPIDResponse(inputHandler->GetPIDResponse());
    
-   // fill all histograms for events only
-   Int_t id, ndef = fHistograms.GetEntries();
-   AliRsnMiniOutput *def = 0x0;
-   for (id = 0; id < ndef; id++) {
-      def = (AliRsnMiniOutput*)fHistograms[id];
-      if (!def) continue;
-      if (!def->IsEventOnly()) continue;
-      def->Fill(miniEvent, &fValues);
-   }
+   // fill a mini-event from current
+   // and skip this event if no tracks were accepted
+   FillMiniEvent(check);
    
-   // fill all histograms for mother only (if MC is present)
+   // fill MC based histograms on mothers,
+   // which do need the original event
    if (fUseMC) {
       if (fRsnEvent.IsESD() && fMCEvent)
-         FillTrueMotherESD(miniEvent);
+         FillTrueMotherESD(fMiniEvent);
       else if (fRsnEvent.IsAOD() && fRsnEvent.GetAODList())
-         FillTrueMotherAOD(miniEvent);
+         FillTrueMotherAOD(fMiniEvent);
    }
    
-   // loop on daughters
-   // and store only those that pass at least one cut
-   Int_t it, ic, ncuts = fTrackCuts.GetEntries(), ntot = fRsnEvent.GetAbsoluteSum();
-   AliRsnDaughter cursor;
-   AliRsnMiniParticle miniParticle;
-   for (it = 0; it < ntot; it++) {
-      fRsnEvent.SetDaughter(cursor, it);
-      miniParticle.CopyDaughter(&cursor);
-      for (ic = 0; ic < ncuts; ic++) {
-         AliRsnCutSet *cuts = (AliRsnCutSet*)fTrackCuts[ic];
-         if (cuts->IsSelected(&cursor)) miniParticle.SetCutBit(ic);
-      }
-      AliDebugClass(3, Form("-- Track %d: REC px = %f -- py = %f -- pz = %f", it, miniParticle.Px(0), miniParticle.Py(0), miniParticle.Pz(0)));
-      AliDebugClass(3, Form("-- Track %d: SIM px = %f -- py = %f -- pz = %f", it, miniParticle.Px(1), miniParticle.Py(1), miniParticle.Pz(1)));
-      AliDebugClass(3, Form("-- Track %d: Charge = %c -- PDG = %d -- mother = %d -- motherPDG = %d", it, miniParticle.Charge(), miniParticle.PDG(), miniParticle.Mother(), miniParticle.MotherPDG()));
-      AliDebugClass(2, Form("-- Track %d: Cutbit = %u", it, miniParticle.CutBits()));
-      if (miniParticle.CutBits()) {
-         miniEvent->AddParticle(miniParticle);
-      }
+   // if the event is not empty, store it
+   if (fMiniEvent->IsEmpty()) {
+      AliDebugClass(2, Form("Rejecting empty event #%d", fEvNum));
+   } else {
+      Int_t id = fEvBuffer->GetEntries();
+      AliDebugClass(2, Form("Adding event #%d with ID = %d", fEvNum, id));
+      fMiniEvent->ID() = id;
+      fEvBuffer->Fill();
    }
-   AliDebugClass(1, Form("Event %d: selected tracks = %d", fEvNum, miniEvent->Particles().GetEntriesFast()));
    
-   // store event
-   fEvBuffer->Fill();
-   
-   // process single event
-   ProcessEvents(miniEvent);
-   
-   // post outputs
+   // post data for computed stuff
    PostData(1, fOutput);
 }
 
@@ -336,47 +313,170 @@ void AliRsnMiniAnalysisTask::UserExec(Option_t *)
 void AliRsnMiniAnalysisTask::FinishTaskOutput()
 {
 //
-// At the end of execution, when the temporary tree is filled,
-// perform mixing with all found events
+// This function is called at the end of the loop on available events,
+// and then the buffer will be full with all the corresponding mini-events,
+// each one containing all tracks selected by each of the available track cuts.
+// Here a loop is done on each of these events, and both single-event and mixing are computed
 //
 
-   Int_t i1, i2, imix, nEvents = fEvBuffer->GetEntries();
-   AliRsnMiniEvent *evMix = 0x0, evMain;
-   fEvBuffer->SetBranchAddress("events", &evMix);
+   // security code: reassign the buffer to the mini-event cursor
+   fEvBuffer->SetBranchAddress("events", &fMiniEvent);
+   
+   // prepare variables
+   Int_t ievt, nEvents = (Int_t)fEvBuffer->GetEntries();
+   Int_t idef, nDefs   = fHistograms.GetEntries();
+   Int_t imix, iloop, ifill;
+   AliRsnMiniOutput *def = 0x0;
+   AliRsnMiniOutput::EComputation compType;
+   
+   // loop on events, and for each one fill all outputs
+   // using the appropriate procedure depending on its type
+   // only mother-related histograms are filled in UserExec,
+   // since they require direct access to MC event
+   for (ievt = 0; ievt < nEvents; ievt++) {
+      // get next entry
+      fEvBuffer->GetEntry(ievt);
+      // fill
+      for (idef = 0; idef < nDefs; idef++) {
+         def = (AliRsnMiniOutput*)fHistograms[idef];
+         if (!def) continue;
+         compType = def->GetComputation();
+         // execute computation in the appropriate way
+         switch (compType) {
+            case AliRsnMiniOutput::kEventOnly:
+               //AliDebugClass(1, Form("Event %d, def '%s': event-value histogram filling", ievt, def->GetName()));
+               ifill = 1;
+               def->FillEvent(fMiniEvent, &fValues);
+               break;
+            case AliRsnMiniOutput::kTruePair:
+               //AliDebugClass(1, Form("Event %d, def '%s': true-pair histogram filling", ievt, def->GetName()));
+               ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
+               break;
+            case AliRsnMiniOutput::kTrackPair:
+               //AliDebugClass(1, Form("Event %d, def '%s': pair-value histogram filling", ievt, def->GetName()));
+               ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
+               break;
+            case AliRsnMiniOutput::kTrackPairRotated1:
+               //AliDebugClass(1, Form("Event %d, def '%s': rotated (1) background histogram filling", ievt, def->GetName()));
+               ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
+               break;
+            case AliRsnMiniOutput::kTrackPairRotated2:
+               //AliDebugClass(1, Form("Event %d, def '%s': rotated (2) background histogram filling", ievt, def->GetName()));
+               ifill = def->FillPair(fMiniEvent, fMiniEvent, &fValues);
+               break;
+            default:
+               // other kinds are processed elsewhere
+               ifill = 0;
+               AliDebugClass(2, Form("Computation = %d", (Int_t)compType));
+         }
+         // message
+         AliDebugClass(1, Form("Event %6d: def = '%15s' -- fills = %5d", ievt, def->GetName(), ifill));
+      }
+   }
+   
+   // if no mixing is required, stop here and post the output
+   if (fNMix < 1) {
+      AliDebugClass(2, "Stopping here, since no mixing is required");
+      PostData(1, fOutput);
+      return;
+   }
    
    // initialize mixing counter
-   fNMixed.Set(nEvents);
-   TString msg();
-   
-   // loop on events
-   for (i1 = 0; i1 < nEvents; i1++) {
-      if (fNMixed[i1] >= fNMix) continue;
-      fEvBuffer->GetEntry(i1);
-      evMain = (*evMix);
-      for (i2 = 1; i2 < nEvents; i2++) {
-         imix = i1 + i2;
+   Int_t    nmatched[nEvents];
+   TString *smatched = new TString[nEvents];
+   for (ievt = 0; ievt < nEvents; ievt++) {
+      smatched[ievt] = "|";
+      nmatched[ievt] = 0;
+   }
+   
+   // search for good matchings
+   for (ievt = 0; ievt < nEvents; ievt++) {
+      if (nmatched[ievt] >= fNMix) continue;
+      fEvBuffer->GetEntry(ievt);
+      AliRsnMiniEvent evMain(*fMiniEvent);
+      for (iloop = 1; iloop < nEvents; iloop++) {
+         imix = ievt + iloop;
          if (imix >= nEvents) imix -= nEvents;
-         if (imix == i1) continue;
-         if (fNMixed[i1] >= fNMix) break;
-         if (fNMixed[imix] >= fNMix) continue;
+         if (imix == ievt) continue;
+         // text next entry
          fEvBuffer->GetEntry(imix);
-         // exit if events are not matched
-         if (TMath::Abs(evMain.Vz() - evMix->Vz()) > fMaxDiffVz) continue;
-         if (TMath::Abs(evMain.Mult() - evMix->Mult()) > fMaxDiffMult) continue;
-         if (TMath::Abs(evMain.Angle() - evMix->Angle()) > fMaxDiffAngle) continue;
-         // found a match: increment counter for both events
-         fNMixed[i1]++;
-         fNMixed[imix]++;
-         //cout << "Mixed " << i1 << " with " << imix << endl;
-         ProcessEvents(&evMain, evMix);
-         // if mixed enough times, stop
+         // skip if events are not matched
+         if (!EventsMatch(&evMain, fMiniEvent)) continue;
+         // check that the array of good matches for mixed does not already contain main event
+         if (smatched[imix].Contains(Form("|%d|", ievt))) continue;
+         // check that the found good events has not enough matches already
+         if (nmatched[imix] >= fNMix) continue;
+         // add new mixing candidate
+         smatched[ievt].Append(Form("%d|", imix));
+         nmatched[ievt]++;
+         nmatched[imix]++;
+         if (nmatched[ievt] >= fNMix) break;
       }
+      AliDebugClass(1, Form("Matches for event %5d = %d [%s] (missing are declared above)", evMain.ID(), nmatched[ievt], smatched[ievt].Data()));
    }
    
-   // message
-   for (Int_t ii = 0; ii < nEvents; ii++) cout << Form("Event #%6d mixed %3d times", ii, fNMixed[ii]) << endl;
+   // perform mixing
+   TObjArray *list = 0x0;
+   TObjString *os = 0x0;
+   for (ievt = 0; ievt < nEvents; ievt++) {
+      ifill = 0;
+      fEvBuffer->GetEntry(ievt);
+      AliRsnMiniEvent evMain(*fMiniEvent);
+      list = smatched[ievt].Tokenize("|");
+      TObjArrayIter next(list);
+      while ( (os = (TObjString*)next()) ) {
+         imix = os->GetString().Atoi();
+         fEvBuffer->GetEntry(imix);
+         for (idef = 0; idef < nDefs; idef++) {
+            def = (AliRsnMiniOutput*)fHistograms[idef];
+            if (!def) continue;
+            if (!def->IsTrackPairMix()) continue;
+            ifill += def->FillPair(&evMain, fMiniEvent, &fValues, kTRUE);
+            if (!def->IsSymmetric()) {
+               AliDebugClass(2, "Reflecting non symmetric pair");
+               ifill += def->FillPair(fMiniEvent, &evMain, &fValues, kFALSE);
+            }
+         }
+      }
+      delete list;
+   }
    
-   // post outputs
+   delete [] smatched;
+      
+   /*
+   OLD   
+   ifill = 0;
+   for (iloop = 1; iloop < nEvents; iloop++) {
+      imix = ievt + iloop;
+      // restart from beginning if reached last event
+      if (imix >= nEvents) imix -= nEvents;
+      // avoid to mix an event with itself
+      if (imix == ievt) continue;
+      // skip all events already mixed enough times
+      if (fNMixed[ievt] >= fNMix) break;
+      if (fNMixed[imix] >= fNMix) continue;
+      fEvBuffer->GetEntry(imix);
+      // skip if events are not matched
+      if (TMath::Abs(evMain.Vz()    - fMiniEvent->Vz()   ) > fMaxDiffVz   ) continue;
+      if (TMath::Abs(evMain.Mult()  - fMiniEvent->Mult() ) > fMaxDiffMult ) continue;
+      if (TMath::Abs(evMain.Angle() - fMiniEvent->Angle()) > fMaxDiffAngle) continue;
+      // found a match: increment counter for both events
+      AliDebugClass(1, Form("Event %d, def '%s': event mixing (%d with %d)", ievt, def->GetName(), ievt, imix));
+      fNMixed[ievt]++;
+      fNMixed[imix]++;
+      // process mixing
+      ifill += def->FillPair(&evMain, fMiniEvent, &fValues);
+      // stop if mixed enough times
+      if (fNMixed[ievt] >= fNMix) break;
+   }
+   break;      
+   // print number of mixings done with each event
+   for (ievt = 0; ievt < nEvents; ievt++) {
+      AliDebugClass(2, Form("Event %6d: mixed %2d times", ievt, fNMixed[ievt]));
+   }
+   */
+   
+   // post computed data
    PostData(1, fOutput);
 }
 
@@ -431,7 +531,7 @@ Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent()
       // --> if this is failed, the event is rejected
       isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
       if (!isSelected) {
-         AliDebugClass(1, "Event does not pass physics selections");
+         AliDebugClass(2, "Event does not pass physics selections");
          fRsnEvent.SetRef(0x0);
          fRsnEvent.SetRefMC(0x0);
          return 0;
@@ -521,7 +621,7 @@ Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent()
    }
    
    // if the above exit point is not taken, the event is accepted
-   AliDebugClass(1, Form("Stats for event %d: %s", fEvNum, msg.Data()));
+   AliDebugClass(2, Form("Stats: %s", msg.Data()));
    if (isSelected) {
       fHEventStat->Fill(3.1);
       return output;
@@ -530,6 +630,52 @@ Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent()
    }
 }
 
+//__________________________________________________________________________________________________
+void AliRsnMiniAnalysisTask::FillMiniEvent(Char_t evType)
+{
+//
+// Refresh cursor mini-event data member to fill with current event.
+// Returns the total number of tracks selected.
+//
+
+   // assign event-related values
+   if (fMiniEvent) delete fMiniEvent;
+   fMiniEvent = new AliRsnMiniEvent;
+   fMiniEvent->Vz()    = fInputEvent->GetPrimaryVertex()->GetZ();
+   fMiniEvent->Angle() = ComputeAngle();
+   fMiniEvent->Mult()  = ComputeCentrality((evType == 'E'));
+   AliDebugClass(2, Form("Event %d: type = %c -- vz = %f -- mult = %f -- angle = %f", fEvNum, evType, fMiniEvent->Vz(), fMiniEvent->Mult(), fMiniEvent->Angle()));
+   
+   // loop on daughters and assign track-related values
+   Int_t ic, ncuts = fTrackCuts.GetEntries();
+   Int_t ip, npart = fRsnEvent.GetAbsoluteSum();
+   Int_t npos = 0, nneg = 0, nneu = 0;
+   AliRsnDaughter cursor;
+   AliRsnMiniParticle miniParticle;
+   for (ip = 0; ip < npart; ip++) {
+      // point cursor to next particle
+      fRsnEvent.SetDaughter(cursor, ip);
+      // copy momentum and MC info if present
+      miniParticle.CopyDaughter(&cursor);
+      miniParticle.Index() = ip;
+      // switch on the bits corresponding to passed cuts
+      for (ic = 0; ic < ncuts; ic++) {
+         AliRsnCutSet *cuts = (AliRsnCutSet*)fTrackCuts[ic];
+         if (cuts->IsSelected(&cursor)) miniParticle.SetCutBit(ic);
+      }
+      // if a track passes at least one track cut, it is added to the pool
+      if (miniParticle.CutBits()) {
+         fMiniEvent->AddParticle(miniParticle);
+         if (miniParticle.Charge() == '+') npos++;
+         else if (miniParticle.Charge() == '-') nneg++;
+         else nneu++;
+      }
+   }
+   
+   // get number of accepted tracks
+   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));
+}
+
 //__________________________________________________________________________________________________
 Double_t AliRsnMiniAnalysisTask::ComputeAngle()
 {
@@ -537,7 +683,15 @@ Double_t AliRsnMiniAnalysisTask::ComputeAngle()
 // Get the plane angle
 //
 
-   AliEventplane *plane = fInputEvent->GetEventplane();
+   AliEventplane *plane = 0x0;
+   
+   if (fInputEvent->InheritsFrom(AliESDEvent::Class()))
+      plane = fInputEvent->GetEventplane();
+   else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) {
+      AliAODEvent *aodEvent = (AliAODEvent*)fInputEvent;
+      plane = aodEvent->GetHeader()->GetEventplaneP();
+   }
+
    if (plane) 
       return plane->GetEventplane("Q");
    else {
@@ -618,7 +772,7 @@ void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
       if (!def->IsMother()) continue;
       for (ip = 0; ip < npart; ip++) {
          AliMCParticle *part = (AliMCParticle*)fMCEvent->GetTrack(ip);
-         if (TMath::Abs(part->Particle()->GetPdgCode()) != def->GetMotherPDG()) continue;
+         if (part->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
          // check that daughters match expected species
          if (part->Particle()->GetNDaughters() < 2) continue;
          label1 = part->Particle()->GetDaughter(0);
@@ -640,7 +794,7 @@ void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
          miniPair.FillRef(def->GetMotherMass());
          // do computations
-         def->Fill(&miniPair, miniEvent, &fValues);
+         def->FillMother(&miniPair, miniEvent, &fValues);
       }
    }
 }
@@ -667,7 +821,7 @@ void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
       if (!def->IsMother()) continue;
       for (ip = 0; ip < npart; ip++) {
          AliAODMCParticle *part = (AliAODMCParticle*)list->At(ip);
-         if (TMath::Abs(part->GetPdgCode()) != def->GetMotherPDG()) continue;
+         if (part->GetPdgCode() != def->GetMotherPDG()) continue;
          // check that daughters match expected species
          if (part->GetNDaughters() < 2) continue;
          label1 = part->GetDaughter(0);
@@ -689,56 +843,54 @@ void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
          miniPair.FillRef(def->GetMotherMass());
          // do computations
-         def->Fill(&miniPair, miniEvent, &fValues);
+         def->FillMother(&miniPair, miniEvent, &fValues);
       }
    }
 }
 
-//________________________________________________________________________
-void AliRsnMiniAnalysisTask::ProcessEvents(AliRsnMiniEvent *evMain, AliRsnMiniEvent *evMix)
+//__________________________________________________________________________________________________
+Bool_t AliRsnMiniAnalysisTask::EventsMatch(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2)
 {
 //
-// This function processes all single tracks.
-// Loops on the daughters defined in current event
-// and on all of them which pass the single-track cuts defined, computes what needed
+// Check if two events are compatible.
+// If the mixing is continuous, this is true if differences in vz, mult and angle are smaller than
+// the specified values. 
+// If the mixing is binned, this is true if the events are in the same bin.
 //
 
-   // check if it is mixed
-   Bool_t isMix;
-   if (evMix) {
-      isMix = kTRUE;
+   if (!event1 || !event2) return kFALSE;
+   Int_t ivz1, ivz2, imult1, imult2, iangle1, iangle2;
+   Double_t dv, dm, da;
+   
+   if (fContinuousMix) {
+      dv = TMath::Abs(event1->Vz()    - event2->Vz()   );
+      dm = TMath::Abs(event1->Mult()  - event2->Mult() );
+      da = TMath::Abs(event1->Angle() - event2->Angle());
+      if (dv > fMaxDiffVz) {
+         //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Vz = %f", event1->ID(), event2->ID(), dv));
+         return kFALSE;
+      }
+      if (dm > fMaxDiffMult ) {
+         //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Mult = %f", event1->ID(), event2->ID(), dm));
+         return kFALSE;
+      }
+      if (da > fMaxDiffAngle) {
+         //AliDebugClass(2, Form("Events #%4d and #%4d don't match due to a too large diff in Angle = %f", event1->ID(), event2->ID(), da));
+         return kFALSE;
+      }
+      return kTRUE;
    } else {
-      evMix = evMain;
-      isMix = kFALSE;
+      ivz1 = (Int_t)(event1->Vz() / fMaxDiffVz);
+      ivz2 = (Int_t)(event2->Vz() / fMaxDiffVz);
+      imult1 = (Int_t)(event1->Mult() / fMaxDiffMult);
+      imult2 = (Int_t)(event2->Mult() / fMaxDiffMult);
+      iangle1 = (Int_t)(event1->Angle() / fMaxDiffAngle);
+      iangle2 = (Int_t)(event2->Angle() / fMaxDiffAngle);
+      if (ivz1 != ivz2) return kFALSE;
+      if (imult1 != imult2) return kFALSE;
+      if (iangle1 != iangle2) return kFALSE;
+      return kTRUE;
    }
-   AliDebugClass(1, Form("Event %d: processing %s", fEvNum, (isMix ? "mixing" : "single event")));
-   
-   // 2 nested loops on tracks
-   // inner starts from next track or from first, depending if is mixing or not
-   Int_t i1, i2, idef, start;
-   Int_t n1 = evMain->Particles().GetEntries();
-   Int_t n2 = evMix ->Particles().GetEntries();
-   Int_t ndef = fHistograms.GetEntries();
-   AliRsnMiniParticle *daughter1 = 0x0, *daughter2 = 0x0;
-   AliRsnMiniOutput *def = 0x0;
-   for (i1 = 0; i1 < n1; i1++) {
-      daughter1 = (AliRsnMiniParticle*)evMain->Particles()[i1];
-      if (!daughter1) continue;
-      if (isMix) start = 0; else start = i1 + 1;
-      for (i2 = start; i2 < n2; i2++) {
-         daughter2 = (AliRsnMiniParticle*)evMix->Particles()[i2];
-         if (!daughter2) continue;
-         // loop on definitions
-         for (idef = 0; idef < ndef; idef++) {
-            def = (AliRsnMiniOutput*)fHistograms[idef];
-            if (!def) continue;
-            if (def->IsEventOnly()) continue;
-            if (def->IsMother()) continue;
-            if (isMix && !def->IsTrackPairMix()) continue;
-            if (!isMix && def->IsTrackPairMix()) continue;
-            // check if definitions here match the current pair
-            def->Fill(daughter1, daughter2, evMain, &fValues);
-         } // end loop on definitions
-      } // end loop on daughter #2
-   } // end loop on daughter #1
 }
+
+