]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/RESONANCES/AliRsnMiniAnalysisTask.cxx
Fix for coverity 24399, 24400, 24401
[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / AliRsnMiniAnalysisTask.cxx
index be893b6313220d59d7ff22bce13acede2b22f1d3..9ae2d8c2313dc1768e4eff18d161cdf34b52a37e 100644 (file)
@@ -16,6 +16,7 @@
 #include <TList.h>
 #include <TTree.h>
 #include <TStopwatch.h>
+#include "TRandom.h"
 
 #include "AliLog.h"
 #include "AliEventplane.h"
@@ -45,9 +46,11 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() :
    AliAnalysisTaskSE(),
    fUseMC(kFALSE),
    fEvNum(0),
+   fTriggerMask(0),
    fUseCentrality(kFALSE),
    fCentralityType("QUALITY"),
    fUseAOD049CentralityPatch(kFALSE),
+   fUseCentralityPatchPbPb2011(0),
    fContinuousMix(kTRUE),
    fNMix(0),
    fMaxDiffMult(10),
@@ -57,6 +60,10 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() :
    fHistograms("AliRsnMiniOutput", 0),
    fValues("AliRsnMiniValue", 0),
    fHEventStat(0x0),
+   fHAEventsVsMulti(0x0),
+   fHAEventVz(0x0),
+   fHAEventMultiCent(0x0),
+   fHAEventPlane(0x0),
    fEventCuts(0x0),
    fTrackCuts(0),
    fRsnEvent(),
@@ -65,7 +72,14 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() :
    fESDtrackCuts(0x0),
    fMiniEvent(0x0),
    fBigOutput(kFALSE),
-   fMixPrintRefresh(-1)
+   fMixPrintRefresh(-1),
+   fMaxNDaughters(-1),
+   fCheckP(kFALSE),
+   fCheckFeedDown(kFALSE),   
+   fOriginDselection(kFALSE),
+   fKeepDfromB(kFALSE),
+   fKeepDfromBOnly(kFALSE),
+   fRejectIfNoQuark(kFALSE)
 {
 //
 // Dummy constructor ALWAYS needed for I/O.
@@ -77,9 +91,11 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name, Bool_t useMC) :
    AliAnalysisTaskSE(name),
    fUseMC(useMC),
    fEvNum(0),
+   fTriggerMask(AliVEvent::kMB),
    fUseCentrality(kFALSE),
    fCentralityType("QUALITY"),
    fUseAOD049CentralityPatch(kFALSE),
+  fUseCentralityPatchPbPb2011(0),
    fContinuousMix(kTRUE),
    fNMix(0),
    fMaxDiffMult(10),
@@ -89,6 +105,10 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name, Bool_t useMC) :
    fHistograms("AliRsnMiniOutput", 0),
    fValues("AliRsnMiniValue", 0),
    fHEventStat(0x0),
+   fHAEventsVsMulti(0x0),
+   fHAEventVz(0x0),
+   fHAEventMultiCent(0x0),
+   fHAEventPlane(0x0),
    fEventCuts(0x0),
    fTrackCuts(0),
    fRsnEvent(),
@@ -97,7 +117,14 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name, Bool_t useMC) :
    fESDtrackCuts(0x0),
    fMiniEvent(0x0),
    fBigOutput(kFALSE),
-   fMixPrintRefresh(-1)
+   fMixPrintRefresh(-1),
+   fMaxNDaughters(-1),
+   fCheckP(kFALSE),
+   fCheckFeedDown(kFALSE),   
+   fOriginDselection(kFALSE),
+   fKeepDfromB(kFALSE),
+   fKeepDfromBOnly(kFALSE),
+   fRejectIfNoQuark(kFALSE)
 {
 //
 // Default constructor.
@@ -114,9 +141,11 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask &cop
    AliAnalysisTaskSE(copy),
    fUseMC(copy.fUseMC),
    fEvNum(0),
+   fTriggerMask(copy.fTriggerMask),
    fUseCentrality(copy.fUseCentrality),
    fCentralityType(copy.fCentralityType),
    fUseAOD049CentralityPatch(copy.fUseAOD049CentralityPatch),
+   fUseCentralityPatchPbPb2011(copy.fUseCentralityPatchPbPb2011),
    fContinuousMix(copy.fContinuousMix),
    fNMix(copy.fNMix),
    fMaxDiffMult(copy.fMaxDiffMult),
@@ -126,6 +155,10 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask &cop
    fHistograms(copy.fHistograms),
    fValues(copy.fValues),
    fHEventStat(0x0),
+   fHAEventsVsMulti(0x0),
+   fHAEventVz(0x0),
+   fHAEventMultiCent(0x0),
+   fHAEventPlane(0x0),
    fEventCuts(copy.fEventCuts),
    fTrackCuts(copy.fTrackCuts),
    fRsnEvent(),
@@ -134,7 +167,14 @@ AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask &cop
    fESDtrackCuts(copy.fESDtrackCuts),
    fMiniEvent(0x0),
    fBigOutput(copy.fBigOutput),
-   fMixPrintRefresh(copy.fMixPrintRefresh)
+   fMixPrintRefresh(copy.fMixPrintRefresh),
+   fMaxNDaughters(copy.fMaxNDaughters),
+   fCheckP(copy.fCheckP),
+   fCheckFeedDown(copy.fCheckFeedDown),   
+   fOriginDselection(copy.fOriginDselection),
+   fKeepDfromB(copy.fOriginDselection),
+   fKeepDfromBOnly(copy.fKeepDfromBOnly),
+   fRejectIfNoQuark(copy.fRejectIfNoQuark)
 {
 //
 // Copy constructor.
@@ -156,9 +196,12 @@ AliRsnMiniAnalysisTask &AliRsnMiniAnalysisTask::operator=(const AliRsnMiniAnalys
    if (this == &copy)
       return *this;
    fUseMC = copy.fUseMC;
+   fEvNum = copy.fEvNum;
+   fTriggerMask = copy.fTriggerMask;
    fUseCentrality = copy.fUseCentrality;
    fCentralityType = copy.fCentralityType;
    fUseAOD049CentralityPatch = copy.fUseAOD049CentralityPatch;
+   fUseCentralityPatchPbPb2011 = copy.fUseCentralityPatchPbPb2011;
    fContinuousMix = copy.fContinuousMix;
    fNMix = copy.fNMix;
    fMaxDiffMult = copy.fMaxDiffMult;
@@ -166,12 +209,24 @@ AliRsnMiniAnalysisTask &AliRsnMiniAnalysisTask::operator=(const AliRsnMiniAnalys
    fMaxDiffAngle = copy.fMaxDiffAngle;
    fHistograms = copy.fHistograms;
    fValues = copy.fValues;
+   fHEventStat = copy.fHEventStat;
+   fHAEventsVsMulti = copy.fHAEventsVsMulti;
+   fHAEventVz = copy.fHAEventVz;
+   fHAEventMultiCent = copy.fHAEventMultiCent;
+   fHAEventPlane = copy.fHAEventPlane;
    fEventCuts = copy.fEventCuts;
    fTrackCuts = copy.fTrackCuts;
    fTriggerAna = copy.fTriggerAna;
    fESDtrackCuts = copy.fESDtrackCuts;
    fBigOutput = copy.fBigOutput;
    fMixPrintRefresh = copy.fMixPrintRefresh;
+   fMaxNDaughters = copy.fMaxNDaughters;
+   fCheckP = copy.fCheckP;
+   fCheckFeedDown = copy.fCheckFeedDown;
+   fOriginDselection = copy.fOriginDselection;
+   fKeepDfromB = copy.fOriginDselection;
+   fKeepDfromBOnly = copy.fKeepDfromBOnly;
+   fRejectIfNoQuark = copy.fRejectIfNoQuark;
    return (*this);
 }
 
@@ -249,6 +304,16 @@ void AliRsnMiniAnalysisTask::UserCreateOutputObjects()
    fHEventStat->GetXaxis()->SetBinLabel(4, "Accepted");
    fOutput->Add(fHEventStat);
 
+   if (fUseCentrality)
+      fHAEventsVsMulti = new TH1F("hAEventsVsMulti", "Accepted events vs Centrality", 100, 0, 100.0);
+   else
+      fHAEventsVsMulti = new TH1F("hAEventsVsMulti", "Accepted events vs Multiplicity",1000, 0, 1000.0);
+   fOutput->Add(fHAEventsVsMulti);
+
+   if(fHAEventVz) fOutput->Add(fHAEventVz);
+   if(fHAEventMultiCent) fOutput->Add(fHAEventMultiCent);
+   if(fHAEventPlane) fOutput->Add(fHAEventPlane);
+
    TIter next(&fTrackCuts);
    AliRsnCutSet *cs;
    while ((cs = (AliRsnCutSet *) next())) {
@@ -575,7 +640,8 @@ Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent()
       output = 'E';
       // ESD specific check: Physics Selection
       // --> if this is failed, the event is rejected
-      isSelected = (((AliInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
+      isSelected = (((AliInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
+
       if (!isSelected) {
          AliDebugClass(2, "Event does not pass physics selections");
          fRsnEvent.SetRef(0x0);
@@ -670,6 +736,11 @@ Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent()
    AliDebugClass(2, Form("Stats: %s", msg.Data()));
    if (isSelected) {
       fHEventStat->Fill(3.1);
+      Double_t multi = ComputeCentrality((output == 'E'));
+      fHAEventsVsMulti->Fill(multi);
+      if(fHAEventVz) fHAEventVz->Fill(multi,fInputEvent->GetPrimaryVertex()->GetZ());
+      if(fHAEventMultiCent) fHAEventMultiCent->Fill(multi,ComputeMultiplicity(output == 'E',fHAEventMultiCent->GetYaxis()->GetTitle()));
+      if(fHAEventPlane) fHAEventPlane->Fill(multi,ComputeAngle());
       return output;
    } else {
       return 0;
@@ -756,17 +827,19 @@ Double_t AliRsnMiniAnalysisTask::ComputeCentrality(Bool_t isESD)
 //
 
    if (fUseCentrality) {
-
-      if ((!fUseMC) && (!isESD) && (fUseAOD049CentralityPatch)) {
-         return ApplyCentralityPatchAOD049();
-      } else {
-         AliCentrality *centrality = fInputEvent->GetCentrality();
+     if ((!fUseMC) && (fUseCentralityPatchPbPb2011)) {
+       return ApplyCentralityPatchPbPb2011();//
+    }
+     if ((!fUseMC) && (!isESD) && (fUseAOD049CentralityPatch)) {
+       return ApplyCentralityPatchAOD049();
+     } else {
+       AliCentrality *centrality = fInputEvent->GetCentrality();
          if (!centrality) {
-            AliError("Cannot compute centrality!");
-            return -1.0;
+          AliError("Cannot compute centrality!");
+          return -1.0;
          }
          return centrality->GetCentralityPercentile(fCentralityType.Data());
-      }
+     }
    } else {
       if (!fCentralityType.CompareTo("TRACKS"))
          return fInputEvent->GetNumberOfTracks();
@@ -802,6 +875,49 @@ Double_t AliRsnMiniAnalysisTask::ComputeCentrality(Bool_t isESD)
    }
 }
 
+//__________________________________________________________________________________________________
+Double_t AliRsnMiniAnalysisTask::ComputeMultiplicity(Bool_t isESD,TString type)
+{
+//
+// Computes event multiplicity according to the string defining
+// what criterion must be used for specific computation.
+//
+
+   type.ToUpper();
+
+   if (!type.CompareTo("TRACKS"))
+      return fInputEvent->GetNumberOfTracks();
+   else if (!type.CompareTo("QUALITY"))
+      if (isESD)
+         return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent *)fInputEvent, kTRUE);
+      else {
+         Double_t count = 0.;
+         Int_t iTrack, ntracksLoop = fInputEvent->GetNumberOfTracks();
+         for (iTrack = 0; iTrack < ntracksLoop; iTrack++) {
+            AliVTrack   *track = (AliVTrack *)fInputEvent->GetTrack(iTrack);
+            AliAODTrack *aodt  = dynamic_cast<AliAODTrack *>(track);
+            if (!aodt) continue;
+            if (!aodt->TestFilterBit(5)) continue;
+            count++;
+         }
+         return count;
+      }
+   else if (!type.CompareTo("TRACKLETS")) {
+      if (isESD) {
+         const AliMultiplicity *mult = ((AliESDEvent *)fInputEvent)->GetMultiplicity();
+         Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
+         for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
+         return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ());
+      } else {
+         AliWarning("Cannot compute multiplicity with SPD tracklets from AOD");
+         return 1E20;
+      }
+   } else {
+      AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", type.Data()));
+      return -1.0;
+   }
+}
+
 //__________________________________________________________________________________________________
 void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
 {
@@ -823,9 +939,11 @@ void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
       if (!def->IsMother()) continue;
       for (ip = 0; ip < npart; ip++) {
          AliMCParticle *part = (AliMCParticle *)fMCEvent->GetTrack(ip);
+         //get mother pdg code
          if (part->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
          // check that daughters match expected species
-         if (part->Particle()->GetNDaughters() < 2) continue;
+         if (part->Particle()->GetNDaughters() < 2) continue; 
+        if (fMaxNDaughters > 0 && part->Particle()->GetNDaughters() > fMaxNDaughters) continue;
          label1 = part->Particle()->GetDaughter(0);
          label2 = part->Particle()->GetDaughter(1);
          daughter1 = (AliMCParticle *)fMCEvent->GetTrack(label1);
@@ -841,6 +959,58 @@ void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent)
             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
          }
          if (!okMatch) continue;
+        if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&   
+                               (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
+                               (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
+        if(fCheckFeedDown){
+               Int_t pdgGranma = 0;
+               Int_t mother = 0;
+               mother = part->GetMother();
+               Int_t istep = 0;
+               Int_t abspdgGranma =0;
+               Bool_t isFromB=kFALSE;
+               Bool_t isQuarkFound=kFALSE;
+               while (mother >=0 ){
+                       istep++;
+                       AliDebug(2,Form("mother at step %d = %d", istep, mother));
+                       AliMCParticle* mcGranma = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(mother));
+                       if (mcGranma){
+                               pdgGranma = mcGranma->PdgCode();
+                               AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
+                               abspdgGranma = TMath::Abs(pdgGranma);
+                               if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
+                                 isFromB=kTRUE;
+                               }
+                               if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
+                               mother = mcGranma->GetMother();
+                       }else{
+                               AliError("Failed casting the mother particle!");
+                               break;
+                       }
+               }
+               if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
+               if(isFromB){
+                 if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
+               }
+               else{ 
+                 if (fKeepDfromBOnly) pdgGranma = -999;
+                 
+               if (pdgGranma == -99999){
+                       AliDebug(2,"This particle does not have a quark in his genealogy\n");
+                       continue;
+               }
+               if (pdgGranma == -9999){
+                       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");   
+                       continue;
+               }       
+               
+               if (pdgGranma == -999){
+                       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");  
+                       continue;
+               }       
+
+             }
+        }
          // assign momenta to computation object
          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
          miniPair.FillRef(def->GetMotherMass());
@@ -875,6 +1045,7 @@ void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
          if (part->GetPdgCode() != def->GetMotherPDG()) continue;
          // check that daughters match expected species
          if (part->GetNDaughters() < 2) continue;
+        if (fMaxNDaughters > 0 && part->GetNDaughters() > fMaxNDaughters) continue;
          label1 = part->GetDaughter(0);
          label2 = part->GetDaughter(1);
          daughter1 = (AliAODMCParticle *)list->At(label1);
@@ -890,7 +1061,58 @@ void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
             p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0));
          }
          if (!okMatch) continue;
-         // assign momenta to computation object
+        if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&   
+                               (TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
+                               (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
+        if(fCheckFeedDown){
+               Int_t pdgGranma = 0;
+               Int_t mother = 0;
+               mother = part->GetMother();
+               Int_t istep = 0;
+               Int_t abspdgGranma =0;
+               Bool_t isFromB=kFALSE;
+               Bool_t isQuarkFound=kFALSE;
+               while (mother >=0 ){
+                       istep++;
+                       AliDebug(2,Form("mother at step %d = %d", istep, mother));
+                       AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(list->At(mother));
+                       if (mcGranma){
+                               pdgGranma = mcGranma->GetPdgCode();
+                               AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
+                               abspdgGranma = TMath::Abs(pdgGranma);
+                               if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
+                                 isFromB=kTRUE;
+                               }
+                               if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
+                               mother = mcGranma->GetMother();
+                       }else{
+                               AliError("Failed casting the mother particle!");
+                               break;
+                       }
+               }
+               if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
+               if(isFromB){
+                 if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
+               }
+               else{ 
+                 if (fKeepDfromBOnly) pdgGranma = -999;
+                 }
+                 
+               if (pdgGranma == -99999){
+                       AliDebug(2,"This particle does not have a quark in his genealogy\n");
+                       continue;
+               }
+               if (pdgGranma == -9999){
+                       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");   
+                       continue;
+               }       
+               
+               if (pdgGranma == -999){
+                       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");  
+                       continue;
+               }       
+        } 
+        // assign momenta to computation object
          miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
          miniPair.FillRef(def->GetMotherMass());
          // do computations
@@ -898,7 +1120,33 @@ void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent)
       }
    }
 }
-
+//___________________________________________________________
+void AliRsnMiniAnalysisTask::SetDselection(UShort_t originDselection)
+{
+       // setting the way the D0 will be selected
+       // 0 --> only from c quarks
+       // 1 --> only from b quarks
+       // 2 --> from both c quarks and b quarks
+               
+       fOriginDselection = originDselection;
+       
+       if (fOriginDselection == 0) {
+               fKeepDfromB = kFALSE;
+               fKeepDfromBOnly = kFALSE;
+       }
+       
+       if (fOriginDselection == 1) {
+               fKeepDfromB = kTRUE;
+               fKeepDfromBOnly = kTRUE;
+       }
+       
+       if (fOriginDselection == 2) {
+               fKeepDfromB = kTRUE;
+               fKeepDfromBOnly = kFALSE;
+       }
+       
+       return; 
+}
 //__________________________________________________________________________________________________
 Bool_t AliRsnMiniAnalysisTask::EventsMatch(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2)
 {
@@ -944,7 +1192,55 @@ Bool_t AliRsnMiniAnalysisTask::EventsMatch(AliRsnMiniEvent *event1, AliRsnMiniEv
    }
 }
 
+//---------------------------------------------------------------------
+Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchPbPb2011(){
+  //This part rejects randomly events such that the centrality gets flat for LHC11h Pb-Pb data
+  //for 0-5% and 10-20% centrality bin
+  
+  if (fCentralityType!="V0M") {
+    AliWarning("Wrong value (not centrality from V0).");
+    return -999.0;
+  }
+  
+  AliCentrality *centrality = fInputEvent->GetCentrality();
+  if (!centrality) {
+    AliWarning("Cannot get centrality from AOD event.");
+    return -999.0;
+  }
+  
+  Double_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));               
+  Double_t rnd_hc = -1., testf = 0.0, ff = 0, N1 = -1., N2 = -1.;
+
+  if(fUseCentralityPatchPbPb2011==510){
+    N1 = 1.9404e+06;
+    N2 = 1.56435e+06; //N2 is the reference 
+    ff = 5.04167e+06 - 1.49885e+06*cent + 2.35998e+05*cent*cent -1.22873e+04*cent*cent*cent;
+  } else {
+    if(fUseCentralityPatchPbPb2011==1020){
+      N2 = 2.0e+05; //N2 is the reference
+      N1 = 3.7e+05;
+      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;
+    } else {
+      AliError(Form("Patch for the requested centrality (%i) is not available", fUseCentralityPatchPbPb2011));
+      return -999.0;
+    }
+  }
+  testf = ( N2 + (N1-ff) ) / N1;
+  rnd_hc = gRandom->Rndm();
 
+  //AliDebugClass(1, Form("Flat Centrality %d", fUseCentralityPatchPbPb2011));
+
+  if (rnd_hc < 0 || rnd_hc > 1 ) 
+    {
+      AliWarning("Wrong Random number generated");
+      return -999.0;
+    }
+  
+  if (rnd_hc < testf)
+    return cent;
+  else
+    return -999.0;
+}
 //---------------------------------------------------------------------
 Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchAOD049()
 {
@@ -1023,3 +1319,93 @@ Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchAOD049()
    }
    return cent;
 }
+
+//----------------------------------------------------------------------------------
+void AliRsnMiniAnalysisTask::SetEventQAHist(TString type,TH2F *histo)
+{
+   if(!histo) {
+      AliWarning(Form("event QA histogram pointer not defined for slot %s",type.Data()));
+      return;
+   }
+
+   type.ToLower();
+
+   if(!type.CompareTo("vz")) fHAEventVz = histo;
+   else if(!type.CompareTo("multicent")) {
+      TString mtype(histo->GetYaxis()->GetTitle());
+      mtype.ToUpper();
+      if(mtype.CompareTo("QUALITY") && mtype.CompareTo("TRACKS") && mtype.CompareTo("TRACKLETS")) {
+         AliWarning(Form("multiplicity vs. centrality histogram y-axis %s unknown, setting to TRACKS",mtype.Data()));
+         histo->GetYaxis()->SetTitle("TRACKS");
+      }
+      fHAEventMultiCent = histo;
+   }
+   else if(!type.CompareTo("eventplane")) fHAEventPlane = histo;
+   else AliWarning(Form("event QA histogram slot %s undefined",type.Data()));
+
+   return;
+}
+
+//----------------------------------------------------------------------------------
+Int_t AliRsnMiniAnalysisTask::CreateValue(AliRsnMiniValue::EType type, Bool_t useMC)
+{
+//
+// Create a new value in the task,
+// and returns its ID, which is needed for setting up histograms.
+// If that value was already initialized, returns its ID and does not recreate it.
+//
+
+   Int_t valID = ValueID(type, useMC);
+   if (valID >= 0 && valID < fValues.GetEntries()) {
+      AliInfo(Form("Value '%s' is already created in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
+   } else {
+      valID = fValues.GetEntries();
+      AliInfo(Form("Creating value '%s' in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID));
+      new (fValues[valID]) AliRsnMiniValue(type, useMC);
+   }
+
+   return valID;
+}
+
+//----------------------------------------------------------------------------------
+Int_t AliRsnMiniAnalysisTask::ValueID(AliRsnMiniValue::EType type, Bool_t useMC)
+{
+//
+// Searches if a value computation is initialized
+//
+
+   const char *name = AliRsnMiniValue::ValueName(type, useMC);
+   TObject *obj = fValues.FindObject(name);
+   if (obj)
+      return fValues.IndexOf(obj);
+   else
+      return -1;
+}
+
+//----------------------------------------------------------------------------------
+AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, AliRsnMiniOutput::EOutputType type, AliRsnMiniOutput::EComputation src)
+{
+//
+// Create a new histogram definition in the task,
+// which is then returned to the user for its configuration
+//
+
+   Int_t n = fHistograms.GetEntries();
+   AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, type, src);
+
+   return newDef;
+}
+
+//----------------------------------------------------------------------------------
+AliRsnMiniOutput *AliRsnMiniAnalysisTask::CreateOutput(const char *name, const char *outType, const char *compType)
+{
+//
+// Create a new histogram definition in the task,
+// which is then returned to the user for its configuration
+//
+
+   Int_t n = fHistograms.GetEntries();
+   AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, outType, compType);
+
+   return newDef;
+}