]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardFlowTaskQC.cxx
index 011bd57071fd695d8a4ac14d9acc4088968237dc..f3315bc8b8af11040212c9adcad9b2d95294a8de 100644 (file)
@@ -35,8 +35,9 @@
 #include "AliCentrality.h"
 #include "AliESDEvent.h"
 #include "AliVTrack.h"
-#include "AliESDtrackCuts.h"
+#include "AliESDtrack.h"
 #include "AliAODTrack.h"
+#include "AliAnalysisFilter.h"
 
 ClassImp(AliForwardFlowTaskQC)
 #if 0
@@ -50,13 +51,13 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC()
     fFMDCut(-1),         // FMD sigma cut
     fSPDCut(-1),         // SPD sigma cut
     fFlowFlags(0),       // Flow flags
-    fEtaGap(-1),         // Eta gap value
+    fEtaGap(0.),         // Eta gap value
     fBinsForward(),      // List with forward flow hists
     fBinsCentral(),      // List with central flow hists
     fSumList(0),        // Event sum list
     fOutputList(0),     // Result output list
     fAOD(0),            // AOD input event
-    fESDTrackCuts(0),    // ESD track cuts
+    fTrackCuts(0),    // ESD track cuts
     fMaxMoment(0),       // Max flow moment
     fVtx(1111),                 // Z vertex coordinate
     fCent(-1),          // Centrality
@@ -78,14 +79,14 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name)
     fCentAxis(),        // Axis to control centrality/multiplicity binning
     fFMDCut(-1),        // FMD sigma cut
     fSPDCut(-1),        // SPD sigma cut
-    fFlowFlags(0x0),       // Flow flags
-    fEtaGap(2.),        // Eta gap value
+    fFlowFlags(0x0),    // Flow flags
+    fEtaGap(0.),        // Eta gap value
     fBinsForward(),     // List with forward flow hists
     fBinsCentral(),     // List with central flow hists
     fSumList(0),        // Event sum list           
     fOutputList(0),     // Result output list       
     fAOD(0),           // AOD input event          
-    fESDTrackCuts(0),   // ESD track cuts
+    fTrackCuts(0),      // ESD track cuts
     fMaxMoment(4),      // Max flow moment
     fVtx(1111),         // Z vertex coordinate      
     fCent(-1),          // Centrality               
@@ -119,7 +120,7 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o)
     fSumList(o.fSumList),              // Event sum list           
     fOutputList(o.fOutputList),        // Result output list       
     fAOD(o.fAOD),                     // AOD input event          
-    fESDTrackCuts(o.fESDTrackCuts),    // ESD track cuts
+    fTrackCuts(o.fTrackCuts),          // ESD track cuts
     fMaxMoment(o.fMaxMoment),          // Flow moments
     fVtx(o.fVtx),                      // Z vertex coordinate      
     fCent(o.fCent),                   // Centrality
@@ -154,7 +155,7 @@ AliForwardFlowTaskQC::operator=(const AliForwardFlowTaskQC& o)
   fSumList        = o.fSumList;
   fOutputList     = o.fOutputList;
   fAOD            = o.fAOD;
-  fESDTrackCuts   = o.fESDTrackCuts;
+  fTrackCuts      = o.fTrackCuts;
   fMaxMoment      = o.fMaxMoment;
   fVtx            = o.fVtx;
   fCent           = o.fCent;
@@ -190,12 +191,7 @@ void AliForwardFlowTaskQC::UserCreateOutputObjects()
   //
   InitVertexBins();
   InitHists();
-  if (fFlowFlags & kTPC) {
-    fESDTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
-    fESDTrackCuts->SetPtRange(0.2, 5.0);
-    fESDTrackCuts->SetEtaRange(-0.8, 0.8);
-    fESDTrackCuts->SetMinNClustersTPC(70);
-  }
+  if ((fFlowFlags & kTracks) && !fTrackCuts) AliFatal("No track cuts set!");
   PrintFlowSetup();
 
   PostData(1, fSumList);
@@ -211,11 +207,13 @@ void AliForwardFlowTaskQC::InitVertexBins()
     Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
     if ((fFlowFlags & kFMD)) {
       fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "FMD", fFlowFlags, fFMDCut, fEtaGap));
-      if (!(fFlowFlags & k3Cor)) fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-FMD", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
+      if (!(fFlowFlags & k3Cor)) 
+       fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-FMD", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
     }
     else if ((fFlowFlags & kVZERO)) {
       fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "VZERO", fFlowFlags, 0, fEtaGap));
-      if ((fFlowFlags & kEtaGap)) fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-VZERO", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
+      if ((fFlowFlags & kEtaGap) && !(fFlowFlags & kTracks)) 
+       fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-VZERO", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
     }
   }
 }
@@ -396,20 +394,17 @@ void AliForwardFlowTaskQC::FillVtxBinList(const TList& list, TH2D& h, Int_t vtx,
   Int_t nVtxBins = fVtxAxis->GetNbins();
   
   while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
+    i++;
     // If no tracks do things normally
-    if (!(fFlowFlags & kTPC) && !bin->FillHists(h, fCent, kFillBoth|flags|kReset)) return;
+    if (!(fFlowFlags & kTracks) || (flags & kMC)) {
+      if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
+    }
     // if tracks things are more complicated
-    else if ((fFlowFlags & kTPC)) {
-      TObjArray* trList = GetTracks();
-      if (!trList) return;
-      Bool_t useEvent = bin->FillTracks(trList, kFillRef|kReset|flags);
-      // If esd input trList is a new object owned by this task and should be cleaned up
-      if (AliForwardUtil::CheckForAOD() == 2) delete trList;
-      if (!useEvent) return;
-      if (!bin->FillHists(h, fCent, kFillDiff|kReset|flags)) return;
+    else if ((fFlowFlags & kTracks)) {
+      if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
+      if (!bin->FillHists(h, fCent, kFillDiff|kReset|flags)) continue;
     }
     bin->CumulantsAccumulate(fCent);
-    i++;
   }
 
   return;
@@ -434,10 +429,15 @@ void AliForwardFlowTaskQC::FillVtxBinListEtaGap(const TList& list, TH2D& href,
   Int_t nVtxBins = fVtxAxis->GetNbins();
 
   while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
-    if (!bin->FillHists(href, fCent, kFillRef|flags|kReset)) return;
-    bin->FillHists(hdiff, fCent, kFillDiff|kReset);
-    bin->CumulantsAccumulate(fCent);
     i++;
+    if (!(fFlowFlags & kTracks) || (flags & kMC)) {
+      if(!bin->FillHists(href, fCent, kFillRef|flags|kReset)) continue;
+    }
+    else if ((fFlowFlags & kTracks)) {
+      if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
+    }
+    if (!bin->FillHists(hdiff, fCent, kFillDiff|kReset|flags)) continue;
+    bin->CumulantsAccumulate(fCent);
   }
 
   return;
@@ -464,9 +464,9 @@ void AliForwardFlowTaskQC::FillVtxBinList3Cor(const TList& list, TH2D& hcent,
   TH2D& h = CombineHists(hcent, hfwd);
 
   while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
-    if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) return;
-    bin->CumulantsAccumulate3Cor(fCent);
     i++;
+    if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
+    bin->CumulantsAccumulate3Cor(fCent);
   }
 
   return;
@@ -552,7 +552,7 @@ TH2D& AliForwardFlowTaskQC::CombineHists(TH2D& hcent, TH2D& hfwd)
   return fHistdNdedp3Cor;
 }
 //_____________________________________________________________________
-TObjArray* AliForwardFlowTaskQC::GetTracks() const
+Bool_t AliForwardFlowTaskQC::FillTracks(VertexBin* bin, UShort_t mode) const
 {
   // 
   //  Get TPC tracks to use for reference flow.
@@ -560,25 +560,14 @@ TObjArray* AliForwardFlowTaskQC::GetTracks() const
   //  Return: TObjArray with tracks
   //
   TObjArray* trList = 0;
-  // Get input type
-  UShort_t input = AliForwardUtil::CheckForAOD();
-  switch (input) {
-    // If AOD input, simply get the track array from the event
-    case 1: trList = static_cast<TObjArray*>(fAOD->GetTracks());
-           break;
-    case 2: {
-    // If ESD input get event, apply track cuts
-             AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
-             if (!esd) return 0;
-             // Warning! trList is now a new array, we need to delete it after use
-             // this is not a very good implementation!
-             trList = fESDTrackCuts->GetAcceptedTracks(esd, kTRUE);
-             break;
-           }
-    default: AliFatal("Neither ESD or AOD input. This should never happen");
-           break;
-  }
-  return trList;
+  AliESDEvent* esdEv = 0;
+  if (AliForwardUtil::CheckForAOD() == 1) // AOD tracks
+    trList = static_cast<TObjArray*>(fAOD->GetTracks());
+  else 
+    esdEv = dynamic_cast<AliESDEvent*>(InputEvent());
+  
+  Bool_t useEvent = bin->FillTracks(trList, esdEv, fTrackCuts, mode);
+  return useEvent;
 }
 //_____________________________________________________________________
 void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
@@ -601,7 +590,7 @@ void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
   fOutputList->SetName("Results");
   fOutputList->SetOwner();
 
-  if ((fFlowFlags & kEtaGap)) {
+  if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks)) {
     TParameter<Double_t>* etaGap = new TParameter<Double_t>("EtaGap", fEtaGap);
     fOutputList->Add(etaGap);
   }
@@ -865,6 +854,11 @@ Bool_t AliForwardFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm)
 // _____________________________________________________________________
 AliVVZERO* AliForwardFlowTaskQC::GetVZERO() const
 {
+  //
+  //  Get VZERO object from ESD or AOD
+  //
+  //  Return: VZERO data object
+  //
   AliVVZERO* vzero = 0;
   // Get input type
   UShort_t input = AliForwardUtil::CheckForAOD();
@@ -905,10 +899,11 @@ void AliForwardFlowTaskQC::FillVZEROHist(AliVVZERO* vzero)
                      0.521922, 0.785915, 0.703658, 0.832479, 0.77461, 0.73129, 0.778697, 0.710265, 
                      0.89686, 0.967688, 0.974225, 0.873445, 0.811096, 0.828493, 0.889609, 0.586056, 
                      1.15877, 0.954656, 0.914557, 0.979028, 1.04907, 0.748518, 0.928043, 0.98175 };
+
   for (Int_t i = 0; i < 64; i++) {
     if (i % 8 == 0) {
       ring++;
-      bin = (ring < 5 ? ring+1 : 15-ring);
+      bin = (ring < 5 ? 11-ring : ring-3);
       eta = fHistdNdedpV0.GetXaxis()->GetBinCenter(bin);
       fHistdNdedpV0.SetBinContent(bin, 0, 1);
     }
@@ -1058,7 +1053,7 @@ void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist, TAxis* centAx
 
   // Get bin numbers and binning defined
   Int_t nHBins = GetBinNumberSin();
-  Int_t nEtaBins = 48
+  Int_t nEtaBins = 24
   if ((fFlags & k3Cor)) {
     if ((fFlags & kFMD)) nEtaBins = 24;
     else if ((fFlags & kVZERO)) nEtaBins = 19;
@@ -1074,7 +1069,7 @@ void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist, TAxis* centAx
 
   Int_t nRefBins = nEtaBins; // needs to be something as default
   if ((fFlags & kStdQC)) {
-    if ((fFlags & kSymEta) && !((fFlags & kTPC) && (fFlags & kSPD))) nRefBins = 1;
+    if ((fFlags & kSymEta) && !((fFlags & kTracks) && (fFlags & kSPD))) nRefBins = 1;
     else nRefBins = 2;
   } else if ((fFlags & kEtaGap )) {
     nRefBins = 2;
@@ -1207,7 +1202,6 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D& dNdetadphi, Double_t cen
     if ((mode & kFillRef))  fCumuRef->Reset();
     if ((mode & kFillDiff)) fCumuDiff->Reset();
   }
-
   // Then we loop over the input and calculate sum cos(k*n*phi)
   // and fill it in the reference and differential histograms
   Int_t nBadBins = 0;
@@ -1228,7 +1222,7 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D& dNdetadphi, Double_t cen
             TMath::Abs(eta) < fEtaGap) break;
        // Backward and forward eta gap break for reference flow
        if ((fFlags & kEtaGap) && (mode & kFillRef) && TMath::Abs(eta) > TMath::Abs(limit)) break;
-       if ((fFlags & kStdQC) && (fFlags & kMC)) {
+       if ((fFlags & kStdQC) && (fFlags & kMC) && !(fFlags & kTracks)) {
          if (!(fFlags & kSPD) && TMath::Abs(eta) < 1.75) break; 
          if ((fFlags & kSPD) && TMath::Abs(eta) > 2.00) break;
        }
@@ -1244,9 +1238,8 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D& dNdetadphi, Double_t cen
       nInAvg++;
       if (weight == 0) continue;
       if (weight > max) max = weight;
-      
       // Fill into Cos() and Sin() hists
-      if ((mode & kFillRef)) {
+      if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
        fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
         fdNdedpRefAcc->Fill(eta, phi, weight);
       }
@@ -1260,7 +1253,7 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D& dNdetadphi, Double_t cen
        Double_t cosnPhi = weight*TMath::Cos(n*phi);
        Double_t sinnPhi = weight*TMath::Sin(n*phi);
         // fill ref
-       if ((mode & kFillRef)) {
+       if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
          fCumuRef->Fill(eta, cosBin, cosnPhi);
          fCumuRef->Fill(eta, sinBin, sinnPhi);
        }
@@ -1289,7 +1282,8 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D& dNdetadphi, Double_t cen
   return useEvent;
 }
 //_____________________________________________________________________
-Bool_t AliForwardFlowTaskQC::VertexBin::FillTracks(TObjArray* trList, UShort_t mode) 
+Bool_t AliForwardFlowTaskQC::VertexBin::FillTracks(TObjArray* trList, AliESDEvent* esd,
+                                                   AliAnalysisFilter* trFilter, UShort_t mode) 
 {
   // 
   //  Fill reference and differential eta-histograms
@@ -1299,6 +1293,10 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillTracks(TObjArray* trList, UShort_t m
   //   mode: filling mode: kFillRef/kFillDiff/kFillBoth
   //
   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
+  if (!trList && !esd) {
+    AliError("FillTracks: No AOD track list or ESD event - something might be wrong!");
+    return kFALSE;
+  }
 
   // Fist we reset histograms
   if ((mode & kReset)) {
@@ -1306,42 +1304,54 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillTracks(TObjArray* trList, UShort_t m
     if ((mode & kFillDiff)) fCumuDiff->Reset();
   }
 
-  UShort_t input = AliForwardUtil::CheckForAOD();
   // Then we loop over the input and calculate sum cos(k*n*phi)
   // and fill it in the reference and differential histograms
-  Int_t nTr = trList->GetEntries();
+  Int_t nTr = 0;
+  if (trList) nTr = trList->GetEntries();
+  if (esd) nTr = esd->GetNumberOfTracks();
   if (nTr == 0) return kFALSE;
   AliVTrack* tr = 0;
-  AliAODTrack* aodTr = 0;
-  // Cuts for AOD tracks (have already been applied to ESD tracks)
-  const Double_t pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70;
+  // Cuts for AOD tracks (have already been applied to ESD tracks) - except dEdx
+//  const tpcdEdx = 10;
   for (Int_t i = 0; i < nTr; i++) { // track loop
-    tr = (AliVTrack*)trList->At(i);
+    tr = (trList ? (AliVTrack*)trList->At(i) : (AliVTrack*)esd->GetTrack(i));
     if (!tr) continue;
-    if (input == 1) { // If AOD input
-      // A dynamic cast would be more safe here, but this is faster...
-      aodTr = (AliAODTrack*)tr;
+    if (esd) {
+      AliESDtrack* esdTr = (AliESDtrack*)tr;
+      if (!trFilter->IsSelected(esdTr)) continue;
+    }
+    else if (trList) { // If AOD input
+      Double_t pTMin = 0, pTMax = 0, etaMin = 0, etaMax = 0, minNCl = 0;
+      UInt_t bit = 0;
+      if ((fFlags & kTPC) == kTPC)    pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 128;
+      if ((fFlags & kHybrid) == kHybrid) pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 272;
+
+      AliAODTrack* aodTr = (AliAODTrack*)tr;
       if (aodTr->GetID() > -1) continue;
-      if (!aodTr->TestFilterBit(128) || !aodTr->Pt() > pTMax || aodTr->Pt() < pTMin || 
+      if (!aodTr->TestFilterBit(bit) || aodTr->Pt() > pTMax || aodTr->Pt() < pTMin || 
        aodTr->Eta() > etaMax || aodTr->Eta() < etaMin || aodTr->GetTPCNcls() < minNCl) continue;
     }
+
+//    if (tr->GetTPCsignal() < tpcdEdx) continue;
     // Track accepted
     Double_t eta = tr->Eta();
-//    if ((fFlags & kSPD) && TMath::Abs(eta) < 0.4) continue;
+    if (((fFlags & kSPD) || (fFlags & kEtaGap)) && TMath::Abs(eta) < fEtaGap) continue;
     Double_t phi = tr->Phi();
+    Double_t weight = 1.;
+
     if ((mode & kFillRef)) {
-      fCumuRef->Fill(eta, 0.);// mult goes in underflowbin - no visual, but not needed?
-      fdNdedpRefAcc->Fill(eta, phi);
+      fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
+      fdNdedpRefAcc->Fill(eta, phi, weight);
     }
     if ((mode & kFillDiff)) {
-      fCumuDiff->Fill(eta, 0);
-      fdNdedpDiffAcc->Fill(eta, phi);
+      fCumuDiff->Fill(eta, 0., weight);
+      fdNdedpDiffAcc->Fill(eta, phi, weight);
     }
     for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
       Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
       Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
-      Double_t cosnPhi = TMath::Cos(n*phi);
-      Double_t sinnPhi = TMath::Sin(n*phi);
+      Double_t cosnPhi = weight*TMath::Cos(n*phi);
+      Double_t sinnPhi = weight*TMath::Sin(n*phi);
       // fill ref
       if ((mode & kFillRef)) {
        fCumuRef->Fill(eta, cosBin, cosnPhi);
@@ -1370,14 +1380,16 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
   // Fill out NUA hists
   for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
     Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
-    if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
-    if ((fFlags & kTPC) && (fFlags && kSPD)) eta = -eta;
+    if (fCumuRef->GetBinContent(etaBin, 0) <= 3) continue;
+    if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) eta = -eta;
     for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
       fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
     }
   }
   for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
     Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
+    Double_t refetaBin = fCumuRef->GetXaxis()->FindBin(eta);
+    if (fCumuRef->GetBinContent(refetaBin, 0) <= 3) continue;
     if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
     for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
       fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
@@ -1403,12 +1415,11 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
     for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
       Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
       Double_t refEta = eta;
-      if ((fFlags & kTPC) && (fFlags && kSPD)) refEta = -eta;
+      if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) refEta = -eta;
       Int_t refEtaBinA = fCumuRef->GetXaxis()->FindBin(refEta);
       if ((fFlags & kEtaGap)) refEta = -eta;
       Int_t refEtaBinB = fCumuRef->GetXaxis()->FindBin(refEta);
       if (refEtaBinA != prevRefEtaBin) {
-       prevRefEtaBin = refEtaBinA;
        // Reference flow
        multA  = fCumuRef->GetBinContent(refEtaBinA, 0);
        dQnReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n));
@@ -1420,7 +1431,7 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
        dQnReB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberCos(n));
        dQnImB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberSin(n));
 
-       if (multA <= 3 || multB <= 3) return
+       if (multA <= 3 || multB <= 3) continue
        // The reference flow is calculated 
        // 2-particle
        if ((fFlags & kStdQC)) {
@@ -1462,6 +1473,7 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
          cumuRef->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
          cumuRef->Fill(eta, cent, k3pWeight, multA*(multA-1.)*(multA-2.));
        } // End of QC{4}
+       prevRefEtaBin = refEtaBinA;
       } // End of reference flow
       // For each etaBin bin the necessary values for differential flow is calculated
       Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
@@ -1478,7 +1490,7 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
 
       // Differential flow calculations for each eta bin is done:
       // 2-particle differential flow
-      if ((fFlags & kStdQC) && !(fFlags & kTPC)) {
+      if ((fFlags & kStdQC) && (!(fFlags & kTracks) || ((fFlags & kTracks) && (fFlags & kMC) && !(fFlags & kSPD) && TMath::Abs(eta) < 0.75))) {
        mq = mp;
        qnRe = pnRe;
        qnIm = pnIm;
@@ -1976,7 +1988,7 @@ void AliForwardFlowTaskQC::VertexBin::CalculateReferenceFlow(CumuHistos& cumu2h,
          // The next line covers both cases.
          qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
          // Extra NUA term from 2n cosines and sines
-         Double_t den = 1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB);
+         Double_t den = 1+(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB);
          if (den != 0) qc2 /= den;
          else qc2 = 0;
        }
@@ -2124,7 +2136,7 @@ void AliForwardFlowTaskQC::VertexBin::CalculateDifferentialFlow(CumuHistos& cumu
          // Old nua
          qc2Prime -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB;
          // Extra NUA term from 2n cosines and sines
-         qc2Prime /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
+         qc2Prime /= (1.+(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
        }
        if (!TMath::IsNaN(qc2Prime)) {
          quality->Fill((n-2)*qualityFactor+3, Int_t(cent));
@@ -2287,7 +2299,10 @@ void AliForwardFlowTaskQC::VertexBin::Calculate3CorFlow(CumuHistos& cumu2h, TH2I
            sin2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
            multB += fCumuNUARef->GetBinContent(b, cBin, 0);
          }
-         if (multA == 0 || multB == 0) AliFatal("Empty NUA values for 3Cor!");
+         if (multA == 0 || multB == 0) {
+           AliWarning(Form("Empty NUA values for 3Cor! (%s)", cumuRef->GetName()));
+           continue;
+         }
          cosP1nPhiA /= multA;
          sinP1nPhiA /= multA;
          cos2nPhiA /= multA;
@@ -2308,7 +2323,7 @@ void AliForwardFlowTaskQC::VertexBin::Calculate3CorFlow(CumuHistos& cumu2h, TH2I
          // Old nua
          qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
          // Extra NUA term from 2n cosines and sines
-         qc2 /= (1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
+         qc2 /= (1+(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
        }
        if (qc2 <= 0) { 
          if (fDebug > 0) 
@@ -2354,8 +2369,8 @@ void AliForwardFlowTaskQC::VertexBin::Calculate3CorFlow(CumuHistos& cumu2h, TH2I
          qc2PrimeA -= cosP1nPsi*cosP1nPhiA + sinP1nPsi*sinP1nPhiA;
          qc2PrimeB -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB; // Is this OK?
          // Extra NUA term from 2n cosines and sines
-         qc2PrimeA /= (1.-(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
-         qc2PrimeB /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
+         if (cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA != -1.) qc2PrimeA /= (1.+(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
+         if (cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB != -1.) qc2PrimeB /= (1.+(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
        }
        if (!TMath::IsNaN(qc2PrimeA) && !TMath::IsNaN(qc2PrimeB) && qc2 != 0) {
        if (qc2PrimeA*qc2PrimeB >= 0) {
@@ -2800,16 +2815,17 @@ void AliForwardFlowTaskQC::PrintFlowSetup() const
   for (Int_t n  = 2; n <= fMaxMoment; n++) printf("v%d ", n);
   printf("\n");
   TString type = "Standard QC{2} and QC{4} calculations.";
-  if ((fFlowFlags & kTPC)) type.ReplaceAll(".", " with TPC tracks for reference.");
   if ((fFlowFlags & kEtaGap)) type = "QC{2} with a rapidity gap.";
   if ((fFlowFlags & k3Cor))   type = "QC{2} with 3 correlators.";
+  if ((fFlowFlags & kTPC) == kTPC)    type.ReplaceAll(".", " with TPC tracks for reference.");
+  if ((fFlowFlags & kHybrid) == kHybrid) type.ReplaceAll(".", " with hybrid tracks for reference.");
   Printf("QC calculation type               :\t%s", type.Data());
   Printf("Symmetrize ref. flow wrt. eta = 0 :\t%s", ((fFlowFlags & kSymEta) ? "true" : "false"));
   Printf("Apply NUA correction terms        :\t%s", ((fFlowFlags & kNUAcorr) ? "true" : "false"));
   Printf("Satellite vertex flag             :\t%s", ((fFlowFlags & kSatVtx) ? "true" : "false"));
   Printf("FMD sigma cut:                    :\t%f", fFMDCut);
   Printf("SPD sigma cut:                    :\t%f", fSPDCut);
-  if ((fFlowFlags & kEtaGap)) 
+  if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks)
     Printf("Eta gap:                          :\t%f", fEtaGap);
   Printf("=======================================================");
 }
@@ -2834,7 +2850,9 @@ const Char_t* AliForwardFlowTaskQC::GetQCType(UShort_t flags, Bool_t prependUS)
   else if ((flags & k3Cor))   type = "3Cor";
   else type = "UNKNOWN";
   if (prependUS) type.Prepend("_");
-  if ((flags & kTPC)) type.Append("TPCTr");
+  if ((flags & kTPC) == kTPC)    type.Append("TPCTr");
+  if ((flags & kHybrid) == kHybrid) type.Append("HybTr");
+  
   return type.Data();
 }
 //_____________________________________________________________________