]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Lots of small improvements and optimizations to flow task.
authorhansena <hansena@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Mar 2012 11:32:17 +0000 (11:32 +0000)
committerhansena <hansena@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Mar 2012 11:32:17 +0000 (11:32 +0000)
PWGLF/FORWARD/analysis2/AddTaskForwardFlow.C
PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx
PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.h
PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx
PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.h
PWGLF/FORWARD/analysis2/MakeFlow.C
PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C

index 9767dbb82a96fe7c26fc7b95b20b1a1713a1b394..15408d9c641cf85b09982ec5ff7d589f891d1106 100644 (file)
@@ -82,10 +82,18 @@ void AddTaskForwardFlow(TString type = "",
     task = new AliForwardFlowTaskQC("QCumulants");
   mgr->AddTask(task); 
 
+  // Set which harmonics to do
   task->SetDoHarmonics(v1, v2, v3, v4, v5, v6);
+  // Set non-default axis for vertices
+  TAxis* a = new TAxis(20, -10, 10);
+  task->SetVertexAxis(a);
+  // Set debug flag
+  task->SetDebugLevel(0);
+  // Set up adding flow to MC input
   if (mc) {
     AliForwardMCFlowTaskQC* mcTask = 
       static_cast<AliForwardMCFlowTaskQC*>task;
+    mcTask->SetUseImpactParameter(true);
     mcTask->AddFlow(addFlow);
     mcTask->AddFlowType(addFType);
     mcTask->AddFlowOrder(addFOrder);
index 1feaac37fe18b9cf96051d468ee3da64516600bf..785b2e38b6bdf67cf661a1b72f925bbd022958f5 100644 (file)
@@ -33,12 +33,13 @@ ClassImp(AliForwardFlowTaskQC)
 
 AliForwardFlowTaskQC::AliForwardFlowTaskQC()
   : AliAnalysisTaskSE(),
+    fVtxAxis(),         // Axis to contorl vertex binning
     fBinsFMD(),         // List with FMD flow histos
     fBinsSPD(),         // List with SPD flow histos
     fSumList(0),       // Event sum list
     fOutputList(0),    // Result output list
     fAOD(0),           // AOD input event
-    fZvertex(1111),    // Z vertex coordinate
+    fVtx(1111),        // Z vertex coordinate
     fCent(-1),         // Centrality
     fHistCent(),        // Histo for centrality
     fHistVertexSel(),   // Histo for selected vertices
@@ -52,12 +53,13 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC()
 //_____________________________________________________________________
 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) 
   : AliAnalysisTaskSE(name),
+    fVtxAxis(),         // Axis to contorl vertex binning
     fBinsFMD(),         // List with FMD flow histos
     fBinsSPD(),         // List with SPD flow histos
     fSumList(0),        // Event sum list           
     fOutputList(0),     // Result output list       
     fAOD(0),           // AOD input event          
-    fZvertex(1111),     // Z vertex coordinate      
+    fVtx(1111),     // Z vertex coordinate      
     fCent(-1),          // Centrality               
     fHistCent(),        // Histo for centrality
     fHistVertexSel(),   // Histo for selected vertices
@@ -71,6 +73,7 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name)
   //
   for (Int_t n = 0; n <= 6; n++) fv[n] = kTRUE;
 
+  fVtxAxis       = new TAxis(20, -10, 10);
   fHistCent      = new TH1D("hCent", "Centralities", 100, 0, 100);
   fHistVertexSel = new TH1D("hVertexSel", "Selected vertices", 40, -20, 20);
   fHistVertexAll = new TH1D("hVertexAll", "All vertices", 40, -20, 20);
@@ -81,12 +84,13 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name)
 //_____________________________________________________________________
 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o)
   : AliAnalysisTaskSE(o),
+    fVtxAxis(o.fVtxAxis),              // Axis to contorl vertex binning
     fBinsFMD(),                        // List with FMD flow histos
     fBinsSPD(),                        // List with SPD flow histos
     fSumList(o.fSumList),              // Event sum list           
     fOutputList(o.fOutputList),        // Result output list       
     fAOD(o.fAOD),                     // AOD input event          
-    fZvertex(o.fZvertex),              // Z vertex coordinate      
+    fVtx(o.fVtx),              // Z vertex coordinate      
     fCent(o.fCent),                   // Centrality               
     fHistCent(o.fHistCent),            // Histo for centrality
     fHistVertexSel(o.fHistVertexSel),  // Histo for selected vertices
@@ -108,10 +112,11 @@ AliForwardFlowTaskQC::operator=(const AliForwardFlowTaskQC& o)
   // Assignment operator 
   //
   if (&o == this) return *this;
+  fVtxAxis       = o.fVtxAxis;
   fSumList       = o.fSumList;
   fOutputList    = o.fOutputList;
   fAOD           = o.fAOD;
-  fZvertex       = o.fZvertex;
+  fVtx       = o.fVtx;
   fCent          = o.fCent;
   fHistCent      = o.fHistCent;
   fHistVertexSel = o.fHistVertexSel;
@@ -140,11 +145,11 @@ void AliForwardFlowTaskQC::InitVertexBins()
   // 
   // Init vertexbin objects for FMD and SPD, and add them to the lists
   //
-  for(Int_t n = 1; n <= 6; n++) {
+  for(UShort_t n = 1; n <= 6; n++) {
     if (!fv[n]) continue;
-    for (Int_t v = -10; v < 10; v++) {
-      fBinsFMD.Add(new VertexBin(v, v+1, n, "FMD"));
-      fBinsSPD.Add(new VertexBin(v, v+1, n, "SPD"));
+    for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
+      fBinsFMD.Add(new VertexBin(fVtxAxis->GetBinLowEdge(v), fVtxAxis->GetBinUpEdge(v), n, "FMD"));
+      fBinsSPD.Add(new VertexBin(fVtxAxis->GetBinLowEdge(v), fVtxAxis->GetBinUpEdge(v), n, "SPD", kFALSE));
     }
   }
 
@@ -162,6 +167,8 @@ void AliForwardFlowTaskQC::InitHists()
 
   TList* dList = new TList();
   dList->SetName("Diagnostics");
+  fVtxAxis->SetName("VtxAxis");
+  dList->Add(fVtxAxis);
   dList->Add(fHistCent);
   dList->Add(fHistVertexSel);
   dList->Add(fHistVertexAll);
@@ -186,11 +193,12 @@ void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
   // Parameters:
   //  option: Not used
   //
-
+  
   Analyze();
   
   PostData(1, fSumList);
 
+  return;
 }
 //_____________________________________________________________________
 Bool_t AliForwardFlowTaskQC::Analyze()
@@ -202,40 +210,52 @@ Bool_t AliForwardFlowTaskQC::Analyze()
 
   // Reset data members
   fCent = -1;
-  fZvertex = 1111;
+  fVtx = 1111;
 
   // Get input event
   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
   if (!fAOD) return kFALSE;
 
-  AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
-  if (!aodfmult) return kFALSE;
-  if (!AODCheck(aodfmult)) return kFALSE;
-  TH2D fmddNdetadphi = aodfmult->GetHistogram();
-
-  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
-  if (!aodcmult) return kFALSE;
-  TH2D spddNdetadphi = aodcmult->GetHistogram();
+  const AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
+  const AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
+  if (!aodfmult || !aodcmult) return kFALSE;
+  
+  // Check event for triggers, get centrality, vtx etc.
+  if (!CheckEvent(aodfmult)) return kFALSE;
 
-  // TODO: remove me!
-//  fCent = 0.5;
+  // If everything is OK: get histos
+  const TH2D& fmddNdetadphi = aodfmult->GetHistogram();
+  const TH2D& spddNdetadphi = aodcmult->GetHistogram();
 
   // Run analysis on FMD and SPD
-  TIter nextFMD(&fBinsFMD);
+  Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
+  if (!FillVtxBinList(fBinsFMD, fmddNdetadphi, vtx)) return kFALSE;
+  if (!FillVtxBinList(fBinsSPD, spddNdetadphi, vtx)) return kFALSE;
+
+  return kTRUE;
+}
+//_____________________________________________________________________
+Bool_t AliForwardFlowTaskQC::FillVtxBinList(const TList& list, const TH2D& h, Int_t vtx) const
+{
+  //
+  // Loops over list of VtxBins, fills hists of bins for current vertex
+  // and runs analysis on those bins
+  //
+  // Parameters:
+  //  list: list of VtxBins
+  //  h:    dN/detadphi histogram
+  //  vBin: current vertex bin
+  //
+  // return true on success
+  //
   VertexBin* bin = 0;
-  while ((bin = static_cast<VertexBin*>(nextFMD()))) {
-    if (bin->CheckVertex(fZvertex)) {
-      if (!bin->FillHists(&fmddNdetadphi)) return kFALSE; 
-      bin->CumulantsAccumulate(fCent);
-    }
-  }
+  Int_t i = 0;
+  Int_t nVtxBins = fVtxAxis->GetNbins();
 
-  TIter nextSPD(&fBinsSPD);
-  while ((bin = static_cast<VertexBin*>(nextSPD()))) {
-    if (bin->CheckVertex(fZvertex)) {
-      if (!bin->FillHists(&spddNdetadphi)) return kFALSE;
-      bin->CumulantsAccumulate(fCent);
-    }
+  while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
+    if (!bin->FillHists(h)) return kFALSE; 
+    bin->CumulantsAccumulate(fCent);
+    i++;
   }
 
   return kTRUE;
@@ -294,6 +314,7 @@ void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
 
   PostData(2, fOutputList);
 
+  return;
 }
 //_____________________________________________________________________
 void AliForwardFlowTaskQC::Finalize()
@@ -302,25 +323,34 @@ void AliForwardFlowTaskQC::Finalize()
   // Finalize command, called by Terminate()
   //
 
-
   // Reinitiate vertex bins if Terminate is called separately!
   if (fBinsFMD.GetEntries() == 0) InitVertexBins();
 
   // Iterate over all vertex bins objects and finalize cumulants
   // calculations
-  TIter nextFMD(&fBinsFMD);
+  EndVtxBinList(fBinsFMD);
+  EndVtxBinList(fBinsSPD);
+
+  return;
+} 
+//_____________________________________________________________________
+void AliForwardFlowTaskQC::EndVtxBinList(const TList& list) const
+{
+  //
+  // Loop over VertexBin list and call terminate on each 
+  //
+  // Parameters:
+  //  list VertexBin list
+  //
+  TIter next(&list);
   VertexBin* bin = 0;
-  while ((bin = static_cast<VertexBin*>(nextFMD()))) {
+  while ((bin = static_cast<VertexBin*>(next()))) {
     bin->CumulantsTerminate(fSumList, fOutputList);
   }
-  TIter nextSPD(&fBinsSPD);
-  while ((bin = static_cast<VertexBin*>(nextSPD()))) {
-    bin->CumulantsTerminate(fSumList, fOutputList);
-  }
-
-} 
+  return;
+}
 // _____________________________________________________________________
-Bool_t AliForwardFlowTaskQC::AODCheck(const AliAODForwardMult* aodfm) 
+Bool_t AliForwardFlowTaskQC::CheckEvent(const AliAODForwardMult* aodfm) 
 {
   // 
   // Function to check that and AOD event meets the cuts
@@ -333,21 +363,63 @@ Bool_t AliForwardFlowTaskQC::AODCheck(const AliAODForwardMult* aodfm)
   //
 
   // First check for trigger
-  if (!aodfm->IsTriggerBits(AliAODForwardMult::kOffline)) return kFALSE;
+  if (!CheckTrigger(aodfm)) return kFALSE;
 
   // Then check for centrality
+  if (!GetCentrality(aodfm)) return kFALSE;
+
+  // And finally check for vertex
+  if (!GetVertex(aodfm)) return kFALSE;
+
+  return kTRUE;
+}
+// _____________________________________________________________________
+Bool_t AliForwardFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const 
+{
+  //
+  // Function to look for a trigger string in the event.
+  //
+  // Parameters: 
+  //  AliAODForwardMult: forward mult object with trigger and vertex info
+  //
+  // Returns true if offline trigger is present
+  //
+  return aodfm->IsTriggerBits(AliAODForwardMult::kOffline);
+}
+// _____________________________________________________________________
+Bool_t AliForwardFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm) 
+{
+  //
+  // Function to look get centrality of the event.
+  //
+  // Parameters: 
+  //  AliAODForwardMult: forward mult object with trigger and vertex info
+  //
+  // Returns true if centrality determination is present
+  //
   fCent = (Double_t)aodfm->GetCentrality();
   if (0. >= fCent || fCent >= 100.) return kFALSE;
   fHistCent->Fill(fCent);
 
-  // And finally check for vertex
-  fZvertex = aodfm->GetIpZ();
-  fHistVertexAll->Fill(fZvertex);
-  if (TMath::Abs(fZvertex) >= 10.) return kFALSE;
-  fHistVertexSel->Fill(fZvertex);
-
   return kTRUE;
+}
+// _____________________________________________________________________
+Bool_t AliForwardFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm) 
+{
+  //
+  // Function to look for vertex determination in the event.
+  //
+  // Parameters: 
+  //  AliAODForwardMult: forward mult object with trigger and vertex info
+  //
+  // Returns true if vertex is determined
+  //
+  fVtx = aodfm->GetIpZ();
+  fHistVertexAll->Fill(fVtx);
+  if (fVtx < fVtxAxis->GetXmin() || fVtx > fVtxAxis->GetXmax()) return kFALSE;
+  fHistVertexSel->Fill(fVtx);
 
+  return kTRUE;
 }
 //_____________________________________________________________________
 AliForwardFlowTaskQC::VertexBin::VertexBin()
@@ -355,28 +427,33 @@ AliForwardFlowTaskQC::VertexBin::VertexBin()
     fMoment(0),      // Flow moment for this vertexbin
     fVzMin(0),       // Vertex z-coordinate min
     fVzMax(0),       // Vertex z-coordinate max
-    fType(),         // Data type (FMD/SPD/FMDTR/SPDTR/MC)
+    fType(),         // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
+    fSymEta(1),      // Use forward-backward symmetry, if detector allows it
     fCumuRef(),      // Histogram for reference flow
     fCumuDiff(),     // Histogram for differential flow
     fCumuHist(),     // Sum histogram for cumulants
-    fdNdedpAcc()     // Diagnostics histogram to make acc. maps
+    fdNdedpAcc(),    // Diagnostics histogram to make acc. maps
+    fDebug()         // Debug level
 {
   //
   // Default constructor
   //
 }
 //_____________________________________________________________________
-AliForwardFlowTaskQC::VertexBin::VertexBin(const Int_t vLow, const Int_t vHigh, 
-                                           const Int_t moment, const Char_t* name)
+AliForwardFlowTaskQC::VertexBin::VertexBin(Int_t vLow, Int_t vHigh, 
+                                           UShort_t moment, TString name,
+                                           Bool_t sym)
   : TNamed("", ""),
     fMoment(moment),  // Flow moment for this vertexbin
     fVzMin(vLow),     // Vertex z-coordinate min
     fVzMax(vHigh),    // Vertex z-coordinate max
-    fType(name),      // Data type (FMD/SPD/FMDTR/SPDTR/MC)
+    fType(name),      // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
+    fSymEta(sym),     // Use forward-backward symmetry, if detector allows it
     fCumuRef(),       // Histogram for reference flow
     fCumuDiff(),      // Histogram for differential flow
     fCumuHist(),      // Sum histogram for cumulants
-    fdNdedpAcc()      // Diagnostics histogram to make acc. maps
+    fdNdedpAcc(),     // Diagnostics histogram to make acc. maps
+    fDebug(0)         // Debug level
 {
   //
   // Constructor
@@ -386,9 +463,14 @@ AliForwardFlowTaskQC::VertexBin::VertexBin(const Int_t vLow, const Int_t vHigh,
   //  vHigh: max z-coordinate
   //  moment: flow moment
   //  name: data type name (FMD/SPD/FMDTR/SPDTR/MC)
+  //  sym: data is symmetric in eta
   //
-  SetName(Form("%svertexBin%d_%d_%d", name, moment, vLow, vHigh));
-  SetTitle(Form("%svertexBin%d_%d_%d", name, moment, vLow, vHigh));
+  fType.ToUpper();
+
+  SetName(Form("%svertexBin%d_%d_%d", fType.Data(), moment, vLow, vHigh));
+  SetTitle(Form("%svertexBin%d_%d_%d", fType.Data(), moment, vLow, vHigh));
+
+  fDebug = AliAnalysisManager::GetAnalysisManager()->GetDebugLevel();
 
 }
 //_____________________________________________________________________
@@ -407,6 +489,7 @@ AliForwardFlowTaskQC::VertexBin::operator=(const AliForwardFlowTaskQC::VertexBin
   fCumuDiff     = o.fCumuDiff;
   fCumuHist     = o.fCumuHist;
   fdNdedpAcc    = o.fdNdedpAcc;
+  fDebug        = o.fDebug;
 
   return *this;
 }
@@ -421,23 +504,25 @@ void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist)
   //
 
   // First we try to find an outputlist for this vertexbin
-  TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d", fType, fVzMin, fVzMax));
+  TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d", fType.Data(), fVzMin, fVzMax));
 
   // If it doesn't exist we make one
   if (!list) {
     list = new TList();
-    list->SetName(Form("%svertex_%d_%d", fType, fVzMin, fVzMax));
+    list->SetName(Form("%svertex_%d_%d", fType.Data(), fVzMin, fVzMax));
     outputlist->Add(list);
   }
 
   // We initiate the reference histogram according to an acceptance correction map,
   // so we don't shift the SPD coverage within a reference bin
-  fCumuRef = new TH2D(Form("%s_v%d_%d_%d_ref", fType, fMoment, fVzMin, fVzMax),
-                        Form("%s_v%d_%d_%d_ref", fType, fMoment, fVzMin, fVzMax),
-                        24, -6., 6., 5, 0.5, 5.5);
+  // We start with many bins, to avoid memory problems. After checking for acc maps we
+  // rebin so those with full coverage only have 2 bins.
+  fCumuRef = new TH2D(Form("%s_v%d_%d_%d_ref", fType.Data(), fMoment, fVzMin, fVzMax),
+                        Form("%s_v%d_%d_%d_ref", fType.Data(), fMoment, fVzMin, fVzMax),
+                        48, -6., 6., 5, 0.5, 5.5);
   TFile acc("$ALICE_ROOT/PWGLF/FORWARD/corrections/FlowCorrections/FlowAccMap.root", "READ");
-  TH1D* accMap = (TH1D*)acc.Get(Form("%saccVertex_%d_%d", fType, fVzMin, fVzMax));
-  if (accMap) {
+  TH1D* accMap = (TH1D*)acc.Get(Form("%saccVertex_%d_%d", fType.Data(), fVzMin, fVzMax));
+  if (accMap && !fType.EqualTo("FMD")) {
     Int_t nBins = accMap->GetNbinsX();
     Double_t eta[48] = { 0. };
     Int_t n = 0;
@@ -445,32 +530,35 @@ void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist)
     Double_t prev = -1;
     for (Int_t i = 0; i < nBins; i++) {
       Double_t occ = accMap->GetBinContent(i+1);
-      if (prev != occ && (occ > 0.6 || occ == 0)) {
+      if (prev != occ && (((occ > 0.6 || occ == 0) && i*0.25-6 < 4) || ((occ == 0) && i*0.25-6 >= 4))) {
         eta[n] = i*0.25-6.;
         newOcc[n] = occ;
         n++;
-//        printf("eta: %f \t occ: %f \t Vertex: %d \n", eta[n-1], occ, fVzMin);
+        if (fDebug > 5) AliInfo(Form("eta: %f \t occ: %f \t Vertex: %d \n", eta[n-1], occ, fVzMin));
       }
       prev = occ;
     }
     eta[n] = 6.;
     fCumuRef->GetXaxis()->Set(n, eta);
+  } else {
+    fCumuRef->RebinX(24);
   }
   acc.Close();
 
   fCumuRef->Sumw2();
   //list->Add(fCumuRef);
 
   // We initiate the differential histogram
-  fCumuDiff = new TH2D(Form("%s_v%d_%d_%d_diff", fType, fMoment, fVzMin, fVzMax),
-                        Form("%s_v%d_%d_%d_diff", fType, fMoment, fVzMin, fVzMax),
-                        48, -6., 6., 5, 0.5, 5.5);
+  fCumuDiff = new TH2D(Form("%s_v%d_%d_%d_diff", fType.Data(), fMoment, fVzMin, fVzMax),
+                       Form("%s_v%d_%d_%d_diff", fType.Data(), fMoment, fVzMin, fVzMax),
+                       48, -6., 6., 5, 0.5, 5.5);
   fCumuDiff->Sumw2();
   //list->Add(fCumuDiff);
 
   // Initiate the cumulant sum histogram
-  fCumuHist = new TH3D(Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax),
-                       Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax),
+  fCumuHist = new TH3D(Form("%sv%d_vertex_%d_%d_cumu", fType.Data(), fMoment, fVzMin, fVzMax),
+                       Form("%sv%d_vertex_%d_%d_cumu", fType.Data(), fMoment, fVzMin, fVzMax),
                        48, -6., 6., 20, 0., 100., 26, 0.5, 26.5);
   fCumuHist->Sumw2();
 
@@ -482,35 +570,24 @@ void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist)
   if (!dList) AliFatal("No diagnostics list found, what kind of game are you running here?!?!");
 
   // Acceptance hists are shared over all moments
-  fdNdedpAcc = (TH2D*)dList->FindObject(Form("h%sdNdedpAcc_%d_%d", fType, fVzMin, fVzMax));
+  fdNdedpAcc = (TH2D*)dList->FindObject(Form("h%sdNdedpAcc_%d_%d", fType.Data(), fVzMin, fVzMax));
   if (!fdNdedpAcc) {
-    fdNdedpAcc = new TH2D(Form("h%sdNdedpAcc_%d_%d", fType, fVzMin, fVzMax), 
-                          Form("%s acceptance map for %d cm < v_{z} < %d cm", fType, fVzMin, fVzMax),
+    fdNdedpAcc = new TH2D(Form("h%sdNdedpAcc_%d_%d", fType.Data(), fVzMin, fVzMax), 
+                          Form("%s acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
                           48, -6, 6, 20, 0, TMath::TwoPi());
     fdNdedpAcc->Sumw2();
     dList->Add(fdNdedpAcc);
   }
+  TAxis* axis = (TAxis*)dList->FindObject(Form("axis%s_%d_%d", fType.Data(), fVzMin, fVzMax));
+  if (!axis) {
+    axis = fCumuRef->GetXaxis();
+    axis->SetName(Form("axis%s_%d_%d", fType.Data(), fVzMin, fVzMax));
+    dList->Add(fCumuRef);
+  }
 
-
-}
-//_____________________________________________________________________
-Bool_t AliForwardFlowTaskQC::VertexBin::CheckVertex(Double_t vz)
-{
-  //
-  // We check if this is the correct bin for the current event's vertex
-  //
-  // Parameters:
-  //  vZ: Current event vertex
-  //
-  // Returns false if out of range, true otherwise
-  //
-  if ((Double_t)fVzMin > vz) return kFALSE;
-  if ((Double_t)fVzMax <= vz) return kFALSE;
-
-  return kTRUE; 
 }
 //_____________________________________________________________________
-Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D* dNdetadphi) 
+Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(const TH2D& dNdetadphi) 
 {
   // 
   // Fill reference and differential eta-histograms
@@ -530,24 +607,29 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D* dNdetadphi)
   Double_t max = 0;
   Int_t nInAvg = 0;
   Int_t nBadBins = 0;
-  Int_t nBins = (dNdetadphi->GetNbinsX() * 6) / (fCumuDiff->GetNbinsX() * 5);
+  Int_t nBins = (dNdetadphi.GetNbinsX() * 6) / (fCumuDiff->GetNbinsX() * 5);
   Int_t nInBin = 0;
-  Int_t nCurBin = 0, nPrevBin = 1;
+  Int_t nCurBin = 0, nPrevBin = 0;
 
   // Then we loop over the input and calculate sum cos(k*n*phi)
   // and fill it in the reference and differential histograms
   Double_t eta, phi, weight;
   Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0;
-  for (Int_t etaBin = 1; etaBin <= dNdetadphi->GetNbinsX(); etaBin++) {
-    eta = dNdetadphi->GetXaxis()->GetBinCenter(etaBin);
+  for (Int_t etaBin = 1; etaBin <= dNdetadphi.GetNbinsX(); etaBin++) {
+    eta = dNdetadphi.GetXaxis()->GetBinCenter(etaBin);
     nCurBin = fCumuDiff->GetXaxis()->FindBin(eta);
     // If we have moved to a new bin in the flow hist, and less than half the eta
     // region has been covered by it we cut it away.
+    if (!nPrevBin) nPrevBin = nCurBin;
     if (nCurBin != nPrevBin) {
-      if (nInBin <= nBins/2) {
-       for (Int_t pBin = 1; pBin <= fCumuDiff->GetNbinsY(); pBin++) {
-         fCumuDiff->SetBinContent(nPrevBin, pBin, 0);
-         fCumuDiff->SetBinError(nPrevBin, pBin, 0);
+      if (nInBin < nBins) {
+       for (Int_t qBin = 1; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
+         Double_t removeContent = fCumuDiff->GetBinContent(nPrevBin, qBin);
+         Double_t removeEta = fCumuDiff->GetXaxis()->GetBinCenter(nPrevBin);
+         fCumuRef->Fill(removeEta, qBin, -removeContent);
+         if (fSymEta) fCumuRef->Fill(-1.*removeEta, qBin, -removeContent);
+         fCumuDiff->SetBinContent(nPrevBin, qBin, 0);
+         fCumuDiff->SetBinError(nPrevBin, qBin, 0);
        }
       }
       nInBin = 0;
@@ -557,9 +639,12 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D* dNdetadphi)
       max = 0;
     }
     Bool_t data = kFALSE;
-    for (Int_t phiBin = 1; phiBin <= dNdetadphi->GetNbinsY(); phiBin++) {
-      phi = dNdetadphi->GetYaxis()->GetBinCenter(phiBin);
-      weight = dNdetadphi->GetBinContent(etaBin, phiBin);
+    for (Int_t phiBin = 1; phiBin <= dNdetadphi.GetNbinsY(); phiBin++) {
+//      if (fType == "FMD" && eta < 0 && phiBin == 11) continue;
+//      if (fType == "FMD" && eta > 0 && phiBin == 20) continue;
+
+      phi = dNdetadphi.GetYaxis()->GetBinCenter(phiBin);
+      weight = dNdetadphi.GetBinContent(etaBin, phiBin);
       if (!weight) continue;
       if (!data) data = kTRUE;
       // We calculate the running average Nch per. bin
@@ -589,10 +674,19 @@ Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D* dNdetadphi)
 
       // Fill acc. map
       fdNdedpAcc->Fill(eta, phi, weight);
+
+      if (!fSymEta) continue;
+
+      fCumuRef->Fill(-eta, kHmult, weight);
+      fCumuRef->Fill(-eta, kHQnRe, dQnRe);
+      fCumuRef->Fill(-eta, kHQnIm, dQnIm);
+      fCumuRef->Fill(-eta, kHQ2nRe, dQ2nRe);
+      fCumuRef->Fill(-eta, kHQ2nIm, dQ2nIm);
+
     }
     if (data) {
       nInBin++;
-      if (max > 2*runAvg) nBadBins++;
+//      if (max > 2*runAvg) nBadBins++;
     }
     // If there are too many bad bins we throw the event away!
     if (nBadBins > 3) return kFALSE;
@@ -623,7 +717,7 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
   Double_t multi = 0, multp = 0, mp = 0, mq = 0;
   Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
   Int_t refEtaBin = 0;
-  
+
   // We loop over the data 1 time!
   for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
     eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
@@ -779,23 +873,23 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* o
 
   // Re-find cumulants hist if Terminate is called separately
   if (!fCumuHist) {
-    TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d", fType, fVzMin, fVzMax));
-    fCumuHist = (TH3D*)list->FindObject(Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax));
+    TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d", fType.Data(), fVzMin, fVzMax));
+    fCumuHist = (TH3D*)list->FindObject(Form("%sv%d_vertex_%d_%d_cumu", fType.Data(), fMoment, fVzMin, fVzMax));
   }
 
   // Create result profiles
-  TProfile2D* cumu2 = (TProfile2D*)outlist->FindObject(Form("%sQC2_v%d_unCorr", fType, fMoment)); 
-  TProfile2D* cumu4 = (TProfile2D*)outlist->FindObject(Form("%sQC4_v%d_unCorr", fType, fMoment)); 
+  TProfile2D* cumu2 = (TProfile2D*)outlist->FindObject(Form("%sQC2_v%d_unCorr", fType.Data(), fMoment)); 
+  TProfile2D* cumu4 = (TProfile2D*)outlist->FindObject(Form("%sQC4_v%d_unCorr", fType.Data(), fMoment)); 
   if (!cumu2) {
-    cumu2 = new TProfile2D(Form("%sQC2_v%d_unCorr", fType, fMoment),
-                           Form("%sQC2_v%d_unCorr", fType, fMoment),
+    cumu2 = new TProfile2D(Form("%sQC2_v%d_unCorr", fType.Data(), fMoment),
+                           Form("%sQC2_v%d_unCorr", fType.Data(), fMoment),
              fCumuHist->GetNbinsX(), fCumuHist->GetXaxis()->GetXmin(), fCumuHist->GetXaxis()->GetXmax(), 
              fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
     outlist->Add(cumu2);
   }
   if (!cumu4) {
-    cumu4 = new TProfile2D(Form("%sQC4_v%d_unCorr", fType, fMoment),
-                           Form("%sQC4_v%d_unCorr", fType, fMoment),
+    cumu4 = new TProfile2D(Form("%sQC4_v%d_unCorr", fType.Data(), fMoment),
+                           Form("%sQC4_v%d_unCorr", fType.Data(), fMoment),
              fCumuHist->GetNbinsX(), fCumuHist->GetXaxis()->GetXmin(), fCumuHist->GetXaxis()->GetXmax(), 
              fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
     outlist->Add(cumu4);
@@ -817,6 +911,7 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* o
   for (Int_t cBin = 1; cBin <= fCumuHist->GetNbinsY(); cBin++) {
     Double_t cent = fCumuHist->GetYaxis()->GetBinCenter(cBin);
     Double_t nEv = 0;
+    if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), fMoment, cent));
     // Eta loop
     for (Int_t etaBin = 1; etaBin <= fCumuHist->GetNbinsX(); etaBin++) {
       Double_t eta = fCumuHist->GetXaxis()->GetBinCenter(etaBin);
@@ -832,11 +927,31 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* o
       sinP1nPhi /= mult;
       two = w2Two / w2;
       qc2 = two - TMath::Power(cosP1nPhi, 2) - TMath::Power(sinP1nPhi, 2);
-      if (qc2 <= 0) continue;
+      if (qc2 <= 0) { 
+       if (fDebug > 0) AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping", fType.Data(), fMoment, qc2, eta, cent));
+       continue;
+      }
       vnTwo = TMath::Sqrt(qc2);
  //     if (!TMath::IsNaN(vnTwo*mult)) 
  //       cumu2->Fill(eta, cent, vnTwo, fCumuHist->GetBinContent(0,cBin,0)); 
 
+      // 2-particle differential flow
+      w2pTwoPrime = fCumuHist->GetBinContent(etaBin, cBin, kw2two);
+      w2p = fCumuHist->GetBinContent(etaBin, cBin, kw2);
+      mp = fCumuHist->GetBinContent(etaBin, cBin, kmp);
+      if (!w2p || !mp) continue;
+      cosP1nPsi = fCumuHist->GetBinContent(etaBin, cBin, kpnRe);
+      sinP1nPsi = fCumuHist->GetBinContent(etaBin, cBin, kpnIm);
+
+      cosP1nPsi /= mp;
+      sinP1nPsi /= mp;
+      twoPrime = w2pTwoPrime / w2p;
+      qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi;
+
+      vnTwoDiff = qc2Prime / TMath::Sqrt(qc2);
+      if (!TMath::IsNaN(vnTwoDiff*mp)) cumu2->Fill(eta, cent, vnTwoDiff);
+      if (fDebug > 1) AliInfo(Form("%s: v_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f", fType.Data(), fMoment, vnTwoDiff, eta, cent));
+
       // 4-particle reference flow
       w4Four = fCumuHist->GetBinContent(etaBin, cBin, kW4Four);
       w4 = fCumuHist->GetBinContent(etaBin, cBin, kW4);
@@ -860,27 +975,14 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* o
          + 8.*two*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.))
          - 6.*TMath::Power((TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)),2.);
       
-      if (qc4 >= 0) continue;
+      if (qc4 >= 0) {
+       if (fDebug > 0) AliInfo(Form("%s: QC_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f - skipping", fType.Data(), fMoment, qc2, eta, cent));
+       continue;
+      }
       vnFour = TMath::Power(-qc4, 0.25);
  //     if (!TMath::IsNaN(vnFour*mult)) 
  //         cumu4->Fill(eta, cent, vnFour, fCumuHist->GetBinContent(0,cBin,0));
 
-      // 2-particle differential flow
-      w2pTwoPrime = fCumuHist->GetBinContent(etaBin, cBin, kw2two);
-      w2p = fCumuHist->GetBinContent(etaBin, cBin, kw2);
-      mp = fCumuHist->GetBinContent(etaBin, cBin, kmp);
-      if (!w2p || !mp) continue;
-      cosP1nPsi = fCumuHist->GetBinContent(etaBin, cBin, kpnRe);
-      sinP1nPsi = fCumuHist->GetBinContent(etaBin, cBin, kpnIm);
-
-      cosP1nPsi /= mp;
-      sinP1nPsi /= mp;
-      twoPrime = w2pTwoPrime / w2p;
-      qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi;
-
-      vnTwoDiff = qc2Prime / TMath::Sqrt(qc2);
-      if (!TMath::IsNaN(vnTwoDiff*mp)) cumu2->Fill(eta, cent, vnTwoDiff, fCumuHist->GetBinContent(0,cBin,0));
-
       // 4-particle differential flow
       w4pFourPrime = fCumuHist->GetBinContent(etaBin, cBin, kw4four);
       w4p = fCumuHist->GetBinContent(etaBin, cBin, kw4);
@@ -923,7 +1025,8 @@ void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* o
                 * (sinP1nPsi*cosP1nPhi+cosP1nPsi*sinP1nPhi);
 
       vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
-      if (!TMath::IsNaN(vnFourDiff*mp)) cumu4->Fill(eta, cent, vnFourDiff, fCumuHist->GetBinContent(0,cBin,0));
+      if (!TMath::IsNaN(vnFourDiff*mp) && vnFourDiff > 0) cumu4->Fill(eta, cent, vnFourDiff);
+      if (fDebug > 1) AliInfo(Form("%s: v_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f", fType.Data(), fMoment, vnFourDiff, eta, cent));
     } // End of eta loop
     // Number of events:
     nEv += fCumuHist->GetBinContent(0,cBin,0);
index 2ea15d3acb28a1b1f413a496dcded79a22a9b8bc..5f7bbda45e9c044748fdf64f2ce347899745b90c 100644 (file)
  * @ingroup pwglf_forward_flow
  */
 #include "AliAnalysisTaskSE.h"
+#include "TString.h"
 class AliAODForwardMult;
 class TH1D;
 class TH2D;
 class TH3D;
+class TAxis;
 
  /**
  * @defgroup pwglf_forward_tasks_flow Flow tasks 
@@ -79,22 +81,65 @@ public:
    */
   virtual void Terminate(Option_t *option);
   /* @} */
+  /**
+   * Loops of vertex bins in list and runs analysis on those for current vertex
+   *
+   * @param list List of vertex bins
+   * @param h dN/detadphi histogram
+   * @param vtx Current vertex bin
+   *
+   * @return true on success
+   */
+  Bool_t FillVtxBinList(const TList& list, const TH2D& h, Int_t vtx) const;
+  /**
+   * Loops over VertexBin list and calls terminate on each
+   *
+   * @param list VertexBin list
+   */
+  void EndVtxBinList(const TList& list) const;
   /**
    * Returns the outputlist
    * 
    * @return TList* 
-  */
+   */
   TList* GetOutputList() { return fOutputList; }
- /**
-  * Check AODForwardMult object for trigger, vertex and centrality
-  * returns true if event is OK
-  * 
-  * @param const aodfm
-  * 
-  * @return Bool_t 
-  */
-  Bool_t AODCheck(const AliAODForwardMult* aodfm);
-   /*
+  /**
+   * Check AODevent object for trigger, vertex and centrality
+   * returns true if event is OK
+   *
+   * @param aodfm AliAODForwardMultObject
+   * 
+   * @return Bool_t 
+   */
+  Bool_t CheckEvent(const AliAODForwardMult* aodfm);
+  /**
+   * Check trigger from AODForwardMult object
+   * returns true if offline trigger is present
+   *
+   * @param aodfm AliAODForwardMultObject
+   * 
+   * @return Bool_t 
+   */
+  virtual Bool_t CheckTrigger(const AliAODForwardMult* aodfm) const;
+  /**
+   * Check for centrality in AliAODForwardMult object, 
+   * if present return true - also sets fCent value
+   *
+   * @param aodfm AliAODForwardMultObject
+   * 
+   * @return Bool_t 
+   */
+  virtual Bool_t GetCentrality(const AliAODForwardMult* aodfm);
+  /*
+   * Check for vertex in AliAODForwardMult
+   * returns true if in range of fVtxAXis, also sets fVtx value
+   *
+   * @param aodfm AliAODForwardMultObject
+   * 
+   * @return Bool_t 
+   */
+  virtual Bool_t GetVertex(const AliAODForwardMult* aodfm);
+  /**
    * Set which harmonics to calculate. @f$ v_{1}@f$ to @f$ v_{4}@f$ is
    * available and calculated as default
    * 
@@ -111,6 +156,14 @@ public:
                      Bool_t v3 = kTRUE, Bool_t v4 = kTRUE,
                      Bool_t v5 = kTRUE, Bool_t v6 = kTRUE) { 
     fv[1] = v1; fv[2] = v2; fv[3] = v3; fv[4] = v4; fv[5] = v5; fv[6] = v6;}
+   /*
+   * Set non-default vertex binning and range
+   *
+   * @param axis Use this vtx axis
+   *
+   * @return void
+   */
+  void SetVertexAxis(TAxis* axis) { fVtxAxis = axis; }
  /**
   * Nested class to handle cumulant calculations in vertex bins
   */
@@ -121,102 +174,98 @@ public:
      * Constructor
      */
     VertexBin();
-   /**
-    * Constructor
-    * 
-    * @param vLow Min vertex z-coordinate
-    * @param vHigh Max vertex z-coordinate
-    * @param moment Flow moment
-    * @param type Data type (FMD/SPD/FMDTR/SPDTR/MC)
-    */
-    VertexBin(const Int_t vLow, const Int_t vHigh, 
-              const Int_t moment, const Char_t* type);
-   /**
-    * Copy constructor 
-    * 
-    * @param o Object to copy from
-    * 
-    * @return VertexBin
-    */
+    /**
+     * Constructor
+     * 
+     * @param vLow Min vertex z-coordinate
+     * @param vHigh Max vertex z-coordinate
+     * @param moment Flow moment
+     * @param type Data type (FMD/SPD/FMDTR/SPDTR/MC)
+     * @parma sym Data is symmetric in eta
+     */
+    VertexBin(Int_t vLow, Int_t vHigh, 
+              UShort_t moment, TString type,
+              Bool_t sym = kTRUE);
+    /**
+     * Copy constructor 
+     * 
+     * @param o Object to copy from
+     * 
+     * @return VertexBin
+     */
     VertexBin(const VertexBin& o);
-   /**
-    * Assignment operator 
-    * 
-    * @param VertexBin&
-    * 
-    * @return VertexBin& 
-    */
+    /**
+     * Assignment operator 
+     
+     * @param VertexBin&
+     
+     * @return VertexBin& 
+     */
     VertexBin& operator=(const VertexBin&);
-   /**
-    * Destructor 
-    */
+    /**
+     * Destructor 
+     */
     ~VertexBin(){}
-   /**
-    * Add vertex bin output to list
-    * 
-    * @param list Histograms are added to this list
-    * 
-    * @return void 
-    */
+    /**
+     * Add vertex bin output to list
+     
+     * @param list Histograms are added to this list
+     
+     * @return void 
+     */
     virtual void AddOutput(TList* list);
-   /**
-    * Check if vertex vZ is in range for this bin
-    * 
-    * @param vZ z-coordinate of vertex to check
-    * 
-    * @return Bool_t 
-    */
-    Bool_t CheckVertex(Double_t vZ);
-  /**
-    * Fill reference and differential flow histograms for analysis
-    *
-    * @param dNdetadphi 2D data histogram
-    *
-    * @return false if bad event (det. hotspot)
-    */
-    Bool_t FillHists(TH2D* dNdetadphi);
-   /**
-    * Do cumulants calculations for current event with 
-    * centrality cent
-    * 
-    * @param cent Event centrality
-    * 
-    * @return void 
-    */
+    /**
+     * Fill reference and differential flow histograms for analysis
+     *
+     * @param dNdetadphi 2D data histogram
+     *
+     * @return false if bad event (det. hotspot)
+     */
+    Bool_t FillHists(const TH2D& dNdetadphi);
+    /**
+     * Do cumulants calculations for current event with 
+     * centrality cent
+     * 
+     * @param cent Event centrality
+     * 
+     * @return void 
+     */
     void CumulantsAccumulate(Double_t cent);
-   /**
-    * Finish cumulants calculations. Takes input and
-    * output lists in case Terminate is called separately
-    * 
-    * @param inlist List with input histograms
-    * @param outlist List with output histograms
-    * 
-    * @return void 
-    */
+    /**
+     * Finish cumulants calculations. Takes input and
+     * output lists in case Terminate is called separately
+     
+     * @param inlist List with input histograms
+     * @param outlist List with output histograms
+     
+     * @return void 
+     */
     void CumulantsTerminate(TList* inlist, TList* outlist);
 
     protected:
-   /*
-    * Enumeration for ref/diff histograms
-    */
+    /*
+     * Enumeration for ref/diff histograms
+     */
     enum { kHmult = 1, kHQnRe, kHQnIm, kHQ2nRe, kHQ2nIm };
-   /*
-    * Enumeration for cumulant histogram
-    */
+    /*
+     * Enumeration for cumulant histogram
+     */
     enum { kW2Two = 1, kW2, kW4Four, kW4, kQnRe, kQnIm, kM,
        kCosphi1phi2, kSinphi1phi2, kCosphi1phi2phi3m, kSinphi1phi2phi3m, kMm1m2, 
        kw2two, kw2, kw4four, kw4, kpnRe, kpnIm, kmp, 
        kCospsi1phi2, kSinpsi1phi2, kCospsi1phi2phi3m, kSinpsi1phi2phi3m,
        kmpmq, kCospsi1phi2phi3p, kSinpsi1phi2phi3p };
 
-    const Int_t   fMoment;        // flow moment 
-    const Int_t   fVzMin;         // z-vertex min
-    const Int_t   fVzMax;         // z-vertex max
-    const Char_t* fType;          // data type
-    TH2D*         fCumuRef;       // histogram for reference flow
-    TH2D*         fCumuDiff;      // histogram for differential flow
-    TH3D*         fCumuHist;      // histogram for cumulants calculations
-    TH2D*         fdNdedpAcc;     // Diagnostics histogram to make acc. maps
+    const UShort_t fMoment;        // flow moment 
+    const Int_t    fVzMin;         // z-vertex min must be in whole [cm]
+    const Int_t    fVzMax;         // z-vertex max must be in whoe [cm]
+    TString        fType;          // data type
+    const Bool_t   fSymEta;        // Use forward-backward symmetry, if detector allows it
+    TH2D*          fCumuRef;       // histogram for reference flow
+    TH2D*          fCumuDiff;      // histogram for differential flow
+    TH3D*          fCumuHist;      // histogram for cumulants calculations
+    TH2D*          fdNdedpAcc;     // Diagnostics histogram to make acc. maps
+    UShort_t       fDebug;         // Debug flag
 
     ClassDef(VertexBin, 1); // object for cumulants ananlysis in FMD
   };
@@ -251,13 +300,14 @@ protected:
    */
   virtual void Finalize();
 
+  TAxis*         fVtxAxis;        //  Axis to control vertex binning
   TList          fBinsFMD;        //  list with FMD VertexBin objects 
   TList          fBinsSPD;        //  list with SPD VertexBin objects
   TList*         fSumList;        //  sum list
   TList*         fOutputList;     //  Output list
   AliAODEvent*   fAOD;            //  AOD event
   Bool_t         fv[7];           //  Calculate v_{n} flag
-  Float_t       fZvertex;        //  Z vertex bin
+  Float_t       fVtx;        //  Z vertex bin
   Double_t       fCent;           //  Centrality
   TH1D*          fHistCent;       //  Diagnostics hist for centrality
   TH1D*          fHistVertexSel;  //  Diagnostics hist for selected vertices
index 55bc30e3fb8a93f32274c2ba5e089ca7b518034a..e277070115eeda693df73c3d4eff312c4a2e4200 100644 (file)
@@ -31,6 +31,7 @@ AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC() :
   fAlicePt4th3040(),     // Alice QC4 vs. pT data points
   fAlicePt4th4050(),     // Alice QC4 vs. pT data points
   fImpactParToCent(),    // Impact parameter to centrality graph
+  fUseImpactPar(0),      // Use impact par for centrality
   fAddFlow(0),           // Add flow to MC truth
   fAddType(0),           // Add type of flow to MC truth
   fAddOrder(0)           // Add order of flow to MC truth        
@@ -50,6 +51,7 @@ AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name) :
   fAlicePt4th3040(),     // Alice QC4 vs. pT data points
   fAlicePt4th4050(),     // Alice QC4 vs. pT data points
   fImpactParToCent(),    // Impact parameter to centrality graph
+  fUseImpactPar(0),      // Use impact par for centrality
   fAddFlow(0),           // Add flow to MC truth
   fAddType(0),           // Add type of flow to MC truth
   fAddOrder(0)           // Add order of flow to MC truth        
@@ -92,8 +94,10 @@ AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name) :
   Int_t nPointsCumulant4th4050ALICE = sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);   
   fAlicePt4th4050 = new TGraph(nPointsCumulant4th4050ALICE, xCumulant4th4050ALICE, yCumulant4th4050ALICE);
 
-  Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,11.565,12.575,13.515,16.679};
-  Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
+//  Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,11.565,12.575,13.515,16.679};
+//  Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
+  Double_t impactParam[] = {0., 3.72, 5.23, 7.31, 8.88, 10.20, 11.38, 12.47, 13.50, 14.51, 16.679};
+  Double_t centrality[] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 100.};
 
   Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
   fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
@@ -111,6 +115,7 @@ AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o)
   fAlicePt4th3040(o.fAlicePt4th3040),    // Alice QC4 vs. pT data points
   fAlicePt4th4050(o.fAlicePt4th4050),    // Alice QC4 vs. pT data points
   fImpactParToCent(o.fImpactParToCent),  // Impact parameter to centrality graph
+  fUseImpactPar(o.fUseImpactPar),        // Use impact par for centrality
   fAddFlow(o.fAddFlow),                  // Add flow to MC truth
   fAddType(o.fAddType),                  // Add type of flow to MC truth
   fAddOrder(o.fAddOrder)                 // Add order of flow to MC truth        
@@ -134,6 +139,7 @@ AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
   fAlicePt4th3040  = o.fAlicePt4th3040;
   fAlicePt4th4050  = o.fAlicePt4th4050;
   fImpactParToCent = o.fImpactParToCent;
+  fUseImpactPar    = o.fUseImpactPar;
   fAddFlow         = o.fAddFlow;
   fAddType         = o.fAddType;
   fAddOrder        = o.fAddOrder;
@@ -147,12 +153,12 @@ void AliForwardMCFlowTaskQC::InitVertexBins()
   //
   AliForwardFlowTaskQC::InitVertexBins();
 
-  for(Int_t n = 1; n <= 6; n++) {
+  for(UShort_t n = 1; n <= 6; n++) {
     if (!fv[n]) continue;
-    for (Int_t v = -10; v < 10; v++) {
-      fBinsFMDTR.Add(new VertexBin(v, v+1, n, "FMDTR"));
-      fBinsSPDTR.Add(new VertexBin(v, v+1, n, "SPDTR"));
-      fBinsMC.Add(new VertexBin(v, v+1, n, "MC"));
+      for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
+      fBinsFMDTR.Add(new VertexBin(fVtxAxis->GetBinLowEdge(v), fVtxAxis->GetBinUpEdge(v), n, "FMDTR"));
+      fBinsSPDTR.Add(new VertexBin(fVtxAxis->GetBinLowEdge(v), fVtxAxis->GetBinUpEdge(v), n, "SPDTR"));
+      fBinsMC.Add(new VertexBin(fVtxAxis->GetBinLowEdge(v), fVtxAxis->GetBinUpEdge(v), n, "MC"));
     }
   }
 
@@ -190,41 +196,22 @@ Bool_t AliForwardMCFlowTaskQC::Analyze()
   if (!AliForwardFlowTaskQC::Analyze()) return kFALSE;
 
   // Run analysis on trackrefs from FMD and SPD
-  AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
-  if (!aodfmult) return kFALSE;
-  TH2D fmdTRdNdetadphi = aodfmult->GetHistogram();
-
-  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
-  if (!aodcmult) return kFALSE;
-  TH2D spdTRdNdetadphi = aodcmult->GetHistogram();
+  const AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
+  const AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
+  if (!aodfmult || !aodcmult) return kFALSE;
   
-  TIter nextFMDTR(&fBinsFMDTR);
-  VertexBin* bin = 0;
-  while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
-    if (bin->CheckVertex(fZvertex)) {
-      if (!bin->FillHists(&fmdTRdNdetadphi)) return kFALSE;
-      bin->CumulantsAccumulate(fCent);
-    }
-  }
+  // if objects are present, get histograms
+  const TH2D& fmdTRdNdetadphi = aodfmult->GetHistogram();
+  const TH2D& spdTRdNdetadphi = aodcmult->GetHistogram();
 
-  TIter nextSPDTR(&fBinsSPDTR);
-  while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
-    if (bin->CheckVertex(fZvertex)) {
-      if (!bin->FillHists(&spdTRdNdetadphi)) return kFALSE;
-      bin->CumulantsAccumulate(fCent);
-    }
-  }
+  // Run analysis on tr refs
+  Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
+  if (!FillVtxBinList(fBinsFMDTR, fmdTRdNdetadphi, vtx)) return kFALSE;
+  if (!FillVtxBinList(fBinsSPDTR, spdTRdNdetadphi, vtx)) return kFALSE;
 
   // Run analysis on MC branch
   if (!LoopAODMC()) return kFALSE;
-
-  TIter nextMC(&fBinsMC);
-  while ((bin = static_cast<VertexBin*>(nextMC()))) {
-    if (bin->CheckVertex(fZvertex)) {
-      if (!bin->FillHists(&fdNdedpMC)) return kFALSE;
-      bin->CumulantsAccumulate(fCent);
-    }
-  }
+  if (!FillVtxBinList(fBinsMC, fdNdedpMC, vtx)) return kFALSE;
 
   return kTRUE;
 }
@@ -236,48 +223,30 @@ void AliForwardMCFlowTaskQC::Finalize()
   //
   AliForwardFlowTaskQC::Finalize();
 
-  TIter nextFMDTR(&fBinsFMDTR);
-  VertexBin* bin = 0;
-  while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
-    bin->CumulantsTerminate(fSumList, fOutputList);
-  }
-  TIter nextSPDTR(&fBinsSPDTR);
-  while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
-    bin->CumulantsTerminate(fSumList, fOutputList);
-  }
-  TIter nextMC(&fBinsMC);
-  while ((bin = static_cast<VertexBin*>(nextMC()))) {
-    bin->CumulantsTerminate(fSumList, fOutputList);
-  }
-
-  TProfile2D* fmdHist = 0;
-  TProfile2D* spdHist = 0;
-  TProfile2D* mcHist = 0;
-
-  for (Int_t i = 2; i <= 4; i += 2) {
-    for (Int_t n = 1; n <= 6; n++) {
-      if (!fv[n]) continue;
-      fmdHist = (TProfile2D*)fOutputList->FindObject(Form("FMDQC%d_v%d_unCorr", i, n))
-                 ->Clone(Form("FMDQC%d_v%d_Correction", i, n));
-      spdHist = (TProfile2D*)fOutputList->FindObject(Form("SPDQC%d_v%d_unCorr", i, n))
-                 ->Clone(Form("SPDQC%d_v%d_Correction", i, n));
-      mcHist = (TProfile2D*)fOutputList->FindObject(Form("MCQC%d_v%d_unCorr", i, n));
-     
-      if (!fmdHist || !spdHist || !mcHist) {
-       AliError(Form("Histogram missing, correction object not created for v%d", n));
-       continue;
-      }
+  EndVtxBinList(fBinsFMDTR);
+  EndVtxBinList(fBinsSPDTR);
+  EndVtxBinList(fBinsMC);
 
-      fmdHist->Divide(mcHist);
-      spdHist->Divide(mcHist);
-      fmdHist->SetTitle(Form("FMD QC{%d} v_{%d} Correction Object", i, n));
-      fmdHist->SetTitle(Form("SPD QC{%d} v_{%d} Correction Object", i, n));
-
-      fOutputList->Add(fmdHist);
-      fOutputList->Add(spdHist);
-    }
+  return;
+}
+//_____________________________________________________________________
+Bool_t AliForwardMCFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm)
+{
+  // 
+  // Function to use centrality parametrization from impact parameter
+  // if flag is not set call AliForwardFlowTaskQC::GetCentrality
+  //
+  // Parameters:
+  //  AliAODForwardMult: forward mult object with trigger and vertex info
+  //
+  // Returns true when centrality is set.
+  //
+  if (fUseImpactPar) {
+    fCent = GetCentFromB();
+    fHistCent->Fill(fCent);
+    return kTRUE;
   }
-
+  else  return AliForwardFlowTaskQC::GetCentrality(aodfm);
 }
 //_____________________________________________________________________
 Bool_t AliForwardMCFlowTaskQC::LoopAODMC()  
@@ -319,7 +288,7 @@ Bool_t AliForwardMCFlowTaskQC::LoopAODMC()
     Double_t pT = particle->Pt();
     Double_t eta = particle->Eta();
     Double_t phi = particle->Phi();
-    if (TMath::Abs(eta) < 6.) {
+    if (eta > -4. && eta < 5.) {
       // Add flow if it is in the argument
       if (fAddFlow.Length() > 1) {
         if (fAddFlow.Contains("pt"))
@@ -420,7 +389,6 @@ Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
 {
   //
   // Get centrality from MC impact parameter.
-  // Values taken from: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
   //
   Double_t cent = -1.;
   Double_t b = -1.;
index 84cf9496f8bd8eddea34c2b37f281c0d500d8321..75cdec0701d6354832b91659a9d819eaff3234ab 100644 (file)
@@ -52,18 +52,23 @@ public:
    * Destructor 
    */
   virtual ~AliForwardMCFlowTaskQC() {}
-  /** 
-   * @{ 
-   * @name Task interface methods 
-   */
   /**
-   * Loop over AliAODMCParticle branch object and fill d^2N/detadphi histograms
-   * add flow if arguments are set
+   * Check for centrality in AliAODForwardMult object, 
+   * if present return true - also sets fCent value
+   * can be used to get centrality from impact parameter
+   *
+   * @param aodfm AliAODForwardMultObject
    * 
-   * @return true on success
+   * @return Bool_t 
    */
-  Bool_t LoopAODMC();
-    /*
+  virtual Bool_t GetCentrality(const AliAODForwardMult* aodfm);
+  /*
+   * Set use parametrization from impact parameter for centrality
+   *
+   * @param use Use impact par
+   */
+  void SetUseImpactParameter(Bool_t use) { fUseImpactPar = use; }
+  /*
    * Set string to add flow to MC truth particles
    *
    * @param type String
@@ -111,6 +116,13 @@ protected:
    * Finalize analysis
    */
   void Finalize();
+  /**
+   * Loop over AliAODMCParticle branch object and fill d^2N/detadphi histograms
+   * add flow if arguments are set
+   * 
+   * @return true on success
+   */
+  Bool_t LoopAODMC();
   /**
    * Add pt dependent flow factor
    *
@@ -146,6 +158,7 @@ protected:
   TGraph*       fAlicePt4th3040;    //  Parametrization of ALICE QC4 vs. pT data
   TGraph*       fAlicePt4th4050;    //  Parametrization of ALICE QC4 vs. pT data
   TGraph*       fImpactParToCent;   //  Parametrization of b to centrality datapoints
+  Bool_t       fUseImpactPar;      //  Flag to use impact parameter for cent
   TString       fAddFlow;           //  Add flow string
   Int_t         fAddType;           //  Add flow type #
   Int_t         fAddOrder;          //  Add flow order
index 424847aea845f80801f4990845fc777692c15d33..0591a7916373847dfcf80c4c896d784616fa53d8 100644 (file)
  * @par Outputs: 
  * - 
  *
- * @ingroup pwg2_forward_flow
+ * @ingroup pwglf_forward_flow
  */
-void MakeFlow(TString data      = "", 
-             Int_t   nEvents   = 0
+void MakeFlow(const char* data  = "", 
+             Int_t   nEvents   = -1
              TString type      = "", 
               Bool_t  mc        = kFALSE,
              const char* name  = 0,
@@ -43,12 +43,32 @@ void MakeFlow(TString data      = "",
     gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/trains/TrainSetup.C+");
     gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C+");
 
-    MakeFlowTrain t(name, type.Data(), mc, addFlow.Data(), addFType, addFOrder, false);
-    t.SetDataDir(data.Data());
-    t.SetDataSet("");
-    t.SetProofServer(Form("workers=%d", proof));
+    char* tok = strtok(name, " ,:");
+    TString server = "", dataset = "", datadir = "", jobname = "";
+    if (tok[0] == NULL) Fatal("Input name invalid!");
+    if (strcmp(tok, "hehi00.nbi.dk") == 0) {
+      server = tok;
+      dataset = data;
+      datadir = "";
+      tok = strtok(NULL, " ,:");
+      jobname = tok;
+    } else {
+      server = Form("workers=%d", proof);
+      dataset = "";
+      datadir = data;
+      jobname = tok;
+    }
+    printf("Making flow train on server: %s\t with datadir: %s\t on dataset: %s\n", 
+           server.Data(), datadir.Data(), dataset.Data());
+
+    MakeFlowTrain t(jobname.Data(), type.Data(), mc, addFlow.Data(), addFType, addFOrder, false);
+    t.SetProofServer(server.Data());
+    t.SetDataDir(datadir.Data());
+    t.SetDataSet(dataset.Data());
+    t.SetDebugLevel(1);
     t.SetUseGDB(gdb);
-    t.Run(proof > 0 ? "proof" : "local", "full", nEvents, proof > 0);
+    t.Run("proof", "full", nEvents, proof > 0);
+
     return;
   }
 
@@ -58,12 +78,12 @@ void MakeFlow(TString data      = "",
                           gROOT->GetMacroPath()));
 
   // --- Add to chain either AOD ------------------------------------
-  if (data.IsNull()) {
+  if (data == '\0') {
     AliError("You didn't add a data file");
     return;
   }
   gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/MakeChain.C");
-  TChain* chain = MakeChain("AOD", data.Data(), true);
+  TChain* chain = MakeChain("AOD", data, true);
   // If 0 or less events is select, choose all 
   if (nEvents <= 0) nEvents = chain->GetEntries();
 
@@ -78,6 +98,7 @@ void MakeFlow(TString data      = "",
   // --- Add the tasks ---------------------------------------------
   gROOT->LoadMacro("AddTaskForwardFlow.C");
   AddTaskForwardFlow(type, mc, addFlow, addFType, addFOrder);
+  mgr->SetDebugLevel(1);
 
   // --- Run the analysis --------------------------------------------
   TStopwatch t;
index 2f921427972502048668efd465d09f018b599034..23990b3c1d45905a07a126d16a8ae596a1db0923 100644 (file)
@@ -16,8 +16,8 @@
  * t.Run("LOCAL", "FULL", -1, false, false);
  * @endcode 
  *
- * @ingroup pwg2_forward_flow
- * @ingroup pwg2_forward_trains
+ * @ingroup pwglf_forward_flow
+ * @ingroup pwglf_forward_trains
  */
 class MakeFlowTrain : public TrainSetup
 {
@@ -65,7 +65,7 @@ public:
   void Run(const char* mode, const char* oper, 
           Int_t nEvents=-1, Bool_t usePar=false)
   {
-    Exec("AOD", mode, oper, nEvents, false, usePar);
+    Exec("AOD", mode, oper, nEvents, false, usePar, fDebug);
   }
   /** 
    * Run this analysis 
@@ -78,8 +78,9 @@ public:
   void Run(EMode mode, EOper oper, Int_t nEvents=-1, 
           Bool_t usePar=false)
   {
-    Exec(kAOD, mode, oper, nEvents, false, usePar);
+    Exec(kAOD, mode, oper, nEvents, false, usePar, fDebug);
   }
+  void SetDebugLevel(Int_t dbg = 0) { fDebug = dbg; }
 protected:
   /** 
    * Create the tasks 
@@ -118,6 +119,7 @@ protected:
   char* fAddFlow;
   Int_t fAddFType;
   Int_t fAddFOrder;
+  Int_t fDebug;
 };
 //
 // EOF