Mega commit.
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 Apr 2011 13:55:47 +0000 (13:55 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 Apr 2011 13:55:47 +0000 (13:55 +0000)
In general, more diagnostics histograms.

Event inspector can now flag #_SPD_Cluster > 0 events

New implementation of sharing algorithm that zero's candidates
if not used.

dN/deta tasks now outputs summary stacks.

38 files changed:
PWG2/FORWARD/analysis2/AliAODForwardMult.cxx
PWG2/FORWARD/analysis2/AliAODForwardMult.h
PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx
PWG2/FORWARD/analysis2/AliBasedNdetaTask.h
PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx
PWG2/FORWARD/analysis2/AliFMDCorrELossFit.h
PWG2/FORWARD/analysis2/AliFMDCorrector.cxx
PWG2/FORWARD/analysis2/AliFMDCorrector.h
PWG2/FORWARD/analysis2/AliFMDDensityCalculator.cxx
PWG2/FORWARD/analysis2/AliFMDDensityCalculator.h
PWG2/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx
PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx
PWG2/FORWARD/analysis2/AliFMDEventInspector.h
PWG2/FORWARD/analysis2/AliFMDHistCollector.cxx
PWG2/FORWARD/analysis2/AliFMDHistCollector.h
PWG2/FORWARD/analysis2/AliFMDMCCorrector.cxx
PWG2/FORWARD/analysis2/AliFMDMCCorrector.h
PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx
PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx
PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.h
PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx
PWG2/FORWARD/analysis2/AliFMDSharingFilter.h
PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx
PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.h
PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx
PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.h
PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.cxx
PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.h
PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.cxx
PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.h
PWG2/FORWARD/analysis2/AliForwardUtil.cxx
PWG2/FORWARD/analysis2/AliForwardUtil.h
PWG2/FORWARD/analysis2/AliForwarddNdetaTask.cxx
PWG2/FORWARD/analysis2/AliForwarddNdetaTask.h
PWG2/FORWARD/analysis2/DrawdNdeta.C
PWG2/FORWARD/analysis2/ForwardAODConfig.C
PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/OtherData.C

index acb8908dba5d55340ac75c69cb3ac85e36d172ec..e6f6029ab6fc0f609b7abbc02027471b124a706c 100644 (file)
@@ -35,7 +35,8 @@ AliAODForwardMult::AliAODForwardMult()
     fHist(),
     fTriggers(0),
     fIpZ(fgkInvalidIpZ), 
-    fCentrality(-1)
+    fCentrality(-1),                           
+    fNClusters(0)
 {
   // 
   // Constructor 
@@ -49,7 +50,8 @@ AliAODForwardMult::AliAODForwardMult(Bool_t isMC)
          200, -4, 6, 20, 0, 2*TMath::Pi()),
     fTriggers(0),
     fIpZ(fgkInvalidIpZ), 
-    fCentrality(-1)
+    fCentrality(-1),                           
+    fNClusters(0)
 {
   // 
   // Constructor 
@@ -87,8 +89,9 @@ AliAODForwardMult::Clear(Option_t* option)
   //  option   Passed to TH1::Reset 
   // 
   fHist.Reset(option);
-  fTriggers = 0;
-  fIpZ      = fgkInvalidIpZ;
+  fTriggers  = 0;
+  fIpZ       = fgkInvalidIpZ;
+  fNClusters = 0;
 }
 //____________________________________________________________________
 void
@@ -160,13 +163,16 @@ AliAODForwardMult::Browse(TBrowser* b)
   static TObjString ipz;
   static TObjString trg;
   static TObjString cnt;
+  static TObjString ncl;
   ipz = Form("ip_z=%fcm", fIpZ);
   trg = GetTriggerString(fTriggers);
   cnt = Form("%+6.1f%%", fCentrality);
+  ncl = Form("%d clusters", fNClusters);
   b->Add(&fHist);
   b->Add(&ipz);
   b->Add(&trg);
   b->Add(&cnt);
+  b->Add(&ncl);
 }
 
 //____________________________________________________________________
@@ -189,6 +195,7 @@ AliAODForwardMult::GetTriggerString(UInt_t mask)
   if ((mask & kC)           != 0x0) trg.Append("C ");
   if ((mask & kE)           != 0x0) trg.Append("E ");
   if ((mask & kMCNSD)       != 0x0) trg.Append("MCNSD ");
+  if ((mask & kNClusterGt0) != 0x0) trg.Append("NCluster>0 ");
   return trg.Data();
 }
   
@@ -224,6 +231,7 @@ AliAODForwardMult::MakeTriggerHistogram(const char* name)
   ret->GetXaxis()->SetBinLabel(kBinMCNSD,       "NSD (MC truth)");
   ret->GetXaxis()->SetBinLabel(kBinPileUp,      "w/Pileup");
   ret->GetXaxis()->SetBinLabel(kBinOffline,     "w/Offline");
+  ret->GetXaxis()->SetBinLabel(kBinNClusterGt0, "w/N_{cluster}>1");
   ret->GetXaxis()->SetBinLabel(kWithVertex,     "w/Vertex");
   ret->GetXaxis()->SetBinLabel(kWithTrigger,    "w/Selected trigger");
   ret->GetXaxis()->SetBinLabel(kAccepted,       "Accepted by cut");
@@ -251,6 +259,8 @@ AliAODForwardMult::MakeTriggerMask(const char* what)
     else if (s.CompareTo("A")     == 0) trgMask |= AliAODForwardMult::kA;
     else if (s.CompareTo("C")     == 0) trgMask |= AliAODForwardMult::kC;
     else if (s.CompareTo("E")     == 0) trgMask |= AliAODForwardMult::kE;
+    else if (s.CompareTo("NCluster>0") == 0) 
+      trgMask |= AliAODForwardMult::kNClusterGt0;
     else 
       AliWarningGeneral("MakeTriggerMask", 
                        Form("Unknown trigger %s", s.Data()));
@@ -303,6 +313,7 @@ AliAODForwardMult::CheckEvent(Int_t    triggerMask,
     if (IsTriggerBits(kPileUp))         hist->AddBinContent(kBinPileUp);
     if (IsTriggerBits(kMCNSD))          hist->AddBinContent(kBinMCNSD);
     if (IsTriggerBits(kOffline))        hist->AddBinContent(kBinOffline);
+    if (IsTriggerBits(kNClusterGt0))    hist->AddBinContent(kBinNClusterGt0);
   }
   // Check if we have an event of interest. 
   if (!IsTriggerBits(triggerMask|kB)) return false;
index eb6c96eb0c92d51169517312a7684439c61db2eb..19ab1687cd6d3e8ce4ff9fc311b06bec568ee30a 100644 (file)
@@ -127,7 +127,9 @@ public:
     /** true NSD from MC */
     kMCNSD    = 0x400,    
     /** Offline MB triggered */
-    kOffline  = 0x800
+    kOffline  = 0x800,
+    /** At least one SPD cluster */ 
+    kNClusterGt0 = 0x1000
   };
   /** 
    * Bin numbers in trigger histograms 
@@ -144,6 +146,7 @@ public:
     kBinPileUp, 
     kBinMCNSD,
     kBinOffline,
+    kBinNClusterGt0,
     kWithTrigger, 
     kWithVertex, 
     kAccepted
@@ -201,13 +204,27 @@ public:
    */
   void SetTriggerBits(UInt_t bits) { fTriggers |= bits; } // Set trigger bits
   /** 
-   * Check if bit(s) are set in the trigger mask 
+   * Check if all bit(s) are set in the trigger mask.  Note, this is
+   * an @e and between the bits.  If you need an @e or you should use
+   * the member function IsTriggerOrBits
    * 
    * @param bits Bits to test for 
    * 
-   * @return 
+   * @return true if all enabled bits in the argument is also set in
+   * the trigger word
    */
   Bool_t IsTriggerBits(UInt_t bits) const;
+  /** 
+   * Check if any of bit(s) are enabled in the trigger word.  This is
+   * an @e or between the selected bits.  If you need and @a and you
+   * should use the member function IsTriggerBits;
+   * 
+   * @param bits Bits to check for 
+   * 
+   * @return true if any of the enabled bits in the arguments are also
+   * enabled in the trigger mask
+   */
+  Bool_t IsTriggerOrBits(UInt_t bits) const;
   /** 
    * Whether we have any trigger bits 
    */
@@ -316,7 +333,18 @@ public:
    * @return 
    */
   Bool_t  HasCentrality() const { return !(fCentrality  < 0); }
-  
+  /** 
+   * Get the number of SPD clusters seen in @f$ |\eta|<1@f$ 
+   * 
+   * @return Number of SPD clusters seen
+   */
+  UShort_t GetNClusters() const { return fNClusters; }
+  /** 
+   * Set the number of SPD clusters seen in @f$ |\eta|<1@f$ 
+   * 
+   * @param n Number of SPD clusters 
+   */
+  void SetNClusters(UShort_t n) { fNClusters = n; }
   /** 
    * Get the name of the object 
    * 
@@ -385,14 +413,15 @@ public:
    */
   static UInt_t MakeTriggerMask(const char* what);
 protected: 
-  Bool_t  fIsMC;     // Whether this is from MC 
-  TH2D    fHist;     // Histogram of d^2N_{ch}/(deta dphi) for this event
-  UInt_t  fTriggers; // Trigger bit mask 
-  Float_t fIpZ;      // Z coordinate of the interaction point
-  Float_t fCentrality; // Event centrality 
+  Bool_t   fIsMC;       // Whether this is from MC 
+  TH2D     fHist;       // Histogram of d^2N_{ch}/(deta dphi) for this event
+  UInt_t   fTriggers;   // Trigger bit mask 
+  Float_t  fIpZ;        // Z coordinate of the interaction point
+  Float_t  fCentrality; // Event centrality 
+  UShort_t fNClusters;  // Number of SPD clusters in |eta|<1
 
   static const Float_t fgkInvalidIpZ; // Invalid IpZ value 
-  ClassDef(AliAODForwardMult,2); // AOD forward multiplicity 
+  ClassDef(AliAODForwardMult,3); // AOD forward multiplicity 
 };
 
 //____________________________________________________________________
@@ -408,6 +437,12 @@ AliAODForwardMult::IsTriggerBits(UInt_t bits) const
 { 
   return HasTrigger() && ((fTriggers & bits) == bits); 
 }
+//____________________________________________________________________
+inline Bool_t 
+AliAODForwardMult::IsTriggerOrBits(UInt_t bits) const 
+{ 
+  return HasTrigger() && ((fTriggers & bits) != 0);
+}
 
 #endif
 // Local Variables:
index bb7cab406fbbfab315e7d21a43b7ed6e3adf2e6a..8595c84f103476ddb5e3ae7dfabda5cc7b6bc163 100644 (file)
@@ -13,6 +13,7 @@
 #include "AliForwardUtil.h"
 #include "AliAODForwardMult.h"
 #include <TFile.h>
+#include <TStyle.h>
 
 //____________________________________________________________________
 AliBasedNdetaTask::AliBasedNdetaTask()
@@ -153,7 +154,6 @@ AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
   //    high High cut
   //
   CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
-  AliInfo(Form("Adding centrality bin %p [%3d,%3d] @ %d", bin, low, high, at));
   fListOfCentralities->AddAtAndExpand(bin, at);
 }
 
@@ -548,10 +548,67 @@ AliBasedNdetaTask::Terminate(Option_t *)
   // Loop over centrality bins 
   TIter next(fListOfCentralities);
   CentralityBin* bin = 0;
-  while ((bin = static_cast<CentralityBin*>(next()))) 
+  gStyle->SetPalette(1);
+  THStack* dndetaStack        = new THStack("dndeta", "dN/d#eta");
+  THStack* dndetaStackRebin   = new THStack(Form("dndeta_rebin%02d", fRebin), 
+                                           "dN_{ch}/d#eta");
+  THStack* dndetaMCStack      = new THStack("dndetaMC", "dN_{ch}/d#eta");
+  THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin), 
+                                           "dN_{ch}/d#eta");
+  
+  while ((bin = static_cast<CentralityBin*>(next()))) {
     bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, fTriggerEff,
             fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges, 
-            fTriggerMask, GetColor(), GetMarker());
+            fTriggerMask, GetMarker());
+    if (fCentAxis && bin->IsAllBin()) continue;
+    TH1* dndeta      = bin->GetResult(0, false, "");
+    TH1* dndetaSym   = bin->GetResult(0, true,  "");
+    TH1* dndetaMC    = bin->GetResult(0, false, "MC");
+    TH1* dndetaMCSym = bin->GetResult(0, true,  "MC");
+    if (dndeta)      dndetaStack->Add(dndeta);
+    if (dndetaSym)   dndetaStack->Add(dndetaSym);
+    if (dndetaMC)    dndetaMCStack->Add(dndetaMC);
+    if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym);
+    if (fRebin > 1) { 
+      dndeta      = bin->GetResult(fRebin, false, "");
+      dndetaSym   = bin->GetResult(fRebin, true,  "");
+      dndetaMC    = bin->GetResult(fRebin, false, "MC");
+      dndetaMCSym = bin->GetResult(fRebin, true,  "MC");
+      if (dndeta)      dndetaStackRebin->Add(dndeta);
+      if (dndetaSym)   dndetaStackRebin->Add(dndetaSym);
+      if (dndetaMC)    dndetaMCStackRebin->Add(dndetaMC);
+      if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym);
+    }
+  }
+  // Output the stack
+  fOutput->Add(dndetaStack);
+
+  // If available output rebinned stack 
+  if (!dndetaStackRebin->GetHists() || 
+      dndetaStackRebin->GetHists()->GetEntries() <= 0) {
+    AliWarning("No rebinned histograms found");
+    delete dndetaStackRebin;
+    dndetaStackRebin = 0;
+  }
+  if (dndetaStackRebin) fOutput->Add(dndetaStackRebin);
+
+  // If available, output track-ref stack
+  if (!dndetaMCStack->GetHists() || 
+      dndetaMCStack->GetHists()->GetEntries() <= 0) {
+    AliWarning("No MC histograms found");
+    delete dndetaMCStack;
+    dndetaMCStack = 0;
+  }
+  if (dndetaMCStack) fOutput->Add(dndetaMCStack);
+
+  // If available, output rebinned track-ref stack
+  if (!dndetaMCStackRebin->GetHists() || 
+      dndetaMCStackRebin->GetHists()->GetEntries() <= 0) {
+    AliWarning("No rebinned MC histograms found");
+    delete dndetaMCStackRebin;
+    dndetaMCStackRebin = 0;
+  }
+  if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin);
 
   // Output collision energy string 
   if (fSNNString) fOutput->Add(fSNNString->Clone());
@@ -901,6 +958,17 @@ AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
   return *this;
 }
 //____________________________________________________________________
+Int_t
+AliBasedNdetaTask::CentralityBin::GetColor() const 
+{
+  if (IsAllBin()) return kRed+2;
+  Float_t  fc       = (fLow+double(fHigh-fLow)/2) / 100;
+  Int_t    nCol     = gStyle->GetNumberOfColors();
+  Int_t    icol     = TMath::Min(nCol-1,int(fc * nCol + .5));
+  Int_t    col      = gStyle->GetColorPalette(icol);
+  return col;
+}
+//____________________________________________________________________
 const char* 
 AliBasedNdetaTask::CentralityBin::GetListName() const
 {
@@ -1109,6 +1177,38 @@ AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t,
   return scaler;
 }
 
+//________________________________________________________________________
+const char* 
+AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin,
+                                               Bool_t sym, 
+                                               const char* postfix) const
+{
+  static TString n;
+  n = Form("dndeta%s%s",GetName(), postfix);
+  if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
+  if (sym)       n.Append("_mirror");
+  return n.Data();
+}
+//________________________________________________________________________
+TH1* 
+AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin,
+                                           Bool_t sym, 
+                                           const char* postfix) const
+{
+  if (!fOutput) { 
+    AliWarning(Form("No output list defined in %s [%3d,%3d]", GetName(), 
+                   fLow, fHigh));
+    return 0;
+  }
+  TString  n = GetResultName(rebin, sym, postfix);
+  TObject* o = fOutput->FindObject(n.Data());
+  if (!o) { 
+    // AliWarning(Form("Object %s not found in output list", n.Data()));
+    return 0;
+  }
+  return static_cast<TH1*>(o);
+}
+
 //________________________________________________________________________
 void 
 AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,  
@@ -1120,7 +1220,6 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,
                                             bool        symmetrice, 
                                             Int_t       rebin, 
                                             bool        cutEdges, 
-                                            Int_t       color, 
                                             Int_t       marker)
 {
   // 
@@ -1170,9 +1269,13 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,
   copy->Scale(1., "width");
 
   // --- Set some histogram attributes -------------------------------
-  SetHistogramAttributes(dndeta,  color, marker, Form("ALICE %s", GetName()));
-  SetHistogramAttributes(accNorm, color, marker, Form("ALICE %s normalisation", 
-                                                     GetName())); 
+  TString post;
+  if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix);
+  SetHistogramAttributes(dndeta,  GetColor(), marker, 
+                        Form("ALICE %s%s", GetName(), post.Data()));
+  SetHistogramAttributes(accNorm, GetColor(), marker, 
+                        Form("ALICE %s normalisation%s", 
+                             GetName(), post.Data())); 
 
   // --- Make symmetric extensions and rebinnings --------------------
   if (symmetrice)   fOutput->Add(Symmetrice(dndeta));
@@ -1196,7 +1299,6 @@ AliBasedNdetaTask::CentralityBin::End(TList*      sums,
                                      Bool_t      corrEmpty, 
                                      Bool_t      cutEdges, 
                                      Int_t       triggerMask,
-                                     Int_t       color, 
                                      Int_t       marker) 
 {
   // 
@@ -1256,85 +1358,17 @@ AliBasedNdetaTask::CentralityBin::End(TList*      sums,
 
   // --- Make result and store ---------------------------------------
   MakeResult(fSum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
-            scaler, symmetrice, rebin, cutEdges, color, marker);
+            scaler, symmetrice, rebin, cutEdges, marker);
 
   // --- Process result from TrackRefs -------------------------------
-  if (fSumMC) 
+  if (fSumMC) 
     MakeResult(fSumMC, "MC", rootProj, corrEmpty, 
               (scheme & kShape) ? shapeCorr : 0,
-              scaler, symmetrice, rebin, cutEdges, color+2, marker);
-  }
+              scaler, symmetrice, rebin, cutEdges, marker+2);
 
   // Temporary stuff 
-  if (!IsAllBin()) return;
-
-  TFile* forward = TFile::Open("forward.root", "READ");
-  if (!forward)  { 
-    AliWarning(Form("No forward.root file found"));
-    return;
-  }
-
-  TH1D* shapeCorrProj = 0;
-  if (shapeCorr) {
-    shapeCorrProj = static_cast<const TH2D*>(shapeCorr)->ProjectionX();
-    shapeCorrProj->Scale(1. / shapeCorr->GetNbinsY());
-    shapeCorrProj->SetDirectory(0);
-    fOutput->Add(shapeCorrProj);
-  }
-
-  TList* official = static_cast<TList*>(forward->Get("official"));
-  if (official) { 
-    TH1F*  histEta  = static_cast<TH1F*>(official->FindObject("fHistEta"));
-    if (histEta)  {
-      TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(), 
-                           histEta->GetXaxis()->GetXmin(), 
-                           histEta->GetXaxis()->GetXmax());
-      for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
-       oEta->SetBinContent(i, histEta->GetBinContent(i));
-       oEta->SetBinError(i, histEta->GetBinError(i));
-      }
-      if (shapeCorrProj) oEta->Divide(shapeCorrProj);
-      oEta->Scale(1./ntotal, "width");
-      oEta->SetDirectory(0);
-      oEta->SetMarkerStyle(marker+4);
-      oEta->SetMarkerColor(color+5);
-      fOutput->Add(oEta);
-      fOutput->Add(Rebin(oEta, rebin, false));
-    }
-    else 
-      AliWarning(Form("Couldn't find histogram fHistEta in list %s", 
-                     official->GetName()));
-  }
-  else 
-    AliWarning(Form("Couldn't find list 'official' in %s",forward->GetName()));
-
-  TList* tracks = static_cast<TList*>(forward->Get("tracks"));
-  if (tracks) { 
-    TH1F*  histEta  = static_cast<TH1F*>(tracks->FindObject("fHistEta"));
-    if (histEta)  {
-      TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(), 
-                           histEta->GetXaxis()->GetXmin(), 
-                           histEta->GetXaxis()->GetXmax());
-      for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
-       oEta->SetBinContent(i, histEta->GetBinContent(i));
-       oEta->SetBinError(i, histEta->GetBinError(i));
-      }
-      if (shapeCorrProj) oEta->Divide(shapeCorrProj);
-      oEta->Scale(1./ntotal, "width");
-      oEta->SetDirectory(0);
-      oEta->SetMarkerStyle(marker);
-      oEta->SetMarkerColor(color+5);
-      fOutput->Add(oEta);
-      fOutput->Add(Rebin(oEta, rebin, false));
-    }
-    else 
-      AliWarning(Form("Couldn't find histogram fHistEta in list %s", 
-                     tracks->GetName()));
-  }
-  else 
-    AliWarning(Form("Couldn't find list 'tracks' in %s",forward->GetName()));
+  // if (!IsAllBin()) return;
 
-  forward->Close();
 }
 
 //
index 281937bb28c2caa15ffb82a1c87bdca1491b482b..3efcee3a2add666bace920abda65caf339f1c923 100644 (file)
@@ -463,7 +463,6 @@ protected:
      * @param symmetrice Whether to make symmetric extensions 
      * @param rebin      Whether to rebin
      * @param cutEdges   Whether to cut edges when rebinning 
-     * @param color      Color of markers 
      * @param marker     Marker style 
      */
     virtual void MakeResult(const TH2D* sum,  
@@ -475,7 +474,6 @@ protected:
                            bool        symmetrice, 
                            Int_t       rebin, 
                            bool        cutEdges, 
-                           Int_t       color, 
                            Int_t       marker);
     /** 
      * End of processing 
@@ -491,7 +489,6 @@ protected:
      * @param corrEmpty   Whether to correct for empty bins
      * @param cutEdges    Whether to cut edges when rebinning
      * @param triggerMask Trigger mask 
-     * @param color       Base colour for markers 
      * @param marker      Marker style 
      */
     virtual void End(TList*      sums, 
@@ -505,7 +502,6 @@ protected:
                     Bool_t      corrEmpty, 
                     Bool_t      cutEdges, 
                     Int_t       triggerMask,
-                    Int_t       color,
                     Int_t       marker);
     /**
      * @{
@@ -540,6 +536,14 @@ protected:
      */
     TH1I* GetTrigggers() { return fTriggers; }
     /** @} */
+
+    Int_t GetColor() const;
+    TList* GetResults() const { return fOutput; }
+    const char* GetResultName(Int_t rebin, Bool_t sym, 
+                             const char* postfix="") const;
+    TH1* GetResult(Int_t rebin, Bool_t sym, 
+                  const char* postfix="") const;
+
   protected:
     /** 
      * Create sum histogram 
index e3f42187e82e9bc47faea97be94799d620f50c54..3991420b80e6e1e823fe8be84b3006fdbdf9e80e 100644 (file)
@@ -512,6 +512,17 @@ AliFMDCorrELossFit::ELossFit::Draw(Option_t* option)
 //____________________________________________________________________
 #define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
 
+//____________________________________________________________________
+Double_t
+AliFMDCorrELossFit::ELossFit::GetLowerBound(Double_t f, 
+                                           Bool_t includeSigma) const
+{
+  // 
+  // Return 
+  //    Delta - f * (xi + sigma)
+  return fDelta - f * (fXi + (includeSigma ? fSigma : 0));
+}
+
 //____________________________________________________________________
 void 
 AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu, 
@@ -804,24 +815,24 @@ AliFMDCorrELossFit::FindFit(UShort_t  d, Char_t r, Int_t etabin) const
     return 0;
   }
   if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) { 
-    AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", 
-                 etabin, 1, fEtaAxis.GetNbins(), d, r));
+    // AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", 
+    //              etabin, 1, fEtaAxis.GetNbins(), d, r));
     return 0;
   }
   if (etabin > ringArray->GetEntriesFast()) { 
-    AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", 
-                 etabin, 1, ringArray->GetEntriesFast(), d, r));
+    // AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", 
+    //                      etabin, 1, ringArray->GetEntriesFast(), d, r));
     return 0;
   }
   else if (etabin >= ringArray->GetEntriesFast()) { 
-    AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, " 
-                   "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r,
-                   etabin-1));
+    // AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, " 
+    //             "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r,
+    //             etabin-1));
     etabin--;
   }
   else if (!ringArray->At(etabin)) { 
-    AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d", 
-                   etabin, d, r, etabin+1));
+    // AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d", 
+    //                     etabin, d, r, etabin+1));
     etabin++;
   }
   return static_cast<ELossFit*>(ringArray->At(etabin));
@@ -898,6 +909,38 @@ AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r)
   return static_cast<TObjArray*>(fRings.At(idx));
 }
 
+//____________________________________________________________________
+Double_t
+AliFMDCorrELossFit::GetLowerBound(UShort_t  d, Char_t r, Int_t etabin,
+                                 Double_t f, Bool_t showErrors, 
+                                 Bool_t includeSigma) const
+{
+  ELossFit* fit = GetFit(d, r, etabin);
+  if (!fit) { 
+    if (showErrors) {
+      AliWarning(Form("No fit for FMD%d%c @ etabin=%d", d, r, etabin));
+    }
+    return -1024;
+  }
+  return fit->GetLowerBound(f, includeSigma);
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDCorrELossFit::GetLowerBound(UShort_t  d, Char_t r, Double_t eta,
+                                 Double_t f, Bool_t showErrors, 
+                                 Bool_t includeSigma) const
+{
+  Int_t bin = FindEtaBin(eta);
+  if (bin <= 0) { 
+    if (showErrors)
+      AliError(Form("eta=%f out of bounds for FMD%d%c", eta, d, r));
+    return -1024;
+  }
+  return GetLowerBound(d, r, bin, f, showErrors, includeSigma);
+}
+
+//____________________________________________________________________
 namespace { 
   TH1D* MakeHist(const TAxis& axis, const char* name, const char* title, 
                 Int_t color)
@@ -916,6 +959,8 @@ namespace {
   }
 }
 
+  
+
 #define IDX2RING(I) (i == 0 || i == 1 || i == 3 ? 'I' : 'O')
 #define IDX2DET(I)  (i == 0 ? 1 : (i == 1 || i == 2 ? 2 : 3))
 //____________________________________________________________________
index 4ca91e5be8ef7cc27eafca31f169c929e7580443..185ad3ffc3aeb7e199e671c901512f5ff0c662ee 100644 (file)
@@ -331,6 +331,15 @@ public:
      * @return 
      */
     const Char_t* GetName() const;
+    /** 
+     * Calculate the lower bound 
+     * 
+     * @param f             Width factor
+     * @param includeSigma  Whether to include sigma
+     * 
+     * @return @f$ \Delta - f (\xi + \sigma)@f$
+     */
+    Double_t GetLowerBound(Double_t f, Bool_t includeSigma) const;
     /** 
      * Calculate the quality 
      */
@@ -526,6 +535,42 @@ public:
   ELossFit* GetFit(UShort_t d, Char_t r, Int_t etabin) const;
   /* @} */
 
+  /** 
+   * @{ 
+   * @name Lower bounds on fits 
+   */
+  /** 
+   * Get the lower validity bound of the fit
+   * 
+   * @param d            Detector
+   * @param r            Ring
+   * @param etaBin       Eta bin (1-based)
+   * @param f            Factor on xi (and sigma)
+   * @param showErrors   Show errors
+   * @param includeSigma Whether to include sigma 
+   * 
+   * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$ 
+   */
+  Double_t GetLowerBound(UShort_t d, Char_t r, Int_t etaBin, 
+                        Double_t f, Bool_t showErrors=true,
+                        Bool_t includeSigma=true) const;
+  /** 
+   * Get the lower validity bound of the fit
+   * 
+   * @param d            Detector
+   * @param r            Ring
+   * @param eta          Eta value
+   * @param f            Factor on xi (and sigma)
+   * @param showErrors   Show errors
+   * @param includeSigma Whether to include sigma 
+   * 
+   * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$ 
+   */
+  Double_t GetLowerBound(UShort_t d, Char_t r, Double_t eta, 
+                        Double_t f, Bool_t showErrors=true,
+                        Bool_t includeSigma=true) const;
+
+  
   /**                                          
    * @{ 
    * @name Miscellaneous
index d8fb3cc4dbaad8b1411fafec18f2d262e15c58a6..62a6f29f3c055eba89a4594263feff7672aef4d3 100644 (file)
@@ -15,6 +15,7 @@
 #include "AliLog.h"
 #include <TH2D.h>
 #include <TROOT.h>
+#include <THStack.h>
 #include <iostream>
 #include <iomanip>
 
@@ -178,8 +179,8 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
       if (fUseSecondaryMap) {
         TH2D*  bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
         if (!bg) {
-          AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
-                          d, r, uvb));
+          AliWarning(Form("No secondary correction for FMDM%d%c "
+                         "in vertex bin %d", d, r, uvb));
           continue;
         }
         // Divide by primary/total ratio
@@ -198,18 +199,14 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
       if (fUseAcceptance) {
         TH2D*  ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
         if (!ac) {
-          AliWarning(Form("No acceptance correction for FMD%d%c in vertex bin %d",
-                       d, r, uvb));
+          AliWarning(Form("No acceptance correction for FMD%d%c in "
+                         "vertex bin %d", d, r, uvb));
           continue;
         }
         // Divide by the acceptance correction
         h->Divide(ac);
       }
 
-
-
-      
-
       if (fUseMergingEfficiency) {
        if (!fcm.GetMergingEfficiency()) { 
          AliWarning("No merging efficiencies");
@@ -307,8 +304,19 @@ AliFMDCorrector::ScaleHistograms(const TList* dir, Int_t nEvents)
 
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
-  while ((o = static_cast<RingHistos*>(next())))
+  THStack* sums = new THStack("sums", "Sums of ring results");
+  while ((o = static_cast<RingHistos*>(next()))) {
     o->ScaleHistograms(d, nEvents);
+    TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1, 
+                                        o->fDensity->GetNbinsY(),"e");
+    sum->Scale(1., "width");
+    sum->SetTitle(o->GetName());
+    sum->SetDirectory(0);
+    sum->SetYTitle("#sum N_{ch,primary}");
+    sums->Add(sum);
+  }
+  d->Add(sums);
+
 }
 //____________________________________________________________________
 void
@@ -372,11 +380,13 @@ AliFMDCorrector::RingHistos::RingHistos(UShort_t d, Char_t r)
   //    d detector
   //    r ring 
   //
-  fDensity = new TH2D(Form("FMD%d%c_Primary_Density", d, r)
-                     Form("in FMD%d%c", d, r)
+  fDensity = new TH2D("primaryDensity"
+                     "#sum N_{ch,primary}/(#Delta#eta#Delta#phi)"
                      200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
                      0, 2*TMath::Pi());
   fDensity->SetDirectory(0);
+  fDensity->SetMarkerColor(Color());
+  fDensity->Sumw2();
   fDensity->SetXTitle("#eta");
   fDensity->SetYTitle("#phi [radians]");
   fDensity->SetZTitle("Primary N_{ch} density");
@@ -437,7 +447,7 @@ AliFMDCorrector::RingHistos::Output(TList* dir)
 
 //____________________________________________________________________
 void
-AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t /*nEvents*/)
+AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
 { 
   // 
   // Scale the histograms to the total number of events 
@@ -448,8 +458,8 @@ AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t /*nEvents*/)
   TList* l = GetOutputList(dir);
   if (!l) return; 
 
-  //TH1* density = GetOutputHist(l,Form("%s_Primary_Density", fName.Data()));
-  //if (density) density->Scale(1./nEvents);
+  TH1* density = GetOutputHist(l,"primaryDensity");
+  if (density) density->Scale(1./nEvents);
 }
 
 //____________________________________________________________________
index 4471cbd7103d05bd784817c56866e4e3768ff91c..6665a929ad5e4a0c969a8eeeed8ed4a5cde735e5 100644 (file)
@@ -167,7 +167,7 @@ public:
    * 
    * @param option Not used 
    */
-  void Print(Option_t* option="") const;
+  virtual void Print(Option_t* option="") const;
 protected:
   /** 
    * Internal data structure to keep track of the histograms
index 8320b9c6b71577d9bda5a420d92d64702cb39a9e..d6654c70e7b8e7e53387767a34aaf5ab8d656462 100644 (file)
@@ -12,6 +12,7 @@
 #include "AliLog.h"
 #include <TH2D.h>
 #include <TProfile.h>
+#include <THStack.h>
 #include <TROOT.h>
 #include <iostream>
 #include <iomanip>
@@ -26,11 +27,14 @@ AliFMDDensityCalculator::AliFMDDensityCalculator()
   : TNamed(), 
     fRingHistos(),
     fMultCut(0),
+    fNXi(1),
+    fIncludeSigma(true),
     fSumOfWeights(0),
     fWeightedSum(0),
     fCorrections(0),
     fMaxParticles(5),
     fUsePoisson(false),
+    fUsePhiAcceptance(false),
     fAccI(0),
     fAccO(0),
     fFMD1iMax(0),
@@ -38,6 +42,8 @@ AliFMDDensityCalculator::AliFMDDensityCalculator()
     fFMD2oMax(0),
     fFMD3iMax(0),
     fFMD3oMax(0),
+    fMaxWeights(0),
+    fLowCuts(0),
     fEtaLumping(1), 
     fPhiLumping(1),    
     fDebug(0)
@@ -53,11 +59,14 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
   : TNamed("fmdDensityCalculator", title), 
     fRingHistos(), 
     fMultCut(0),
+    fNXi(1),
+    fIncludeSigma(true),
     fSumOfWeights(0),
     fWeightedSum(0),
     fCorrections(0),
     fMaxParticles(5),
     fUsePoisson(false),
+    fUsePhiAcceptance(false),
     fAccI(0),
     fAccO(0),
     fFMD1iMax(0),
@@ -65,6 +74,8 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
     fFMD2oMax(0),
     fFMD3iMax(0),
     fFMD3oMax(0),
+    fMaxWeights(0),
+    fLowCuts(0),
     fEtaLumping(5), 
     fPhiLumping(1),
     fDebug(0)
@@ -97,6 +108,17 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
 
   fAccI = GenerateAcceptanceCorrection('I');
   fAccO = GenerateAcceptanceCorrection('O');
+
+  fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use", 
+                        1, 0, 1, 1, 0, 1);
+  fMaxWeights->SetXTitle("#eta");
+  fMaxWeights->SetDirectory(0);
+
+  fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1);
+  fLowCuts->SetXTitle("#eta");
+  fLowCuts->SetDirectory(0);
+
+  for (Int_t i = 0; i < 5; i++) fMultCuts[i] = 0;
 }
 
 //____________________________________________________________________
@@ -105,11 +127,14 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const
   : TNamed(o), 
     fRingHistos(), 
     fMultCut(o.fMultCut),
+    fNXi(o.fNXi),
+    fIncludeSigma(o.fIncludeSigma),
     fSumOfWeights(o.fSumOfWeights),
     fWeightedSum(o.fWeightedSum),
     fCorrections(o.fCorrections),
     fMaxParticles(o.fMaxParticles),
     fUsePoisson(o.fUsePoisson),
+    fUsePhiAcceptance(o.fUsePhiAcceptance),
     fAccI(o.fAccI),
     fAccO(o.fAccO),
     fFMD1iMax(o.fFMD1iMax),
@@ -117,6 +142,8 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const
     fFMD2oMax(o.fFMD2oMax),
     fFMD3iMax(o.fFMD3iMax),
     fFMD3oMax(o.fFMD3oMax),
+    fMaxWeights(o.fMaxWeights),
+    fLowCuts(o.fLowCuts),
     fEtaLumping(o.fEtaLumping), 
     fPhiLumping(o.fPhiLumping),
     fDebug(o.fDebug)
@@ -157,9 +184,12 @@ AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
   TNamed::operator=(o);
 
   fMultCut      = o.fMultCut;
+  fNXi          = o.fNXi;
+  fIncludeSigma = o.fIncludeSigma;
   fDebug        = o.fDebug;
   fMaxParticles = o.fMaxParticles;
   fUsePoisson   = o.fUsePoisson;
+  fUsePhiAcceptance  = o.fUsePhiAcceptance;
   fAccI         = o.fAccI;
   fAccO         = o.fAccO;
   fFMD1iMax     = o.fFMD1iMax;
@@ -167,6 +197,8 @@ AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
   fFMD2oMax     = o.fFMD2oMax;
   fFMD3iMax     = o.fFMD3iMax;
   fFMD3oMax     = o.fFMD3oMax;
+  fMaxWeights   = o.fMaxWeights;
+  fLowCuts      = o.fLowCuts;
   fEtaLumping   = o.fEtaLumping;
   fPhiLumping   = o.fPhiLumping;
   fRingHistos.Delete();
@@ -220,10 +252,27 @@ AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
   
   return static_cast<RingHistos*>(fRingHistos.At(idx));
 }
-    
+
+//____________________________________________________________________
+void 
+AliFMDDensityCalculator::SetMultCuts(Double_t fmd1i, 
+                                    Double_t fmd2i, 
+                                    Double_t fmd2o, 
+                                    Double_t fmd3i, 
+                                    Double_t fmd3o) 
+{
+  fMultCuts[0] = fmd1i;
+  fMultCuts[1] = fmd2i;
+  fMultCuts[2] = fmd2o;
+  fMultCuts[3] = fmd3i;
+  fMultCuts[4] = fmd3o;
+  fMultCut = (fmd1i+fmd2i+fmd2o+fmd3i+fmd3o) / 5;
+}
+
 //____________________________________________________________________
 Double_t
-AliFMDDensityCalculator::GetMultCut() const
+AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
+                                   Bool_t errors) const
 {
   // 
   // Get the multiplicity cut.  If the user has set fMultCut (via
@@ -233,11 +282,35 @@ AliFMDDensityCalculator::GetMultCut() const
   // Return:
   //    Lower cut on multiplicity
   //
+  Int_t    idx = (d == 1 ? 0 : 2*(d - 2) + 1 + ((r=='I' || r=='i') ? 0 : 1));
+  if (fMultCuts[idx] > 0) return fMultCuts[idx];
   if (fMultCut > 0) return fMultCut;
 
   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
   AliFMDCorrELossFit* fits = fcm.GetELossFit();
-  return fits->GetLowCut();
+  if (fNXi < 0) return fits->GetLowCut();
+
+  return fits->GetLowerBound(d, r, ieta, fNXi, errors, fIncludeSigma);
+}
+    
+//____________________________________________________________________
+Double_t
+AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
+                                   Bool_t /*errors*/) const
+{
+  // 
+  // Get the multiplicity cut.  If the user has set fMultCut (via
+  // SetMultCut) then that value is used.  If not, then the lower
+  // value of the fit range for the enery loss fits is returned.
+  // 
+  // Return:
+  //    Lower cut on multiplicity
+  //
+  AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
+  AliFMDCorrELossFit* fits = fcm.GetELossFit();
+  Int_t iEta = fits->FindEtaBin(eta);
+  
+  return GetMultCut(d, r, iEta);
 }
   
 //____________________________________________________________________
@@ -276,22 +349,27 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
       
       for (UShort_t s=0; s<ns; s++) { 
        for (UShort_t t=0; t<nt; t++) {
-         Float_t mult = fmd.Multiplicity(d,r,s,t);
-         Float_t phi  = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
-         Float_t eta  = fmd.Eta(d,r,s,t);
+         Float_t  mult = fmd.Multiplicity(d,r,s,t);
+         Float_t  phi  = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
+         Float_t  eta  = fmd.Eta(d,r,s,t);
+         Double_t cut  = 1024;
+         if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
          rh->fTotalStrips->Fill(eta, phi);
          
-         if (mult == AliESDFMD::kInvalidMult || mult > 20) { 
+         if (mult == AliESDFMD::kInvalidMult || mult > 20) {
            rh->fEmptyStrips->Fill(eta,phi);
+           rh->fEvsM->Fill(mult,0);
            continue;
          }
-         
-         Float_t n   = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
+
+         Double_t n   = 0;
+         if (cut > 0 && mult > cut) 
+           n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
          
          rh->fEvsN->Fill(mult,n);
          rh->fEtaVsN->Fill(eta, n);
          
-         Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
+         Double_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
          fCorrections->Fill(c);
          if (c > 0) n /= c;
          rh->fEvsM->Fill(mult,n);
@@ -302,7 +380,7 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
          else                  rh->fEmptyStrips->Fill(eta,phi);
            
          h->Fill(eta,phi,n);
-         rh->fDensity->Fill(eta,phi,n);
+         if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
        } // for t
       } // for s 
 
@@ -311,20 +389,32 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
       for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) { 
        for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { 
          Double_t eLossV   = h->GetBinContent(ieta, iphi);
-         Float_t  eta      = h->GetXaxis()->GetBinCenter(ieta);
-         Float_t  phi      = h->GetYaxis()->GetBinCenter(iphi);
+         Double_t eta      = h->GetXaxis()->GetBinCenter(ieta);
+         Double_t phi      = h->GetYaxis()->GetBinCenter(iphi);
          Int_t    jeta     = rh->fEmptyStrips->GetXaxis()->FindBin(eta);
          Int_t    jphi     = rh->fEmptyStrips->GetYaxis()->FindBin(phi);
          Double_t empty    = rh->fEmptyStrips->GetBinContent(jeta, jphi);
          Double_t total    = rh->fTotalStrips->GetBinContent(jeta, jphi);
          Double_t hits     = rh->fBasicHits->GetBinContent(ieta,iphi);
+
          // Mean in region of interest 
          Double_t poissonM = (total <= 0 || empty <= 0 ? 0 : 
                               -TMath::Log(empty / total));
+
+         // Note, that given filled=total-empty, and 
+         //
+         //     m = -log(empty/total)
+         //       = -log(1 - filled/total)
+         // 
+         //     v = m / (1 - exp(-m))
+         //       = -total/filled * (log(total-filled)-log(total))
+         //       = -total / (total-empty) * log(empty/total)
+         //       = total (log(total)-log(empty)) / (total-empty)
+         //  
          Double_t poissonV = hits;
          if(poissonM > 0)
            // Correct for counting statistics and weight by counts 
-           poissonV = (hits * poissonM) / (1 - TMath::Exp(-1*poissonM));
+           poissonV *= poissonM / (1 - TMath::Exp(-1*poissonM));
          Double_t poissonE = TMath::Sqrt(hits);
          if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
          
@@ -333,6 +423,7 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
          if (fUsePoisson) {
            h->SetBinContent(ieta,iphi,poissonV);
            h->SetBinError(ieta,iphi,poissonE);
+           rh->fDensity->Fill(eta, phi, poissonV);
          }
        }
       }
@@ -383,12 +474,39 @@ AliFMDDensityCalculator::CacheMaxWeights()
   fFMD3iMax.Set(nEta);
   fFMD3oMax.Set(nEta);
   
+  fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
+  fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
+  fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
+  fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
+  fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
+  fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
+
+  AliInfo(Form("Get eta axis with %d bins from %f to %f",
+              nEta, eta.GetXmin(), eta.GetXmax()));
+  fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
+  fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
+  fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
+  fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
+  fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
+  fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
+  
   for (Int_t i = 0; i < nEta; i++) {
-    fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
-    fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
-    fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
-    fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
-    fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
+    Double_t w[5];
+    w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
+    w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
+    w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
+    w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
+    w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
+    Double_t l[5];
+    l[0] = GetMultCut(1, 'I', i+1, false);
+    l[1] = GetMultCut(2, 'I', i+1, false);
+    l[2] = GetMultCut(2, 'O', i+1, false);
+    l[3] = GetMultCut(3, 'I', i+1, false);
+    l[4] = GetMultCut(3, 'O', i+1, false);
+    for (Int_t j = 0; j < 5; j++) { 
+      if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
+      if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
+    }
   }
 }
 
@@ -481,7 +599,7 @@ AliFMDDensityCalculator::NParticles(Float_t  mult,
   // Return:
   //    The number of particles 
   //
-  if (mult <= GetMultCut()) return 0;
+  // if (mult <= GetMultCut()) return 0;
   if (lowFlux) return 1;
   
   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
@@ -540,7 +658,8 @@ AliFMDDensityCalculator::Correction(UShort_t d,
   //
   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 
-  Float_t correction = AcceptanceCorrection(r,t);
+  Float_t correction = 1; 
+  if (fUsePhiAcceptance) correction = AcceptanceCorrection(r,t);
   if (lowFlux) { 
     TH1D* dblHitCor = 0;
     if (fcm.GetDoubleHit()) 
@@ -670,8 +789,18 @@ AliFMDDensityCalculator::ScaleHistograms(const TList* dir, Int_t nEvents)
 
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
-  while ((o = static_cast<RingHistos*>(next())))
+  THStack* sums = new THStack("sums", "sums of ring signals");
+  while ((o = static_cast<RingHistos*>(next()))) {
     o->ScaleHistograms(d, nEvents);
+    TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1, 
+                                        o->fDensity->GetNbinsY(),"e");
+    sum->Scale(1., "width");
+    sum->SetTitle(o->GetName());
+    sum->SetDirectory(0);
+    sum->SetYTitle("#sum N_{ch,incl}");
+    sums->Add(sum);
+  }
+  d->Add(sums);
 }
 
 //____________________________________________________________________
@@ -685,6 +814,7 @@ AliFMDDensityCalculator::DefineOutput(TList* dir)
   //    dir List to write in
   //  
   TList* d = new TList;
+  d->SetOwner();
   d->SetName(GetName());
   dir->Add(d);
   d->Add(fWeightedSum);
@@ -692,6 +822,8 @@ AliFMDDensityCalculator::DefineOutput(TList* dir)
   d->Add(fCorrections);
   d->Add(fAccI);
   d->Add(fAccO);
+  d->Add(fMaxWeights);
+  d->Add(fLowCuts);
 
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
@@ -713,10 +845,19 @@ AliFMDDensityCalculator::Print(Option_t* option) const
   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
   ind[gROOT->GetDirLevel()] = '\0';
   std::cout << ind << ClassName() << ": " << GetName() << '\n'
+           << std::boolalpha 
            << ind << " Multiplicity cut:       " << fMultCut << '\n'
+           << ind << " # of (xi+sigma) factor: " << fNXi << '\n'
+           << ind << " Include sigma in cut:   " << fIncludeSigma << '\n'
+           << ind << " Low cut method:         " 
+           << (fMultCut > 0 ? "fixed" : 
+               (fNXi >= 0 ? "xi+sigma" : "fit range")) << '\n'
            << ind << " Max(particles):         " << fMaxParticles << '\n'
+           << ind << " Poisson method:         " << fUsePoisson << '\n'
+           << ind << " Use phi acceptance:     " << fUsePhiAcceptance << '\n'
            << ind << " Eta lumping:            " << fEtaLumping << '\n'
            << ind << " Phi lumping:            " << fPhiLumping << '\n'
+           << std::noboolalpha
            << std::flush;
   TString opt(option);
   opt.ToLower();
@@ -791,70 +932,61 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
   //    d detector
   //    r ring 
   //
-  fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()), 
-                  Form("#Delta E/#Delta E_{mip} vs uncorrected inclusive "
-                       "N_{ch} in %s", fName.Data()), 
-                  2500, -.5, 24.5, 2500, -.5, 24.5);
-  fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()), 
-                  Form("#Delta E/#Delta E_{mip} vs corrected inclusive "
-                       "N_{ch} in %s", fName.Data()), 
-                  2500, -.5, 24.5, 2500, -.5, 24.5);
+  fEvsN = new TH2D("elossVsNnocorr", 
+                  "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
+                  250, -.5, 24.5, 250, -.5, 24.5);
   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
   fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
   fEvsN->Sumw2();
   fEvsN->SetDirectory(0);
-  fEvsM->SetXTitle("#Delta E/#Delta E_{mip}");
-  fEvsM->SetYTitle("Inclusive N_{ch} (corrected)");
-  fEvsM->Sumw2();
+
+  fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
+  fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
   fEvsM->SetDirectory(0);
 
-  fEtaVsN = new TProfile(Form("%s_Eta_N_nocorr", fName.Data()),
-                        Form("Average inclusive N_{ch} vs #eta (uncorrected) "
-                             "in %s", fName.Data()), 200, -4, 6);
-  fEtaVsM = new TProfile(Form("%s_Eta_N_corr", fName.Data()),
-                        Form("Average inclusive N_{ch} vs #eta (corrected) "
-                             "in %s", fName.Data()), 200, -4, 6);
+  fEtaVsN = new TProfile("etaVsNnocorr",
+                        "Average inclusive N_{ch} vs #eta (uncorrected)",
+                        200, -4, 6);
   fEtaVsN->SetXTitle("#eta");
   fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
   fEtaVsN->SetDirectory(0);
   fEtaVsN->SetLineColor(Color());
   fEtaVsN->SetFillColor(Color());
-  fEtaVsM->SetXTitle("#eta");
+
+  fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
+  fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
   fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
   fEtaVsM->SetDirectory(0);
-  fEtaVsM->SetLineColor(Color());
-  fEtaVsM->SetFillColor(Color());
 
 
-  fCorr = new TProfile(Form("%s_corr", fName.Data()),
-                        Form("Average correction in %s", fName.Data()), 
-                      200, -4, 6);
+  fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
   fCorr->SetXTitle("#eta");
   fCorr->SetYTitle("#LT correction#GT");
   fCorr->SetDirectory(0);
   fCorr->SetLineColor(Color());
   fCorr->SetFillColor(Color());
 
-  fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()), 
-                     Form("Inclusive N_{ch} density in %s", fName.Data()), 
+  fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
                      200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
                      0, 2*TMath::Pi());
   fDensity->SetDirectory(0);
+  fDensity->Sumw2();
+  fDensity->SetMarkerColor(Color());
   fDensity->SetXTitle("#eta");
   fDensity->SetYTitle("#phi [radians]");
   fDensity->SetZTitle("Inclusive N_{ch} density");
 
-  fELossVsPoisson = new TH2D(Form("%s_eloss_vs_poisson", fName.Data()),
-                            Form("N_{ch} from energy loss vs from Poission %s",
-                                 fName.Data()), 100, 0, 20, 100, 0, 20);
+  fELossVsPoisson = new TH2D("elossVsPoisson", 
+                            "N_{ch} from energy loss vs from Poission",
+                            100, 0, 20, 100, 0, 20);
   fELossVsPoisson->SetDirectory(0);
   fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
   fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
   fELossVsPoisson->SetZTitle("Correlation");
 
-  fEmptyVsTotal = new TH2D(Form("%s_empty_vs_total", fName.Data())
-                          Form("# of empty strips vs. total @ # strips in %s", 
-                               fName.Data()), 21, -.5, 20.5, 21, -0.5, 20.5);
+  fEmptyVsTotal = new TH2D("emptyVsTotal"
+                          "# of empty strips vs. total # strips", 
+                          21, -.5, 20.5, 21, -0.5, 20.5);
   fEmptyVsTotal->SetDirectory(0);
   fEmptyVsTotal->SetXTitle("Total # strips");
   fEmptyVsTotal->SetYTitle("Empty # strips");
@@ -1030,7 +1162,7 @@ AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
   TList* l = GetOutputList(dir);
   if (!l) return; 
 
-  TH1* density = GetOutputHist(l,Form("%s_Incl_Density", fName.Data()));
+  TH1* density = GetOutputHist(l,"inclDensity");
   if (density) density->Scale(1./nEvents);
 }
 
index fc5391cc1510ef77d4b56b17f450ce442a07d36f..80309cdb18a2a2816c42dc765c309f532100807b 100644 (file)
@@ -132,13 +132,33 @@ public:
    * number of particles that has hit within a region.
    */
   void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
+  /** 
+   * Set whether to use the phi acceptance correction. 
+   * 
+   * @param u If true, use the phi acceptance (default is false)
+   */
+  void SetUsePhiAcceptance(Bool_t u) { fUsePhiAcceptance = u; }
   /** 
    * Set the lower multiplicity cut.  This overrides the setting in
    * the energy loss fits.
    * 
    * @param cut Cut to use 
    */
-  void SetMultCut(Double_t cut) { fMultCut = cut; }
+  void SetMultCut(Double_t cut) { SetMultCuts(cut,cut,cut,cut,cut); }
+  /** 
+   * Set the lower multiplicity cuts 
+   * 
+   * @param fmd1i Lower mulitplicyt cut for FMD1i
+   * @param fmd2i Lower mulitplicyt cut for FMD2i 
+   * @param fmd2o Lower mulitplicyt cut for FMD2o 
+   * @param fmd3i Lower mulitplicyt cut for FMD3i 
+   * @param fmd3o Lower mulitplicyt cut for FMD3o 
+   */
+  void SetMultCuts(Double_t fmd1i, 
+                  Double_t fmd2i, 
+                  Double_t fmd2o, 
+                  Double_t fmd3i, 
+                  Double_t fmd3o); 
   /** 
    * Set the luming factors used in the Poisson method
    * 
@@ -149,6 +169,29 @@ public:
     fEtaLumping = (eta < 1 ? 1 : eta); 
     fPhiLumping = (phi < 1 ? 1 : phi); 
   }
+  /** 
+   * Set the number of landau width to subtract from the most probably
+   * value to get the low cut.
+   * 
+   * @param n Number of @f$ \xi@f$ 
+   */
+  void SetNXi(Double_t nXi) { fNXi = nXi; } 
+  /** 
+   * Whether to include sigma in the number subtracted from the MPV to
+   * get the low cut
+   * 
+   * @param u If true, then low cut is @f$ \Delta_{mp} - n(\xi+\sigma)@f$ 
+   */
+  void SetIncludeSigma(Bool_t u) { fIncludeSigma = u; }
+  /** 
+   * Get the multiplicity cut.  If the user has set fMultCut (via
+   * SetMultCut) then that value is used.  If not, then the lower
+   * value of the fit range for the enery loss fits is returned.
+   * 
+   * @return Lower cut on multiplicity
+   */
+  Double_t GetMultCut(UShort_t d, Char_t r, Double_t eta, 
+                     Bool_t errors=true) const;
   /** 
    * Get the multiplicity cut.  If the user has set fMultCut (via
    * SetMultCut) then that value is used.  If not, then the lower
@@ -156,7 +199,8 @@ public:
    * 
    * @return Lower cut on multiplicity
    */
-  Double_t GetMultCut() const;
+  Double_t GetMultCut(UShort_t d, Char_t r, Int_t ieta, 
+                     Bool_t errors=true) const;
   /** 
    * Print information 
    * 
@@ -341,11 +385,14 @@ protected:
 
   TList    fRingHistos;    //  List of histogram containers
   Double_t fMultCut;       //  Low cut on scaled energy loss
+  Double_t fNXi;           //  Delta-n(xi+sigma) factor 
+  Bool_t   fIncludeSigma;  //  Whether or not to include sigma in cut
   TH1D*    fSumOfWeights;  //  Histogram
   TH1D*    fWeightedSum;   //  Histogram
   TH1D*    fCorrections;   //  Histogram
   UShort_t fMaxParticles;  //  Maximum particle weight to use 
   Bool_t   fUsePoisson;    //  If true, then use poisson statistics 
+  Bool_t   fUsePhiAcceptance; // Whether to correct for corners 
   TH1D*    fAccI;          //  Acceptance correction for inner rings
   TH1D*    fAccO;          //  Acceptance correction for outer rings
   TArrayI  fFMD1iMax;      //  Array of max weights 
@@ -353,12 +400,14 @@ protected:
   TArrayI  fFMD2oMax;      //  Array of max weights 
   TArrayI  fFMD3iMax;      //  Array of max weights 
   TArrayI  fFMD3oMax;      //  Array of max weights 
+  TH2D*    fMaxWeights;    //  Histogram of max weights
+  TH2D*    fLowCuts;       //  Histogram of low cuts
   Int_t    fEtaLumping;    //  How to lump eta bins for Poisson 
   Int_t    fPhiLumping;    //  How to lump phi bins for Poisson 
   Int_t    fDebug;         //  Debug level 
+  Double_t fMultCuts[5];   //  Per-ring multiplicity cuts
 
-
-  ClassDef(AliFMDDensityCalculator,2); // Calculate Nch density 
+  ClassDef(AliFMDDensityCalculator,3); // Calculate Nch density 
 };
 
 #endif
index fb111a97694cc1d4baad0ab9dc1d789a7b3f7942..02f19642263af167639f2d310099eef9074e6100 100644 (file)
@@ -187,7 +187,8 @@ AliFMDEnergyFitterTask::UserExec(Option_t*)
   AliMCEvent* mcevent = MCEvent();
   if(mcevent) {
     AliHeader* header            = mcevent->Header();
-    AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
+    AliGenHijingEventHeader* hijingHeader = 
+      dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
     if(hijingHeader) {
       Float_t b = hijingHeader->ImpactParameter();
       if(b<fbLow || b>fbHigh) return;
@@ -229,13 +230,14 @@ AliFMDEnergyFitterTask::UserExec(Option_t*)
 
     InitializeSubs();
   }
-  Bool_t   lowFlux  = kFALSE;
-  UInt_t   triggers = 0;
-  UShort_t ivz      = 0;
-  Double_t vz       = 0;
-  Double_t cent     = 0;
-  UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, 
-                                             ivz, vz, cent);
+  Bool_t   lowFlux   = kFALSE;
+  UInt_t   triggers  = 0;
+  UShort_t ivz       = 0;
+  Double_t vz        = 0;
+  Double_t cent      = 0;
+  UShort_t nClusters = 0;
+  UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
+                                              ivz, vz, cent, nClusters);
   if (found & AliFMDEventInspector::kNoEvent)    return;
   if (found & AliFMDEventInspector::kNoTriggers) return;
   if (found & AliFMDEventInspector::kNoSPD)     return;
@@ -277,20 +279,6 @@ AliFMDEnergyFitterTask::Terminate(Option_t*)
     return;
   }
   
-#if 0
-  // Get our histograms from the container 
-  TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
-  TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
-  TH1I* hTriggers    = 0;
-  if (!fEventInspector.FetchHistograms(list, hEventsTr, 
-                                      hEventsTrVtx, hTriggers)) { 
-    AliError(Form("Didn't get histograms from event selector "
-                 "(hEventsTr=%p,hEventsTrVtx=%p)", 
-                 hEventsTr, hEventsTrVtx));
-    list->ls();
-    return;
-  }
-#endif
   AliInfo("Fitting energy loss spectra");
   fEnergyFitter.Fit(list);
 
@@ -298,17 +286,6 @@ AliFMDEnergyFitterTask::Terminate(Option_t*)
   TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults", 
                                                      list->GetName())));
   PostData(2, list2);
-#if 0
-  // Temporary code to save output to a file - should be disabled once 
-  // the plugin stuff works 
-  list->ls();
-  TDirectory* savdir = gDirectory;
-  TFile* tmp = TFile::Open("elossfits.root", "RECREATE");
-  list->Write(list->GetName(), TObject::kSingleKey);
-  tmp->Write();
-  tmp->Close();
-  savdir->cd();
-#endif 
 }
 
 //____________________________________________________________________
index 3b5153827d3e9d14a036725dc899dbebb779c521..ece4bb3f1be785517d1ad55fd1fe707fd2cdbe4e 100644 (file)
@@ -37,6 +37,7 @@ AliFMDEventInspector::AliFMDEventInspector()
   : TNamed(),
     fHEventsTr(0), 
     fHEventsTrVtx(0),
+    fHEventsAccepted(0),
     fHTriggers(0),
     fHType(0),
     fHWords(0),
@@ -49,7 +50,8 @@ AliFMDEventInspector::AliFMDEventInspector()
     fField(999), 
     fCollisionSystem(kUnknown),
     fDebug(0),
-    fCentAxis(0)
+    fCentAxis(0),
+    fVtxAxis(10,-10,10)
 {
   // 
   // Constructor 
@@ -61,6 +63,7 @@ AliFMDEventInspector::AliFMDEventInspector(const char* name)
   : TNamed("fmdEventInspector", name),
     fHEventsTr(0), 
     fHEventsTrVtx(0), 
+    fHEventsAccepted(0),
     fHTriggers(0),
     fHType(0),
     fHWords(0),
@@ -73,7 +76,8 @@ AliFMDEventInspector::AliFMDEventInspector(const char* name)
     fField(999), 
     fCollisionSystem(kUnknown),
     fDebug(0),
-    fCentAxis(0)
+    fCentAxis(0),
+    fVtxAxis(10,-10,10)
 {
   // 
   // Constructor 
@@ -88,6 +92,7 @@ AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
   : TNamed(o), 
     fHEventsTr(o.fHEventsTr), 
     fHEventsTrVtx(o.fHEventsTrVtx), 
+    fHEventsAccepted(o.fHEventsAccepted),
     fHTriggers(o.fHTriggers),
     fHType(o.fHType),
     fHWords(o.fHWords),
@@ -100,7 +105,8 @@ AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
     fField(o.fField), 
     fCollisionSystem(o.fCollisionSystem),
     fDebug(0),
-    fCentAxis(0)
+    fCentAxis(0),
+    fVtxAxis(o.fVtxAxis)
 {
   // 
   // Copy constructor 
@@ -134,6 +140,7 @@ AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
   TNamed::operator=(o);
   fHEventsTr         = o.fHEventsTr;
   fHEventsTrVtx      = o.fHEventsTrVtx;
+  fHEventsAccepted   = o.fHEventsAccepted;
   fHTriggers         = o.fHTriggers;
   fHType             = o.fHType;
   fHWords            = o.fHWords;
@@ -146,6 +153,8 @@ AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
   fEnergy            = o.fEnergy;
   fField             = o.fField;
   fCollisionSystem   = o.fCollisionSystem;
+  fVtxAxis.Set(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(), 
+              o.fVtxAxis.GetXmax());
   if (fList) { 
     fList->SetName(GetName());
     if (fHEventsTr)    fList->Add(fHEventsTr);
@@ -206,12 +215,14 @@ AliFMDEventInspector::Init(const TAxis& vtxAxis)
   // ----- 92 number --------- ---- 1 ---
   TArrayD limits(93);
   for (Int_t i = 0; i < 92; i++) limits[i] = -1.5 + i;
+
+  fVtxAxis.Set(vtxAxis.GetNbins(), vtxAxis.GetXmin(), vtxAxis.GetXmax());
   
   fCentAxis  = new TAxis(limits.GetSize()-1, limits.GetArray());
   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger", 
-                       vtxAxis.GetNbins(), 
-                       vtxAxis.GetXmin(), 
-                       vtxAxis.GetXmax());
+                       4*vtxAxis.GetNbins(), 
+                       2*vtxAxis.GetXmin(), 
+                       2*vtxAxis.GetXmax());
   fHEventsTr->SetXTitle("v_{z} [cm]");
   fHEventsTr->SetYTitle("# of events");
   fHEventsTr->SetFillColor(kRed+1);
@@ -220,19 +231,26 @@ AliFMDEventInspector::Init(const TAxis& vtxAxis)
   // fHEventsTr->Sumw2();
   fList->Add(fHEventsTr);
 
-  fHEventsTrVtx = new TH1I("nEventsTrVtx", 
-                          "Number of events w/trigger and vertex", 
-                          vtxAxis.GetNbins(), 
-                          vtxAxis.GetXmin(), 
-                          vtxAxis.GetXmax());
-  fHEventsTrVtx->SetXTitle("v_{z} [cm]");
-  fHEventsTrVtx->SetYTitle("# of events");
+  fHEventsTrVtx = static_cast<TH1I*>(fHEventsTr->Clone("nEventsTrVtx"));
+  fHEventsTrVtx->SetTitle("Number of events w/trigger and vertex"); 
   fHEventsTrVtx->SetFillColor(kBlue+1);
-  fHEventsTrVtx->SetFillStyle(3001);
   fHEventsTrVtx->SetDirectory(0);
   // fHEventsTrVtx->Sumw2();
   fList->Add(fHEventsTrVtx);
 
+  fHEventsAccepted = new TH1I("nEventsAccepted", 
+                             "Number of events  w/trigger and vertex in range",
+                             2*vtxAxis.GetNbins(), 
+                             2*vtxAxis.GetXmin(), 
+                             2*vtxAxis.GetXmax());
+  fHEventsAccepted->SetXTitle("v_{z} [cm]");
+  fHEventsAccepted->SetYTitle("# of events");
+  fHEventsAccepted->SetFillColor(kGreen+1);
+  fHEventsAccepted->SetFillStyle(3001);
+  fHEventsAccepted->SetDirectory(0);
+  // fHEventsAccepted->Sumw2();
+  fList->Add(fHEventsAccepted);
+                             
       
   fHTriggers = new TH1I("triggers", "Triggers", kOffline+1, 0, kOffline+1);
   fHTriggers->SetFillColor(kRed+1);
@@ -316,7 +334,8 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
                              Bool_t&            lowFlux,
                              UShort_t&          ivz, 
                              Double_t&          vz,
-                             Double_t&          cent)
+                             Double_t&          cent,
+                             UShort_t&          nClusters)
 {
   // 
   // Process the event 
@@ -342,7 +361,7 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
   }
 
   // --- Read trigger information from the ESD and store in AOD object
-  if (!ReadTriggers(event, triggers)) { 
+  if (!ReadTriggers(event, triggers, nClusters)) { 
     if (fDebug > 2) {
       AliWarning("Failed to read triggers from ESD"); }
     return kNoTriggers;
@@ -386,15 +405,16 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
   fHEventsTrVtx->Fill(vz);
 
   // --- Get the vertex bin ------------------------------------------
-  ivz = fHEventsTr->GetXaxis()->FindBin(vz);
-  if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { 
+  ivz = fVtxAxis.FindBin(vz);
+  if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) { 
     if (fDebug > 3) {
       AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
-                     vz, fHEventsTr->GetXaxis()->GetXmin(), 
-                     fHEventsTr->GetXaxis()->GetXmax())); }
+                     vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); 
+    }
     ivz = 0;
     return kBadVertex;
   }
+  fHEventsAccepted->Fill(vz);
   
   // --- Check the FMD ESD data --------------------------------------
   if (!event->GetFMDData()) { 
@@ -439,7 +459,8 @@ AliFMDEventInspector::ReadCentrality(const AliESDEvent* esd,
 
 //____________________________________________________________________
 Bool_t
-AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
+AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers,
+                                  UShort_t& nClusters)
 {
   // 
   // Read the trigger information from the ESD event 
@@ -475,6 +496,7 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
   // then the AliPhysicsSelection object would overcount by a 
   // factor of 2! :-(
   Bool_t offline = ih->IsEventSelected();
+  nClusters = 0;
   if (offline) {
     triggers |= AliAODForwardMult::kOffline;
     triggers |= AliAODForwardMult::kInel;
@@ -489,14 +511,22 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
       // Check if we have one or more tracklets 
       // in the range -1 < eta < 1 to set the INEL>0 
       // trigger flag. 
+      // 
+      // Also count tracklets as a single cluster 
       Int_t n = spdmult->GetNumberOfTracklets();
       for (Int_t j = 0; j < n; j++) { 
        if(TMath::Abs(spdmult->GetEta(j)) < 1) { 
          triggers |= AliAODForwardMult::kInelGt0;
-         break;
+         nClusters++;
        }
       }
+      n = spdmult->GetNumberOfSingleClusters();
+      for (Int_t j = 0; j < n; j++) { 
+       Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
+       if (TMath::Abs(eta) < 1) nClusters++;
+      }
     }
+    if (nClusters > 0) triggers |= AliAODForwardMult::kNClusterGt0;
   }
 
   // Analyse some trigger stuff 
@@ -505,7 +535,8 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
     triggers |= AliAODForwardMult::kNSD;
  
     
-  //Check pileup
+  // Check for multiple vertices (pile-up) with at least 3
+  // contributors and at least 0.8cm from the primary vertex
   Bool_t pileup =  esd->IsPileupFromSPD(3,0.8);
   if (pileup) {
     triggers |= AliAODForwardMult::kPileUp;
@@ -628,6 +659,40 @@ AliFMDEventInspector::ReadVertex(const AliESDEvent* esd, Double_t& vz)
   //    @c true on success, @c false otherwise 
   //
   vz = 0;
+#if 1
+  // This is the code used by the 1st physics people 
+  const AliESDVertex* vertex    = esd->GetPrimaryVertex();
+  if (!vertex  || !vertex->GetStatus()) {
+    if (fDebug > 2) {
+      AliWarning(Form("No primary vertex (%p) or bad status %d", 
+                     vertex, (vertex ? vertex->GetStatus() : -1)));
+    }
+    return false;
+  }
+  const AliESDVertex* vertexSPD = esd->GetPrimaryVertexSPD();
+  if (!vertexSPD || !vertexSPD->GetStatus()) {
+    if (fDebug > 2) {
+      AliWarning(Form("No primary SPD vertex (%p) or bad status %d", 
+                     vertexSPD, (vertexSPD ? vertexSPD->GetStatus() : -1)));
+    }
+    return false;
+  }
+
+  // if vertex is from SPD vertexZ, require more stringent cuts 
+  if (vertex->IsFromVertexerZ()) { 
+    if (vertex->GetDispersion() > fMaxVzErr || 
+       vertex->GetZRes() > 1.25 * fMaxVzErr) {
+      if (fDebug > 2) {
+       AliWarning(Form("Dispersion %f > %f or resolution %f > %f",
+                       vertex->GetDispersion(), fMaxVzErr,
+                       vertex->GetZRes(), 1.25 * fMaxVzErr)); 
+      }
+      return false;
+    }
+  }
+  vz = vertex->GetZ();
+  return true;
+#else 
   // Get the vertex 
   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
   if (!vertex) { 
@@ -655,6 +720,7 @@ AliFMDEventInspector::ReadVertex(const AliESDEvent* esd, Double_t& vz)
   // Get the z coordiante 
   vz = vertex->GetZ();
   return kTRUE;
+#endif 
 }
   
 //____________________________________________________________________
index ec91eb6a7e1b89aea0c392bc8d5ab301a0ff132e..c526202030083e9c8acf50e1b5e4259fb56d403b 100644 (file)
@@ -14,6 +14,7 @@
  * @ingroup pwg2_forward_aod
  */
 #include <TNamed.h>
+#include <TAxis.h>
 class AliESDEvent;
 class TH2D;
 class TH1D;
@@ -135,6 +136,7 @@ public:
    * @param vz        On return, the z position of the interaction
    * @param cent      On return, the centrality (in percent) or < 0 
    *                  if not found
+   * @param nClusters On return, number of SPD clusters in @f$ |\eta|<1@f$ 
    * 
    * @return 0 (or kOk) on success, otherwise a bit mask of error codes 
    */
@@ -143,7 +145,8 @@ public:
                 Bool_t&            lowFlux,
                 UShort_t&          ivz, 
                 Double_t&          vz,
-                Double_t&          cent);
+                Double_t&          cent,
+                UShort_t&          nClusters);
   /** 
    * Define the output histograms.  These are put in a sub list of the
    * passed list.   The histograms are merged before the parent task calls 
@@ -224,10 +227,12 @@ protected:
    * 
    * @param esd        ESD event 
    * @param triggers   On return, contains the trigger bits 
+   * @param nClusters  On return, number of SPD clusters in @f$ |\eta|<1@f$ 
    * 
    * @return @c true on success, @c false otherwise 
    */
-  Bool_t ReadTriggers(const AliESDEvent* esd, UInt_t& triggers);
+  Bool_t ReadTriggers(const AliESDEvent* esd, UInt_t& triggers, 
+                     UShort_t& nClusters);
   /** 
    * Read the vertex information from the ESD event 
    * 
@@ -251,6 +256,7 @@ protected:
 
   TH1I*    fHEventsTr;    //! Histogram of events w/trigger
   TH1I*    fHEventsTrVtx; //! Events w/trigger and vertex 
+  TH1I*    fHEventsAccepted; //! Events w/trigger and vertex in range 
   TH1I*    fHTriggers;    //! Triggers
   TH1I*    fHType;        //! Type (low/high flux) of event
   TH1I*    fHWords;       //! Trigger words 
@@ -264,6 +270,7 @@ protected:
   UShort_t fCollisionSystem; //  Collision system
   Int_t    fDebug;        //  Debug level 
   TAxis*   fCentAxis;     // Centrality axis used in histograms
+  TAxis    fVtxAxis;
   ClassDef(AliFMDEventInspector,3); // Inspect the event 
 };
 
index b4c4d50a470a6df3f314c2c2e4e032d7a8028984..f3b4d4ebc1cde90d4b3f4dc07ff29238565b7876 100644 (file)
@@ -33,6 +33,50 @@ ClassImp(AliFMDHistCollector)
 ; // For Emacs
 #endif 
 
+//____________________________________________________________________
+AliFMDHistCollector::AliFMDHistCollector()
+  : fNCutBins(0), 
+    fCorrectionCut(0), 
+    fFirstBins(), 
+    fLastBins(), 
+    fDebug(0),
+    fList(0),
+    fSumRings(0),
+    fCoverage(0),
+    fMergeMethod(kStraightMean),
+    fFiducialMethod(kByCut)
+{}
+
+//____________________________________________________________________
+AliFMDHistCollector::AliFMDHistCollector(const char* title)
+  : TNamed("fmdHistCollector", title), 
+    fNCutBins(2), 
+    fCorrectionCut(0.5), 
+    fFirstBins(1), 
+    fLastBins(1), 
+    fDebug(0),
+    fList(0),
+    fSumRings(0),
+    fCoverage(0),
+    fMergeMethod(kStraightMean),
+    fFiducialMethod(kByCut)
+{
+}
+//____________________________________________________________________
+AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
+  : TNamed(o), 
+    fNCutBins(o.fNCutBins), 
+    fCorrectionCut(o.fCorrectionCut),
+    fFirstBins(o.fFirstBins), 
+    fLastBins(o.fLastBins), 
+    fDebug(o.fDebug),
+    fList(o.fList),
+    fSumRings(o.fSumRings),
+    fCoverage(o.fCoverage),
+    fMergeMethod(o.fMergeMethod),
+    fFiducialMethod(o.fFiducialMethod)
+{}
+
 //____________________________________________________________________
 AliFMDHistCollector::~AliFMDHistCollector()
 { 
@@ -63,6 +107,7 @@ AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
   fCoverage       = o.fCoverage;
   fMergeMethod    = o.fMergeMethod;
   fFiducialMethod = o.fFiducialMethod;
+
   return *this;
 }
 
@@ -151,7 +196,7 @@ AliFMDHistCollector::Init(const TAxis& vtxAxis,
       // Store the result for later use 
       fFirstBins[(iVz-1)*5+iIdx] = first;
       fLastBins[(iVz-1)*5+iIdx]  = last;
-      TH2D* obg = static_cast<TH2D*>(bg->Clone());
+      TH2D* obg = static_cast<TH2D*>(bg->Clone(Form("secMapFMD%d%c", d, r)));
       obg->SetDirectory(0);
       obg->Reset();
       vtxList->Add(obg);
@@ -212,8 +257,10 @@ AliFMDHistCollector::DefineOutput(TList* dir)
   fList->SetOwner();
   fList->SetName(GetName());
   dir->Add(fList);
+
 }
 
+
 //____________________________________________________________________
 Int_t
 AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
@@ -477,7 +524,8 @@ AliFMDHistCollector::MergeBins(Double_t c,   Double_t e,
 
 //____________________________________________________________________
 Bool_t
-AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
+AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists,
+                            AliForwardUtil::Histos& sums,
                             UShort_t                vtxbin, 
                             TH2D&                   out)
 {
@@ -499,13 +547,32 @@ AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
       TH2D*       h = hists.Get(d,r);
       TH2D*       t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
       Int_t       i = (d == 1 ? 1 : 2*d + (q == 0 ? -2 : -1));
+      TH2D*       o = sums.Get(d, r);
 
-      // Outer rings have better phi segmentation - rebin to same as inner. 
-      if (q == 1) t->RebinY(2);
-
+      // Get valid range 
       Int_t first = 0;
       Int_t last  = 0;
       GetFirstAndLast(d, r, vtxbin, first, last);
+      
+      // Zero outside valid range 
+      for (Int_t iPhi = 0; iPhi <= t->GetNbinsY()+1; iPhi++) { 
+       // Lower range 
+       for (Int_t iEta = 1; iEta < first; iEta++) { 
+         t->SetBinContent(iEta,iPhi,0);
+         t->SetBinError(iEta,iPhi,0);
+       }
+       for (Int_t iEta = last+1; iEta <= t->GetNbinsX(); iEta++) {
+         t->SetBinContent(iEta,iPhi,0);
+         t->SetBinError(iEta,iPhi,0);
+       }
+      }
+      for (Int_t iEta = first; iEta <= last; iEta++)
+       t->SetBinContent(iEta,0,1);
+      // Add to our per-ring sum 
+      o->Add(t);
+      
+      // Outer rings have better phi segmentation - rebin to same as inner. 
+      if (q == 1) t->RebinY(2);
 
       // Now update profile output 
       for (Int_t iEta = first; iEta <= last; iEta++) { 
@@ -518,7 +585,8 @@ AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
        Float_t noc      = overlap >= 0? 0.5 : 1;
        out.SetBinContent(iEta, 0, ooc + noc);
 
-       for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) { 
+       // Should we loop over h or t Y bins - I think it's t
+       for (Int_t iPhi = 1; iPhi <= t->GetNbinsY(); iPhi++) { 
          Double_t c  = t->GetBinContent(iEta,iPhi);
          Double_t e  = t->GetBinError(iEta,iPhi);
          Double_t ee = t->GetXaxis()->GetBinCenter(iEta);
@@ -577,22 +645,26 @@ AliFMDHistCollector::Print(Option_t* /* option */) const
   case kLeastError:         std::cout << "least error\n"; break;
   }
     
-  std::cout << ind << " Bin ranges:\n" << ind << "  v_z bin";
+  std::cout << ind << " Bin ranges:\n" << ind << "  rings  ";
   Int_t nVz = fFirstBins.fN / 5;
-  for (UShort_t iVz = 1; iVz <= nVz; iVz++) 
-    std::cout << " | " << std::setw(7) << iVz;
-  std::cout << '\n' << ind << "  / ring ";
-  for (UShort_t iVz = 1; iVz <= nVz; iVz++) 
-    std::cout << "-+--------";
-  std::cout << std::endl;
-    
   for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
     UShort_t d = 0;
     Char_t   r = 0;
     GetDetRing(iIdx, d, r);
+    std::cout << ind << " | FMD" << d << r << " ";
+  }
+  std::cout << '\n' << ind << "  /vz_bin ";
+  for (Int_t iIdx = 0; iIdx < 5; iIdx++) 
+    std::cout << "-+--------";
+  std::cout << std::endl;
+
+  for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
+    std::cout << " " << std::setw(7) << iVz << "   ";
+    for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
+      UShort_t d = 0;
+      Char_t   r = 0;
+      GetDetRing(iIdx, d, r);
     
-    std::cout << ind << "  FMD" << d << r << "  ";
-    for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
       Int_t first, last;
       GetFirstAndLast(iIdx, iVz, first, last);
       std::cout << " | " << std::setw(3) << first << "-" 
index 56634c88e228ed8f788d5979b62e44b21b0e0c9f..3b0e0898cd9dfa13897ffa1a341339b5df77662d 100644 (file)
@@ -97,54 +97,19 @@ public:
   /** 
    * Constructor 
    */
-  AliFMDHistCollector() 
-    : fNCutBins(0), 
-      fCorrectionCut(0), 
-      fFirstBins(), 
-      fLastBins(), 
-      fDebug(0),
-      fList(0),
-      fSumRings(0),
-      fCoverage(0),
-      fMergeMethod(kStraightMean),
-      fFiducialMethod(kByCut)
-  {}
+  AliFMDHistCollector();
   /** 
    * Constructor 
    * 
    * @param title Name of object
    */
-  AliFMDHistCollector(const char* title)
-    : TNamed("fmdHistCollector", title), 
-      fNCutBins(2), 
-      fCorrectionCut(0.5), 
-      fFirstBins(1), 
-      fLastBins(1), 
-      fDebug(0),
-      fList(0),
-      fSumRings(0),
-      fCoverage(0),
-      fMergeMethod(kStraightMean),
-      fFiducialMethod(kByCut)
-  {}
+  AliFMDHistCollector(const char* title);
   /** 
    * Copy constructor 
    * 
    * @param o Object to copy from 
    */
-  AliFMDHistCollector(const AliFMDHistCollector& o)
-    : TNamed(o), 
-      fNCutBins(o.fNCutBins), 
-      fCorrectionCut(o.fCorrectionCut),
-      fFirstBins(o.fFirstBins), 
-      fLastBins(o.fLastBins), 
-      fDebug(o.fDebug),
-      fList(o.fList),
-      fSumRings(o.fSumRings),
-      fCoverage(o.fCoverage),
-      fMergeMethod(o.fMergeMethod),
-      fFiducialMethod(o.fFiducialMethod)
-  {}
+  AliFMDHistCollector(const AliFMDHistCollector& o);
 
   /** 
    * Destructor 
@@ -170,13 +135,16 @@ public:
    * Do the calculations 
    * 
    * @param hists    Cache of histograms 
+   * @param sums     Cache to sum ring histograms in 
    * @param vtxBin   Vertex bin (1 based)
    * @param out      Output histogram
    * 
    * @return true on successs 
    */
-  virtual Bool_t Collect(AliForwardUtil::Histos& hists, UShort_t vtxBin, 
-                        TH2D& out);
+  virtual Bool_t Collect(const AliForwardUtil::Histos& hists, 
+                        AliForwardUtil::Histos&       sums, 
+                        UShort_t                      vtxBin, 
+                        TH2D&                         out);
   /** 
    * Output diagnostic histograms to directory 
    * 
@@ -368,6 +336,7 @@ protected:
   void MergeBins(Double_t c,   Double_t e, 
                 Double_t oc,  Double_t oe,
                 Double_t& rc, Double_t& re) const;
+  
 
   Int_t       fNCutBins;        // Number of additional bins to cut away
   Float_t     fCorrectionCut;   // Cut-off on secondary corrections 
@@ -379,6 +348,7 @@ protected:
   TH2D*       fCoverage;        // Sum per ring (on y-axis)
   MergeMethod fMergeMethod;     // Merge methiod for overlapping bins 
   FiducialMethod fFiducialMethod; // Fidicual method
+
   ClassDef(AliFMDHistCollector,1); // Calculate Nch density 
 };
 
index d56a10bffb42555be6e42fb8a95f77f11e8b6e18..e12ec02d19710088177e9972e6dd267ba1c02657 100644 (file)
@@ -61,7 +61,7 @@ AliFMDMCCorrector::operator=(const AliFMDMCCorrector& o)
   //    Reference to this object
   //
   AliFMDCorrector::operator=(o);
-
+  fSecondaryForMC = o.fSecondaryForMC;
   return *this;
 }
 
@@ -80,6 +80,9 @@ AliFMDMCCorrector::CorrectMC(AliForwardUtil::Histos& hists,
   // Return:
   //    true on successs 
   //
+  if ((!fUseSecondaryMap || !fSecondaryForMC) && fUseVertexBias) 
+    return kTRUE;
+
   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 
   UShort_t uvb = vtxbin;
@@ -90,7 +93,7 @@ AliFMDMCCorrector::CorrectMC(AliForwardUtil::Histos& hists,
       TH2D*       h  = hists.Get(d,r);
 
 
-      if (fUseSecondaryMap) {
+      if (fUseSecondaryMap && fSecondaryForMC) {
         TH2D*       bg = fcm.GetSecondaryMap()->GetCorrection(d,r,uvb);
         if (!bg) {
           AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
@@ -247,6 +250,23 @@ AliFMDMCCorrector::DefineOutput(TList* dir)
   fComps->SetName("esd_mc_comparison");
   d->Add(fComps);
 }
+//____________________________________________________________________
+void
+AliFMDMCCorrector::Print(Option_t* option) const
+{
+  // 
+  // Print information
+  // Parameters:
+  //    option Not used 
+  //
+  char ind[gROOT->GetDirLevel()+1];
+  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+  ind[gROOT->GetDirLevel()] = '\0';
+  AliFMDCorrector::Print(option);
+  std::cout << std::boolalpha
+            << ind << " Use sec. map on MC:     " << fSecondaryForMC 
+            << std::noboolalpha << std::endl;
+}
 
 //____________________________________________________________________
 //
index a0b2d8e9194af8761873fa0dba80270df06670ed..7954b452d8723dd096647fa93d917b588a6abe17 100644 (file)
@@ -58,7 +58,8 @@ public:
       fFMD2o(0),
       fFMD3i(0),
       fFMD3o(0),
-      fComps(0)
+      fComps(0),
+      fSecondaryForMC(true)
   {}
   /** 
    * Constructor 
@@ -72,7 +73,8 @@ public:
       fFMD2o(0),
       fFMD3i(0),
       fFMD3o(0),
-      fComps(0)
+      fComps(0),
+      fSecondaryForMC(true)
   {}
   /** 
    * Copy constructor 
@@ -86,7 +88,8 @@ public:
       fFMD2o(o.fFMD2o),
       fFMD3i(o.fFMD3i),
       fFMD3o(o.fFMD3o),
-      fComps(0)
+      fComps(0),
+      fSecondaryForMC(o.fSecondaryForMC)
   {}
   /** 
    * Destructor 
@@ -100,6 +103,12 @@ public:
    * @return Reference to this object
    */
   AliFMDMCCorrector& operator=(const AliFMDMCCorrector&);
+  /** 
+   * If set, then do not do the secondary correction for MC data
+   * 
+   * @param use 
+   */
+  void SetSecondaryForMC(Bool_t use) { fSecondaryForMC = use; }
   /** 
    * Initialize this object 
    * 
@@ -133,6 +142,13 @@ public:
    * @param dir List to write in
    */  
   void DefineOutput(TList* dir);
+
+  /**
+   * Print information
+   * 
+   * @param option Not used 
+   */
+  void Print(Option_t* option="") const;
 protected:
   /** 
    * MAke comparison profiles
@@ -160,7 +176,8 @@ protected:
   TProfile2D* fFMD3i; // Comparison
   TProfile2D* fFMD3o; // Comparison
   TList*      fComps; // List of comparisons 
-  
+  Bool_t      fSecondaryForMC;  // Whether to correct MC data 
+
   ClassDef(AliFMDMCCorrector,1); // Calculate Nch density 
 };
 
index edaf4dd1ee2b26c05cb379747f67e129af93d10f..e933a92d9eb2b28b1fd2c68e2060507e73542c42 100644 (file)
@@ -456,7 +456,7 @@ AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
   else 
     return false;
   
-  // Invariant masses 
+  // Rapidity shift
   Double_t m02s = 1 - 2 * p1->Energy() / fEnergy; 
   Double_t m12s = 1 - 2 * p2->Energy() / fEnergy;
   
index db3c03b35e5aa1556bcaf7b12be840235c8171bf..cc2e86664b73d45bd436c823c34b27cfc13a27e6 100644 (file)
@@ -10,7 +10,7 @@
 //    - AliESDFMD object  - copy of input, but with signals merged 
 //
 // Corrections used: 
-//    - None
+//    - ELoss fits
 //
 // Histograms: 
 //    - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of 
@@ -28,6 +28,7 @@
 #include <TH1.h>
 #include <TMath.h>
 #include "AliFMDStripIndex.h"
+#include "AliFMDFloatMap.h"
 #include <AliLog.h>
 #include <TROOT.h>
 #include <iostream>
@@ -47,7 +48,10 @@ AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
     fFMD2o(0),
     fFMD3i(0),
     fFMD3o(0),
-    fSumEta(0)
+    fSumEta(0),
+    fOperComp(0),
+    fThetaVsNr(0), 
+    fOnlyPrimary(false)
 {
   // 
   // Constructor 
@@ -84,7 +88,25 @@ AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
   fSumEta->SetMarkerStyle(22);
   fSumEta->SetFillColor(0);
   fSumEta->SetFillStyle(0);
+
+  fOper     = new AliFMDFloatMap(0,0,0,0);
+  fOperComp = new TH2I("operComp", "Operation vs # track refs", 
+                      kMergedInto, kNone-.5, kMergedInto+.5, 
+                      20, -.5, 19.5);
+  fOperComp->SetXTitle("Operation");
+  fOperComp->SetYTitle("# of track refs in sector");
+  fOperComp->SetZTitle("Observations");
+  fOperComp->GetXaxis()->SetBinLabel(kNone,            "None");
+  fOperComp->GetXaxis()->SetBinLabel(kCandidate,       "Candidate");
+  fOperComp->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
+  fOperComp->GetXaxis()->SetBinLabel(kMergedInto,      "Merged into");
+  fOperComp->SetDirectory(0);
   
+  fThetaVsNr = new TH2D("thetaVsNr", "#theta of track vs # track references",
+                       360, 0, 360, 20, -.5, 19.5);
+  fThetaVsNr->SetXTitle("#theta [degrees]");
+  fThetaVsNr->SetYTitle("# of track references");
+  fThetaVsNr->SetDirectory(0);
 }
 
 //____________________________________________________________________
@@ -95,7 +117,10 @@ AliFMDMCSharingFilter::AliFMDMCSharingFilter(const AliFMDMCSharingFilter& o)
     fFMD2o(o.fFMD2o),
     fFMD3i(o.fFMD3i),
     fFMD3o(o.fFMD3o),
-    fSumEta(o.fSumEta)
+    fSumEta(o.fSumEta),
+    fOperComp(o.fOperComp),
+    fThetaVsNr(o.fThetaVsNr),
+    fOnlyPrimary(o.fOnlyPrimary)
 {
   // 
   // Copy constructor 
@@ -133,6 +158,7 @@ AliFMDMCSharingFilter::operator=(const AliFMDMCSharingFilter& o)
   //    Reference to this 
   //
   AliFMDSharingFilter::operator=(o);
+  fOnlyPrimary = o.fOnlyPrimary;
   return *this;
 }
 
@@ -142,6 +168,8 @@ AliFMDMCSharingFilter::StoreParticle(UShort_t   d,
                                     Char_t     r,
                                     UShort_t   s, 
                                     UShort_t   t, 
+                                    UShort_t   nR,
+                                    Double_t   theta,
                                     AliESDFMD& output) const
 {
   // 
@@ -152,10 +180,15 @@ AliFMDMCSharingFilter::StoreParticle(UShort_t   d,
   //    r       Ring
   //    s       Sector
   //    t       Strip
+  //    nR      Number of references to this particle in this sector
   //    output  Output ESD object
   //
   Double_t old = output.Multiplicity(d,r,s,t);
   if (old == AliESDFMD::kInvalidMult) old = 0;
+  if (fOper) fOperComp->Fill(fOper->operator()(d,r,s,t), nR);
+  if (theta < 0) theta += 2*TMath::Pi();
+  theta *= 180. / TMath::Pi();
+  fThetaVsNr->Fill(theta, nR);
   output.SetMultiplicity(d,r,s,t,old+1);
 }
 
@@ -183,6 +216,10 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input,
   //
   output.Clear();
 
+  // Increase event count - stored in 
+  // underflow bin 
+  fSumEta->AddBinContent(0); 
+
   // Copy eta values to output 
   for (UShort_t ed = 1; ed <= 3; ed++) { 
     UShort_t nq = (ed == 1 ? 1 : 2);
@@ -210,17 +247,28 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input,
     if (!isCharged) continue;
     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
 
+    // Pseudo rapidity and azimuthal angle 
+    Double_t eta = particle->Eta();
+    Double_t phi = particle->Phi();
+    
     // Fill 'dn/deta' histogram 
     if (isPrimary && iTr < nPrim) { 
-      fSumEta->Fill(particle->Eta());
-      primary->Fill(particle->Eta(), particle->Phi());
+      // Avoid under count - used to store event count
+      if (eta >= fSumEta->GetXaxis()->GetXmin()) fSumEta->Fill(eta);
+      primary->Fill(eta, phi);
     }
 
+    // Bail out if we're only processing primaries - perhaps we should
+    // track back to the original primary?
+    if (fOnlyPrimary && !isPrimary) continue;
+
     Int_t    nTrRef  = particle->GetNumberOfTrackReferences();
     Int_t    longest = -1;
     Double_t angle   = 0;
     UShort_t oD = 0, oS = 1024, oT = 1024;
     Char_t   oR = '\0';
+    UShort_t nC = 0;
+    Double_t oTheta = 0;
     for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
       AliTrackReference* ref = particle->GetTrackReference(iTrRef);
       
@@ -231,13 +279,24 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input,
       if (ref->DetectorId() != AliTrackReference::kFMD) 
        continue;
 
+      // Count number of track refs in this sector 
+      nC++;
+
       // Get the detector coordinates 
       UShort_t d, s, t;
       Char_t r;
       AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
-      if (oD > 0 && oR != '\0' && (d != oD || r != oR)) {
+      // If this is a new detector/ring, then reset the other one 
+      if (oD > 0 && oR != '\0' && oS != 1024 && 
+         (d != oD || r != oR || s != oS)) {
        longest = -1;
-       StoreParticle(oD, oR, oS, oT, output);
+       angle   = 0;
+       StoreParticle(oD, oR, oS, oT, nC, oTheta, output);
+       nC = 0;
+       oD = 0;
+       oR = '\0';
+       oS = 1024;
+       oT = 1024;
       }
 
       oD = d;
@@ -256,6 +315,7 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input,
        longest = iTrRef;
        angle   = ang;
       }
+      oTheta = theta;
     } // Loop over track references
     if (longest < 0) continue;
 
@@ -267,7 +327,7 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input,
     Char_t r;
     AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
     
-    StoreParticle(d,r,s,t,output);
+    StoreParticle(d,r,s,t,nC,particle->Theta(),output);
   } // Loop over tracks
   return kTRUE;
 }
@@ -335,6 +395,8 @@ AliFMDMCSharingFilter::DefineOutput(TList* dir)
   cd->Add(fFMD3i);
   cd->Add(fFMD3o);
   dir->Add(fSumEta);
+  cd->Add(fOperComp);
+  cd->Add(fThetaVsNr);
 }
 
 //____________________________________________________________________
@@ -354,7 +416,8 @@ AliFMDMCSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents)
     AliWarning(Form("No mcSumEta histogram found in output list"));
     return;
   }
-  sumEta->Scale(1. / nEvents, "width");
+  Double_t n = nEvents; // sumEta->GetBinContent(0);
+  sumEta->Scale(1. / n, "width");
 }
 
 //____________________________________________________________________
@@ -371,6 +434,9 @@ AliFMDMCSharingFilter::Print(Option_t* option) const
   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
   ind[gROOT->GetDirLevel()] = '\0';
   AliFMDSharingFilter::Print(option);
+  std::cout << std::boolalpha 
+           << ind << " Only primary tracks:    " << fOnlyPrimary 
+           << std::noboolalpha << std::endl;
 }
 
 //____________________________________________________________________
index 894761c60c0b16b1a4f65345bd76fed669832732..09aa4e0724cc5f2ff7e01af49f30c493450b26e7 100644 (file)
@@ -58,7 +58,10 @@ public:
     fFMD2o(0),
     fFMD3i(0),
     fFMD3o(0), 
-    fSumEta(0)
+    fSumEta(0), 
+    fOperComp(0), 
+    fThetaVsNr(0), 
+    fOnlyPrimary(false)
   {}
   /** 
    * Constructor 
@@ -80,6 +83,13 @@ public:
    * @return Reference to this 
    */
   AliFMDMCSharingFilter& operator=(const AliFMDMCSharingFilter& o);
+  
+  /** 
+   * If set, then only process primary tracks 
+   * 
+   * @param use 
+   */
+  void SetOnlyPrimary(Bool_t use) { fOnlyPrimary = use; }
   /** 
    * Filter the input kinematics and track references, using 
    * some of the ESD information
@@ -138,14 +148,16 @@ protected:
    * @param output  Output ESD object
    */
   void StoreParticle(UShort_t d, Char_t r, UShort_t s, UShort_t t, 
-                    AliESDFMD& output) const;
+                    UShort_t nr, Double_t theta, AliESDFMD& output) const;
   TH2D* fFMD1i;  // ESD-MC correlation 
   TH2D* fFMD2i;  // ESD-MC correlation 
   TH2D* fFMD2o;  // ESD-MC correlation 
   TH2D* fFMD3i;  // ESD-MC correlation 
   TH2D* fFMD3o;  // ESD-MC correlation 
   TH1D* fSumEta; // MC dN/deta 
-
+  TH2I* fOperComp; // Operation vs # trackrefs
+  TH2D* fThetaVsNr; // Theta vs # trackrefs
+  Bool_t fOnlyPrimary; // Only process primary tracks 
   ClassDef(AliFMDMCSharingFilter,1); //
 };
 
index 71aba79d3eb04d9d49ff78dc13a0f517021bccad..bfa24faeef56fa4e2c8247650696a14e5b0f0749 100644 (file)
@@ -31,6 +31,7 @@
 #include "AliFMDCorrELossFit.h"
 #include <AliLog.h>
 #include <TROOT.h>
+#include <THStack.h>
 #include <iostream>
 #include <iomanip>
 
@@ -39,6 +40,12 @@ ClassImp(AliFMDSharingFilter)
 ; // This is for Emacs
 #endif 
 
+#define DBG(L,M) \
+  do { if (L>fDebug)break; std::cout << (M) << std::flush;} while(false) 
+#define DBGL(L,M) \
+  do { if (L>fDebug)break; std::cout << (M) << std::endl;} while(false) 
+                             
+
 
 //____________________________________________________________________
 AliFMDSharingFilter::AliFMDSharingFilter()
@@ -47,6 +54,10 @@ AliFMDSharingFilter::AliFMDSharingFilter()
     fLowCut(0.),
     fCorrectAngles(kFALSE), 
     fNXi(1),
+    fIncludeSigma(true),
+    fSummed(0),
+    fHighCuts(0),
+    fOper(0),
     fDebug(0)
 {
   // 
@@ -61,6 +72,10 @@ AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
     fLowCut(0.),
     fCorrectAngles(kFALSE), 
     fNXi(1),
+    fIncludeSigma(true),
+    fSummed(0),
+    fHighCuts(0),
+    fOper(0),
     fDebug(0)
 {
   // 
@@ -83,6 +98,10 @@ AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
     fLowCut(o.fLowCut),
     fCorrectAngles(o.fCorrectAngles), 
     fNXi(o.fNXi),
+    fIncludeSigma(o.fIncludeSigma),
+    fSummed(o.fSummed),
+    fHighCuts(o.fHighCuts),
+    fOper(o.fOper),
     fDebug(o.fDebug)
 {
   // 
@@ -124,6 +143,10 @@ AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
   fCorrectAngles = o.fCorrectAngles;
   fNXi           = o.fNXi;
   fDebug         = o.fDebug;
+  fOper          = o.fOper;
+  fSummed        = o.fSummed;
+  fHighCuts      = o.fHighCuts;
+  fIncludeSigma  = o.fIncludeSigma;
 
   fRingHistos.Delete();
   TIter    next(&o.fRingHistos);
@@ -158,6 +181,43 @@ AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
   return static_cast<RingHistos*>(fRingHistos.At(idx));
 }
 
+//____________________________________________________________________
+void
+AliFMDSharingFilter::Init()
+{
+  // Initialise 
+  AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
+  // Get the high cut.  The high cut is defined as the 
+  // most-probably-value peak found from the energy distributions, minus 
+  // 2 times the width of the corresponding Landau.
+  AliFMDCorrELossFit* fits = fcm.GetELossFit();
+  const TAxis& eAxis = fits->GetEtaAxis();
+
+  UShort_t nEta = eAxis.GetNbins();
+  fHighCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5);
+  fHighCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
+  fHighCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
+  fHighCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
+  fHighCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
+  fHighCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
+
+  UShort_t ybin = 0;
+  for (UShort_t d = 1; d <= 3; d++) {
+    UShort_t nr = (d == 1 ? 1 : 2);
+    for (UShort_t q = 0; q < nr; q++) { 
+      Char_t r = (q == 0 ? 'I' : 'O');
+      ybin++;
+      for (UShort_t e = 1; e <= nEta; e++) { 
+       Double_t eta = eAxis.GetBinCenter(e);
+       Double_t cut = GetHighCut(d, r, eta, false);
+       if (cut <= 0) continue;
+       fHighCuts->SetBinContent(e, ybin, cut);
+      }
+    }
+  }
+}
+
 //____________________________________________________________________
 Bool_t
 AliFMDSharingFilter::Filter(const AliESDFMD& input, 
@@ -178,45 +238,59 @@ AliFMDSharingFilter::Filter(const AliESDFMD& input,
   output.Clear();
   TIter    next(&fRingHistos);
   RingHistos* o      = 0;
-  Double_t    lowCut = GetLowCut();
   while ((o = static_cast<RingHistos*>(next())))
     o->Clear();
 
-  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+  if (fOper) fOper->Reset(0);
+  Int_t nNone      = 0;
+  Int_t nCandidate = 0;
+  Int_t nMerged    = 0;
+  Int_t nSummed    = 0;
 
+  Status status[512];
+  
   for(UShort_t d = 1; d <= 3; d++) {
     Int_t nRings = (d == 1 ? 1 : 2);
     for (UShort_t q = 0; q < nRings; q++) {
-      Char_t      r    = (q == 0 ? 'I' : 'O');
-      UShort_t    nsec = (q == 0 ?  20 :  40);
-      UShort_t    nstr = (q == 0 ? 512 : 256);
-
+      Char_t      r      = (q == 0 ? 'I' : 'O');
+      UShort_t    nsec   = (q == 0 ?  20 :  40);
+      UShort_t    nstr   = (q == 0 ? 512 : 256);
       RingHistos* histos = GetRingHistos(d, r);
       
       for(UShort_t s = 0; s < nsec;  s++) {
-       Bool_t usedThis = kFALSE;
-       Bool_t usedPrev = kFALSE;
-       
+       for (UShort_t t = 0; t < nstr; t++) status[t] = kCandidate;
+#ifdef USE_OLDER_MERGING
+       Bool_t usedThis   = kFALSE;
+       Bool_t usedPrev   = kFALSE;
+#endif 
        for(UShort_t t = 0; t < nstr; t++) {
          output.SetMultiplicity(d,r,s,t,0.);
          Float_t mult = SignalInStrip(input,d,r,s,t);
+         // Get the pseudo-rapidity 
+         Double_t eta = input.Eta(d,r,s,t);
+         Double_t phi = input.Phi(d,r,s,t) * TMath::Pi() / 180.;
+         if (s == 0) output.SetEta(d,r,s,t,eta);
          
          // Keep dead-channel information. 
          if(mult == AliESDFMD::kInvalidMult)
            output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
 
          // If no signal or dead strip, go on. 
-         if (mult == AliESDFMD::kInvalidMult || mult == 0) continue;
+         if (mult == AliESDFMD::kInvalidMult || mult == 0) {
+           if (mult == 0) histos->fSum->Fill(eta,phi,mult);
+           status[t] = kNone;
+           continue;
+         }
 
-         // Get the pseudo-rapidity 
-         Double_t eta = input.Eta(d,r,s,t);
-         
          // Fill the diagnostics histogram 
          histos->fBefore->Fill(mult);
 
          // Get next and previous signal - if any 
          Double_t prevE = 0;
          Double_t nextE = 0;
+         Status   prevStatus = (t == 0      ? kNone : status[t-1]);
+         Status   thisStatus = status[t];
+         Status   nextStatus = (t == nstr-1 ? kNone : status[t+1]);
          if (t != 0) {
            prevE = SignalInStrip(input,d,r,s,t-1);
            if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
@@ -225,20 +299,60 @@ AliFMDSharingFilter::Filter(const AliESDFMD& input,
            nextE = SignalInStrip(input,d,r,s,t+1);
            if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
          }
+         if (t != 0) histos->fNeighborsBefore->Fill(prevE, mult);
          
+#ifdef USE_OLDER_MERGING
          Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
                                                      lowFlux,d,r,s,t, 
                                                      usedPrev,usedThis);
-         if (mergedEnergy > lowCut) histos->Incr();
+         status[t] = (usedPrev ? kMergedWithOther : kNone);
+         if (t != nstr - 1) status[t] = (usedThis ? kMergedWithOther : kNone);
+#else 
+         Double_t mergedEnergy = MultiplicityOfStrip(mult, prevE, nextE, 
+                                                     eta, lowFlux, 
+                                                     d, r, s, t, 
+                                                     prevStatus, 
+                                                     thisStatus, 
+                                                     nextStatus);
+         if (t != 0)      status[t-1] = prevStatus;
+         if (t != nstr-1) status[t+1] = nextStatus;
+         status[t] = thisStatus;
+
+#endif
+         // If we're processing on non-angle corrected data, we
+         // should do the angle correction here
+         if (!fCorrectAngles)
+           mergedEnergy = AngleCorrect(mergedEnergy, eta);
+         if (mergedEnergy > 0) histos->Incr();
+
+         if (t != 0) 
+           histos->fNeighborsAfter->Fill(output.Multiplicity(d,r,s,t-1), 
+                                         mergedEnergy);
+         histos->fBeforeAfter->Fill(mult, mergedEnergy);
          histos->fAfter->Fill(mergedEnergy);
+         histos->fSum->Fill(eta,phi,mergedEnergy);
 
          output.SetMultiplicity(d,r,s,t,mergedEnergy);
-         output.SetEta(d,r,s,t,eta);
        } // for strip
+       for (UShort_t t = 0; t < nstr; t++) {
+         if (fOper) fOper->operator()(d, r, s, t) = status[t];
+         switch (status[t]) { 
+         case kNone:            nNone++;      break;
+         case kCandidate:       nCandidate++; break;
+         case kMergedWithOther: nMerged++;    break;
+         case kMergedInto:      nSummed++;    break;
+         }
+       }
       } // for sector
     } // for ring 
   } // for detector
+  fSummed->Fill(kNone,            nNone);
+  fSummed->Fill(kCandidate,       nCandidate);
+  fSummed->Fill(kMergedWithOther, nMerged);
+  fSummed->Fill(kMergedInto,      nSummed);
 
+  DBGL(1, Form("none=%9d, candidate=%9d, merged=%9d, summed=%9d", 
+              nNone, nCandidate, nMerged, nSummed));
   next.Reset();
   while ((o = static_cast<RingHistos*>(next())))
     o->Finish();
@@ -268,17 +382,26 @@ AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
   //    The energy signal 
   //
   Double_t mult = input.Multiplicity(d,r,s,t);
-  if (mult == AliESDFMD::kInvalidMult || mult == 0) return mult;
-
-  if (fCorrectAngles && !input.IsAngleCorrected()) 
-    mult = AngleCorrect(mult, input.Eta(d,r,s,t));
-  else if (!fCorrectAngles && input.IsAngleCorrected()) 
-    mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
+  // In case of 
+  //  - bad value (invalid or 0) 
+  //  - we want angle corrected and data is 
+  //  - we don't want angle corrected and data isn't 
+  // just return read value  
+  if (mult == AliESDFMD::kInvalidMult               || 
+      mult == 0                                     ||
+      (fCorrectAngles && input.IsAngleCorrected()) || 
+      (!fCorrectAngles && !input.IsAngleCorrected()))
+    return mult;
+
+  // If we want angle corrected data, correct it, 
+  // otherwise de-correct it 
+  if (fCorrectAngles) mult = AngleCorrect(mult, input.Eta(d,r,s,t));
+  else                mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
   return mult;
 }
 //_____________________________________________________________________
 Double_t 
-AliFMDSharingFilter::GetLowCut() const
+AliFMDSharingFilter::GetLowCut(UShort_t, Char_t, Double_t) const
 {
   //
   // Get the low cut.  Normally, the low cut is taken to be the lower
@@ -294,7 +417,8 @@ AliFMDSharingFilter::GetLowCut() const
                        
 //_____________________________________________________________________
 Double_t 
-AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
+AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, 
+                               Double_t eta, Bool_t errors) const
 {
   //
   // Get the high cut.  The high cut is defined as the 
@@ -307,22 +431,130 @@ AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
   // Get the high cut.  The high cut is defined as the 
   // most-probably-value peak found from the energy distributions, minus 
   // 2 times the width of the corresponding Landau.
-  AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
-  Double_t delta = 1024;
-  Double_t xi    = 1024;
-  if (!fit) {
-    AliError(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
+  AliFMDCorrELossFit* fits = fcm.GetELossFit();
+  
+  return fits->GetLowerBound(d, r, eta, fNXi, errors, fIncludeSigma);
+}
+
+//_____________________________________________________________________
+namespace { 
+  const char* status2String(AliFMDSharingFilter::Status s)
+  {
+    switch(s) { 
+    case AliFMDSharingFilter::kNone:            return "none     "; 
+    case AliFMDSharingFilter::kCandidate:       return "candidate"; 
+    case AliFMDSharingFilter::kMergedWithOther: return "merged   "; 
+    case AliFMDSharingFilter::kMergedInto:      return "summed   "; 
+    }
+    return "unknown  ";
   }
-  else {
-    delta = fit->fDelta;
-    xi    = fit->fXi;
+}
+
+//_____________________________________________________________________
+Double_t 
+AliFMDSharingFilter::MultiplicityOfStrip(Double_t thisE,
+                                        Double_t prevE,
+                                        Double_t nextE,
+                                        Double_t eta,
+                                        Bool_t   lowFlux,
+                                        UShort_t d,
+                                        Char_t   r,
+                                        UShort_t s,
+                                        UShort_t t,
+                                        Status&  prevStatus, 
+                                        Status&  thisStatus, 
+                                        Status&  nextStatus) const
+{
+  Double_t lowCut = GetLowCut(d, r, eta);
+
+  DBG(3,Form("FMD%d%c[%2d,%3d]: this=%9f(%s), prev=%9f(%s), next=%9f(%s)", 
+            d, r, s, t,
+            thisE, status2String(thisStatus), 
+            prevE, status2String(prevStatus), 
+            nextE, status2String(nextStatus)));
+
+  // If below cut, then modify zero signal and make sure the next
+  // strip is considered a candidate.
+  if (thisE < lowCut || thisE > 20) { 
+    thisStatus = kNone;
+    DBGL(3,Form(" %9f<%9f || %9f>20, 0'ed", thisE, lowCut, thisE));
+    if (prevStatus == kCandidate) prevStatus = kNone;
+    return 0;  
+  }
+  // It this strip was merged with the previous strip, then 
+  // make the next strip a candidate and zero the value in this strip. 
+  if (thisStatus == kMergedWithOther) { 
+    DBGL(3,Form(" Merged with other, 0'ed"));
+    return 0;
   }
 
-  if (delta > 100) {
-    AliWarning(Form("FMD%d%c, eta=%f, Delta=%f, xi=%f", d, r, eta, delta, xi));
+  // Get the upper cut on the sharing 
+  Double_t highCut = GetHighCut(d, r, eta ,false);
+  if (highCut <= 0) {
+    thisStatus = kNone;
+    return 0;
+  }
+
+  // Only for low-flux  events 
+  if (lowFlux) { 
+    // If this is smaller than the next, and the next is larger 
+    // then the cut, mark both as candidates for sharing. 
+    // They should be be merged on the next iteration 
+    if (thisE < nextE && nextE > highCut) { 
+      thisStatus = kCandidate;
+      nextStatus = kCandidate;
+      return 0;
+    }
+  }
+  
+  // Variable to sum signal in 
+  Double_t totalE = thisE;
+
+  // If the previous signal was between the two cuts, and is still
+  // marked as a candidate , then merge it in here.
+  if (prevStatus == kCandidate && prevE > lowCut && prevE < highCut) {
+    totalE     += prevE;
+    prevStatus =  kMergedWithOther;
+    DBG(3, Form(" Prev candidate %9f<%9f<%9f->%9f", 
+               lowCut, prevE, highCut, totalE));
+  }
+
+  // If the next signal is between the two cuts, then merge it here
+  if (nextE > lowCut && nextE < highCut) { 
+    totalE     += nextE;
+    nextStatus =  kMergedWithOther;
+    DBG(3, Form(" Next %9f<%9f<%9f->%9f", lowCut, nextE, highCut,totalE));
+  }
+  
+  // If the signal is bigger than the high-cut, then return the value 
+  if (totalE > highCut) { 
+    if (totalE > thisE) 
+      thisStatus = kMergedInto;
+    else 
+      thisStatus = kNone;
+    DBGL(3, Form(" %9f>%f9 (was %9f) -> %9f %s", 
+                totalE, highCut, thisE, totalE,status2String(thisStatus)));
+    return totalE;
+  }
+  else {
+    // This is below cut, so flag next as a candidate 
+    DBG(3, Form(" %9f<=%9f, next candidate", totalE, highCut));
+    nextStatus = kCandidate;
+  }
+  
+  // If the total signal is smaller than low cut then zero this and kill this 
+  if (totalE < lowCut)  {
+    DBGL(3, Form(" %9f<%9f (was %9f), zeroed", totalE, lowCut, thisE));
+    thisStatus = kNone;
+    return 0;
   }
-  Double_t highCut = (delta - fNXi * xi);
-  return highCut;
+
+  // If total signal not above high cut or lower than low cut, 
+  // mark this as a candidate for merging into the next, and zero signal 
+  DBGL(3, Form(" %9f<%9f<%9f (was %9f), zeroed, candidate", 
+                  lowCut, totalE, highCut, thisE));
+  thisStatus = kCandidate;
+  return 0;
 }
           
 //_____________________________________________________________________
@@ -361,7 +593,7 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
 
   // IF the multiplicity is very large, or below the cut, reject it, and 
   // flag it as candidate for sharing 
-  Double_t    lowCut = GetLowCut();
+  Double_t    lowCut = GetLowCut(d,r,eta);
   if(mult > 12 || mult < lowCut) {
     usedThis      = kFALSE;
     usedPrev      = kFALSE;
@@ -379,7 +611,12 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
   //analyse and perform sharing on one strip
   
   // Get the high cut 
-  Double_t highCut = GetHighCut(d, r, eta);
+  Double_t highCut = GetHighCut(d, r, eta, false);
+  if (highCut <= 0) {
+    usedThis = kFALSE;
+    usedPrev = kTRUE;
+    return 0;
+  }
 
   // If this signal is smaller than the next, and the next signal is smaller 
   // than then the high cut, and it's a low-flux event, then mark this and 
@@ -406,7 +643,7 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
     totalE += prevE;
 
   // If the next signal is larger than the low cut, and 
-  // the next signal is smaller than the low cut, then add the next signal 
+  // the next signal is smaller than the high cut, then add the next signal 
   // to this, and marked the next signal as used 
   if(nextE > lowCut && nextE < highCut ) {
     totalE      += nextE;
@@ -415,8 +652,8 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
 
   // If we're processing on non-angle corrected data, we should do the 
   // angle correction here 
-  if (!fCorrectAngles) 
-    totalE = AngleCorrect(totalE, eta);
+  // if (!fCorrectAngles) 
+  //   totalE = AngleCorrect(totalE, eta);
 
   // Fill post processing histogram
   // if(totalE > fLowCut)
@@ -494,8 +731,17 @@ AliFMDSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents)
 
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
-  while ((o = static_cast<RingHistos*>(next())))
+  THStack* sums = new THStack("sums", "Sum of ring signals");
+  while ((o = static_cast<RingHistos*>(next()))) {
     o->ScaleHistograms(d, nEvents);
+    TH1D* sum = o->fSum->ProjectionX(o->GetName(), 1, o->fSum->GetNbinsY(),"e");
+    sum->Scale(1., "width");
+    sum->SetTitle(o->GetName());
+    sum->SetDirectory(0);
+    sum->SetYTitle("#sum #Delta/#Delta_{mip}");
+    sums->Add(sum);
+  }
+  d->Add(sums);
 }
 
 //____________________________________________________________________
@@ -513,6 +759,26 @@ AliFMDSharingFilter::DefineOutput(TList* dir)
   TList* d = new TList;
   d->SetName(GetName());
   dir->Add(d);
+
+  fSummed = new TH2I("operations", "Strip operations", 
+                    kMergedInto, kNone-.5, kMergedInto+.5,
+                    51201, -.5, 51200.5);
+  fSummed->SetXTitle("Operation");
+  fSummed->SetYTitle("# of strips");
+  fSummed->SetZTitle("Events");
+  fSummed->GetXaxis()->SetBinLabel(kNone,            "None");
+  fSummed->GetXaxis()->SetBinLabel(kCandidate,       "Candidate");
+  fSummed->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
+  fSummed->GetXaxis()->SetBinLabel(kMergedInto,      "Merged into");
+  fSummed->SetDirectory(0);
+  d->Add(fSummed);
+
+  fHighCuts = new TH2D("highCuts", "High cuts used", 1,0,1, 1,0,1);
+  fHighCuts->SetXTitle("#eta");
+  fHighCuts->SetDirectory(0);
+  d->Add(fHighCuts);
+
+
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
   while ((o = static_cast<RingHistos*>(next()))) {
@@ -533,10 +799,13 @@ AliFMDSharingFilter::Print(Option_t* /*option*/) const
   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
   ind[gROOT->GetDirLevel()] = '\0';
   std::cout << ind << ClassName() << ": " << GetName() << '\n'
+           << std::boolalpha 
+           << ind << " Debug:                  " << fDebug << "\n"
            << ind << " Low cut:                " << fLowCut << '\n'
            << ind << " N xi factor:            " << fNXi    << '\n'
-           << ind << " Use corrected angles:   " 
-           << (fCorrectAngles ? "yes" : "no") << std::endl;
+           << ind << " Include sigma in cut:   " << fIncludeSigma << '\n'
+           << ind << " Use corrected angles:   " << fCorrectAngles 
+           << std::noboolalpha << std::endl;
 }
   
 //====================================================================
@@ -544,6 +813,10 @@ AliFMDSharingFilter::RingHistos::RingHistos()
   : AliForwardUtil::RingHistos(), 
     fBefore(0), 
     fAfter(0), 
+    fBeforeAfter(0),
+    fNeighborsBefore(0),
+    fNeighborsAfter(0),
+    fSum(0),
     fHits(0),
     fNHits(0)
 {
@@ -558,6 +831,10 @@ AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
   : AliForwardUtil::RingHistos(d,r), 
     fBefore(0), 
     fAfter(0),
+    fBeforeAfter(0),
+    fNeighborsBefore(0),
+    fNeighborsAfter(0),
+    fSum(0),
     fHits(0),
     fNHits(0)
 {
@@ -568,28 +845,51 @@ AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
   //    d detector
   //    r ring 
   //
-  fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()), 
-                    Form("Energy loss in %s (reconstruction)", fName.Data()), 
+  fBefore = new TH1D("esdEloss", "Energy loss (reconstruction)", 
                     1000, 0, 25);
-  fAfter  = new TH1D(Form("%s_Ana_Eloss", fName.Data()), 
-                    Form("Energy loss in %s (sharing corrected)",
-                         fName.Data()), 1000, 0, 25);
   fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
   fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
-  fBefore->SetFillColor(kRed+1);
+  fBefore->SetFillColor(Color());
   fBefore->SetFillStyle(3001);
-  // fBefore->Sumw2();
   fBefore->SetDirectory(0);
-  fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
-  fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
-  fAfter->SetFillColor(kBlue+1);
-  fAfter->SetFillStyle(3001);
-  // fAfter->Sumw2();
+
+  fAfter  = static_cast<TH1D*>(fBefore->Clone("anaEloss"));
+  fAfter->SetTitle("Energy loss in %s (sharing corrected)");
+  fAfter->SetFillColor(Color()+2);
   fAfter->SetDirectory(0);
 
-  fHits = new TH1D(Form("%s_Hits", fName.Data()), 
-                  Form("Number of hits in %s", fName.Data()), 
-                  200, 0, 200000);
+  Double_t max = 15;
+  Double_t min = -1;
+  Int_t    n   = int((max-min) / (max / 300));
+  fBeforeAfter = new TH2D("beforeAfter", "Before and after correlation", 
+                         n, min, max, n, min, max);
+  fBeforeAfter->SetXTitle("#Delta E/#Delta E_{mip} before");
+  fBeforeAfter->SetYTitle("#Delta E/#Delta E_{mip} after");
+  fBeforeAfter->SetZTitle("Correlation");
+  fBeforeAfter->SetDirectory(0);
+
+  fNeighborsBefore = static_cast<TH2D*>(fBeforeAfter->Clone("neighborsBefore"));
+  fNeighborsBefore->SetTitle("Correlation of neighbors befire");
+  fNeighborsBefore->SetXTitle("#Delta E_{i}/#Delta E_{mip}");
+  fNeighborsBefore->SetYTitle("#Delta E_{i+1}/#Delta E_{mip}");
+  fNeighborsBefore->SetDirectory(0);
+  
+  fNeighborsAfter = 
+    static_cast<TH2D*>(fNeighborsBefore->Clone("neighborsAfter"));
+  fNeighborsAfter->SetTitle("Correlation of neighbors after");
+  fNeighborsAfter->SetDirectory(0);
+
+  fSum = new TH2D("summed", "Summed signal", 200, -4, 6, 
+                 (fRing == 'I' || fRing == 'i' ? 20 : 40), 0, 2*TMath::Pi());
+  fSum->SetDirectory(0);
+  fSum->Sumw2();
+  fSum->SetMarkerColor(Color());
+  // fSum->SetFillColor(Color());
+  fSum->SetXTitle("#eta");
+  fSum->SetYTitle("#varphi [radians]");
+  fSum->SetZTitle("#sum #Delta/#Delta_{mip}(#eta,#varphi) ");
+  
+  fHits = new TH1D("hits", "Number of hits", 200, 0, 200000);
   fHits->SetDirectory(0);
   fHits->SetXTitle("# of hits");
   fHits->SetFillColor(kGreen+1);
@@ -599,6 +899,10 @@ AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
   : AliForwardUtil::RingHistos(o), 
     fBefore(o.fBefore), 
     fAfter(o.fAfter),
+    fBeforeAfter(o.fBeforeAfter),
+    fNeighborsBefore(o.fNeighborsBefore),
+    fNeighborsAfter(o.fNeighborsAfter),
+    fSum(o.fSum),
     fHits(o.fHits),
     fNHits(o.fNHits)
 {
@@ -631,10 +935,14 @@ AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
   if (fAfter)  delete  fAfter;
   if (fHits)   delete fHits;
   
-  fBefore = static_cast<TH1D*>(o.fBefore->Clone());
-  fAfter  = static_cast<TH1D*>(o.fAfter->Clone());
-  fHits  = static_cast<TH1D*>(o.fHits->Clone());
-  
+  fBefore          = static_cast<TH1D*>(o.fBefore->Clone());
+  fAfter           = static_cast<TH1D*>(o.fAfter->Clone());
+  fBeforeAfter     = static_cast<TH2D*>(o.fBeforeAfter->Clone());
+  fNeighborsBefore = static_cast<TH2D*>(o.fNeighborsBefore->Clone());
+  fNeighborsAfter  = static_cast<TH2D*>(o.fNeighborsAfter->Clone());
+  fHits            = static_cast<TH1D*>(o.fHits->Clone());
+  fSum             = static_cast<TH2D*>(o.fSum->Clone());
+
   return *this;
 }
 //____________________________________________________________________
@@ -643,9 +951,6 @@ AliFMDSharingFilter::RingHistos::~RingHistos()
   // 
   // Destructor 
   //
-  if (fBefore) delete fBefore;
-  if (fAfter)  delete fAfter;
-  if (fHits)   delete fHits;
 }
 //____________________________________________________________________
 void
@@ -673,12 +978,13 @@ AliFMDSharingFilter::RingHistos::ScaleHistograms(const TList* dir,
   TList* l = GetOutputList(dir);
   if (!l) return; 
 
-  TH1D* before = static_cast<TH1D*>(l->FindObject(Form("%s_ESD_ELoss",
-                                                      fName.Data())));
-  TH1D* after  = static_cast<TH1D*>(l->FindObject(Form("%s_Ana_ELoss",
-                                                      fName.Data())));
+  TH1D* before = static_cast<TH1D*>(l->FindObject("esdELoss"));
+  TH1D* after  = static_cast<TH1D*>(l->FindObject("anaELoss"));
+  TH2D* summed = static_cast<TH2D*>(l->FindObject("summed"));
   if (before) before->Scale(1./nEvents);
   if (after)  after->Scale(1./nEvents);
+  if (summed) summed->Scale(1./nEvents);
+  
 }
 
 //____________________________________________________________________
@@ -692,9 +998,15 @@ AliFMDSharingFilter::RingHistos::Output(TList* dir)
   //    dir where to store 
   //
   TList* d = DefineOutputList(dir);
+
   d->Add(fBefore);
   d->Add(fAfter);
+  d->Add(fBeforeAfter);
+  d->Add(fNeighborsBefore);
+  d->Add(fNeighborsAfter);
   d->Add(fHits);
+  d->Add(fSum);
+
   dir->Add(d);
 }
 
index f103bc810354c2505c0b7431dd921268052055e2..14bfe90ad919a4793009fe9271573e374c6b6cc3 100644 (file)
@@ -24,6 +24,7 @@ class AliESDFMD;
 class TAxis;
 class TList;
 class TH2;
+class AliFMDFloatMap;
 
 /**
  * Class to do the sharing correction.  That is, a filter that merges 
@@ -53,6 +54,12 @@ class TH2;
 class AliFMDSharingFilter : public TNamed
 {
 public: 
+  enum Status { 
+    kNone             = 1, 
+    kCandidate        = 2, 
+    kMergedWithOther  = 3, 
+    kMergedInto       = 4
+  };
   /** 
    * Destructor
    */
@@ -82,6 +89,11 @@ public:
    */
   AliFMDSharingFilter& operator=(const AliFMDSharingFilter& o);
 
+  /** 
+   * Initialize 
+   * 
+   */
+  void Init();
   /** 
    * Set the low cut used for sharing 
    * 
@@ -107,14 +119,21 @@ public:
    * otherwise use de-corrected signals.  In the final output, the 
    * signals are always angle corrected. 
    */
-  void UseAngleCorrectedSignals(Bool_t use) { fCorrectAngles = use; }
+  void SetUseAngleCorrectedSignals(Bool_t use) { fCorrectAngles = use; }
   /** 
    * Set the number of landau width to subtract from the most probably
    * value to get the high cut for the merging algorithm.
    * 
    * @param n Number of @f$ \xi@f$ 
    */
-  void SetNXi(Short_t n) { fNXi = n; }
+  void SetNXi(Double_t n) { fNXi = n; }
+  /** 
+   * Whether to include sigma in the number subtracted from the MPV to
+   * get the high cut
+   * 
+   * @param u If true, then high cut is @f$ \Delta_{mp} - n(\xi+\sigma)@f$ 
+   */
+  void SetIncludeSigma(Bool_t u) { fIncludeSigma = u; }
   /** 
    * Filter the input AliESDFMD object
    * 
@@ -213,6 +232,10 @@ protected:
     void ScaleHistograms(const TList* dir, Int_t nEvents);
     TH1D*     fBefore;       // Distribution of signals before filter
     TH1D*     fAfter;        // Distribution of signals after filter
+    TH2D*     fBeforeAfter;  // Correlation of before and after 
+    TH2D*     fNeighborsBefore; // Correlation of neighbors 
+    TH2D*     fNeighborsAfter; // Correlation of neighbors 
+    TH2D*     fSum;          // Summed signal 
     TH1D*     fHits;         // Distribution of hit strips. 
     Int_t     fNHits;        // Number of hit strips per event
     ClassDef(RingHistos,1);
@@ -270,6 +293,18 @@ protected:
                               UShort_t t,
                               Bool_t&  usedPrev, 
                               Bool_t&  usedThis) const;
+  Double_t MultiplicityOfStrip(Double_t thisE,
+                              Double_t prevE,
+                              Double_t nextE,
+                              Double_t eta,
+                              Bool_t   lowFlux,
+                              UShort_t d,
+                              Char_t   r,
+                              UShort_t s,
+                              UShort_t t,
+                              Status&  prevStatus, 
+                              Status&  thisStatus, 
+                              Status&  nextStatus) const;
   /** 
    * Angle correct the signal 
    * 
@@ -288,27 +323,45 @@ protected:
    * @return Angle un-corrected signal 
    */
   Double_t DeAngleCorrect(Double_t mult, Double_t eta) const;
-  /**
+  /** 
    * Get the high cut.  The high cut is defined as the 
    * most-probably-value peak found from the energy distributions, minus 
    * 2 times the width of the corresponding Landau.
+   * 
+   * @param d      Detector 
+   * @param r      Ring 
+   * @param eta    Eta value 
+   * @param errors If false, do not show errors 
+   * 
+   * @return 0 or less on failure, otherwise @f$\Delta_{mp}-n\xi@f$ 
    */
-  virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta) const;
+  virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta,
+                             Bool_t errors=true) const;
   /**
    * Get the low cut.  Normally, the low cut is taken to be the lower
    * value of the fit range used when generating the energy loss fits.
    * However, if fLowCut is set (using SetLowCit) to a value greater
    * than 0, then that value is used.
+   *
+   * @param d    Detector
+   * @param r    Ring 
+   * @param eta  Eta value 
+   * 
+   * @return 
    */
-  virtual Double_t GetLowCut() const;
+  virtual Double_t GetLowCut(UShort_t d, Char_t r, Double_t eta) const;
 
   TList    fRingHistos;    // List of histogram containers
   Double_t fLowCut;        // Low cut on sharing
   Bool_t   fCorrectAngles; // Whether to work on angle corrected signals
-  Short_t  fNXi;           // Number of xi's from Delta to stop merging
+  Double_t fNXi;           // Number of xi's from Delta to stop merging
+  Bool_t   fIncludeSigma;  // Whether to include sigma in cut 
+  TH2*     fSummed;        // Operations histogram 
+  TH2*     fHighCuts;      // High cuts used
+  AliFMDFloatMap* fOper;   // Operation done per strip 
   Int_t    fDebug;         // Debug level 
 
-  ClassDef(AliFMDSharingFilter,1); //
+  ClassDef(AliFMDSharingFilter,2); //
 };
 
 #endif
index 1fda2c5c4ba91329eb75bc29f2a6afcb0824e820..5376027afc06f4d6476c0d37259789a3e70ffa69 100644 (file)
@@ -23,7 +23,8 @@
 #include "AliInputEventHandler.h"
 #include "AliStack.h"
 #include "AliMCEvent.h"
-//#include "AliFMDStripIndex.h"
+#include "AliAODForwardMult.h"
+#include "AliFMDStripIndex.h"
 #include <TH1.h>
 #include <TH3D.h>
 #include <TH2D.h>
@@ -56,10 +57,11 @@ namespace {
 //====================================================================
 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
   : AliAnalysisTaskSE(),
+    fInspector(),
+    fFirstEvent(true),
     fHEvents(0), 
     fHEventsTr(0), 
     fHEventsTrVtx(0),
-    fHEventsVtx(0), 
     fHTriggers(0),
     fPrimaryInnerAll(0),   
     fPrimaryOuterAll(0),   
@@ -89,11 +91,12 @@ AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
 
 //____________________________________________________________________
 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
-  : AliAnalysisTaskSE(name), 
+  : AliAnalysisTaskSE(name),
+    fInspector("eventInspector"), 
+    fFirstEvent(true),
     fHEvents(0), 
     fHEventsTr(0), 
     fHEventsTrVtx(0),
-    fHEventsVtx(0), 
     fHTriggers(0),
     fPrimaryInnerAll(0),   
     fPrimaryOuterAll(0),   
@@ -120,15 +123,17 @@ AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
   //    name Name of task 
   //
   DefineOutput(1, TList::Class());
+  DefineOutput(2, TList::Class());
 }
 
 //____________________________________________________________________
 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o)
   : AliAnalysisTaskSE(o),
+    fInspector(o.fInspector),
+    fFirstEvent(o.fFirstEvent),
     fHEvents(o.fHEvents), 
     fHEventsTr(o.fHEventsTr), 
     fHEventsTrVtx(o.fHEventsTrVtx),
-    fHEventsVtx(o.fHEventsVtx), 
     fHTriggers(o.fHTriggers),
     fPrimaryInnerAll(o.fPrimaryInnerAll),   
     fPrimaryOuterAll(o.fPrimaryOuterAll),   
@@ -169,6 +174,8 @@ AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o)
   // Return:
   //    Reference to this object 
   //
+  fInspector         = o.fInspector;
+  fFirstEvent        = o.fFirstEvent;
   fHEventsTr         = o.fHEventsTr;
   fHEventsTrVtx      = o.fHEventsTrVtx;
   fHTriggers         = o.fHTriggers;
@@ -333,7 +340,8 @@ AliForwardMCCorrectionsTask::UserCreateOutputObjects()
   // 
   //
   fList = new TList;
-  fList->SetName(GetName());
+  fList->SetOwner();
+  fList->SetName(Form("%sSums", GetName()));
 
   fHEvents = new TH1I(GetEventName(false,false),
                      "Number of all events", 
@@ -371,19 +379,6 @@ AliForwardMCCorrectionsTask::UserCreateOutputObjects()
   fHEventsTrVtx->SetDirectory(0);
   fList->Add(fHEventsTrVtx);
   
-  fHEventsVtx = new TH1I(GetEventName(false,true),
-                        "Number of events w/vertex", 
-                        fVtxAxis.GetNbins(), 
-                        fVtxAxis.GetXmin(), 
-                        fVtxAxis.GetXmax());
-  fHEventsVtx->SetXTitle("v_{z} [cm]");
-  fHEventsVtx->SetYTitle("# of events");
-  fHEventsVtx->SetFillColor(kBlue+1);
-  fHEventsVtx->SetFillStyle(3001);
-  fHEventsVtx->SetDirectory(0);
-  fList->Add(fHEventsVtx);
-
-      
   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
   fHTriggers->SetFillColor(kRed+1);
   fHTriggers->SetFillStyle(3001);
@@ -444,6 +439,9 @@ AliForwardMCCorrectionsTask::UserCreateOutputObjects()
   strips->Add(fStripsFMD3o);
   fList->Add(strips);
 
+  fInspector.Init(fVtxAxis);
+  fInspector.DefineOutput(fList);
+
   PostData(1, fList);
 }
 //____________________________________________________________________
@@ -470,52 +468,57 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*)
     AliWarning("No ESD event found for input event");
     return;
   }
-  
-  // Get the particle stack 
-  AliStack* stack = mcEvent->Stack();
 
-  // Get the event generate header 
-  AliHeader*          mcHeader  = mcEvent->Header();
-  AliGenEventHeader*  genHeader = mcHeader->GenEventHeader();
-  
-  // Get the generator vertex 
-  TArrayF mcVertex;
-  genHeader->PrimaryVertex(mcVertex);
-  Double_t mcVtxZ = mcVertex.At(2);
-
-  // Check the MC vertex is in range 
-  Int_t mcVtxBin = fVtxAxis.FindBin(mcVtxZ);
-  if (mcVtxBin < 1 || mcVtxBin > fVtxAxis.GetNbins()) {
-#ifdef VERBOSE 
-    AliWarning(Form("MC v_z=%f is out of rante [%,%f]", 
-                   mcVtxBin, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
-#endif
-    return;
+  //--- Read run information -----------------------------------------
+  if (fFirstEvent && esd->GetESDRun()) {
+    fInspector.ReadRunDetails(esd);
+    
+    AliInfo(Form("Initializing with parameters from the ESD:\n"
+                "         AliESDEvent::GetBeamEnergy()   ->%f\n"
+                "         AliESDEvent::GetBeamType()     ->%s\n"
+                "         AliESDEvent::GetCurrentL3()    ->%f\n"
+                "         AliESDEvent::GetMagneticField()->%f\n"
+                "         AliESDEvent::GetRunNumber()    ->%d\n",
+                esd->GetBeamEnergy(),
+                esd->GetBeamType(),
+                esd->GetCurrentL3(),
+                esd->GetMagneticField(),
+                esd->GetRunNumber()));
+    
+    fFirstEvent = false;
   }
 
-  // UInt_t   triggers;
-  // Bool_t   gotTrigggers = false;
-  Bool_t   gotInel = false;
-  // Double_t vZ;
-  Bool_t   gotVertex = false;
-#if 0
-  // Use event inspector instead 
-  // Get the triggers 
-  UInt_t triggers = 0;
-  Bool_t gotTriggers = AliForwardUtil::ReadTriggers(esd,triggers,fHTriggers);
-  Bool_t gotInel     = triggers & AliAODForwardMult::kInel;
-  
-  // Get the ESD vertex 
-  Double_t vZ = -1000000;
-  Bool_t gotVertex = AliForwardUtil::ReadVertex(esd,vZ);
-#endif
-
+  // Get the particle stack 
+  AliStack* stack = mcEvent->Stack();
 
+  // Some variables 
+  UInt_t   triggers; // Trigger bits
+  Bool_t   lowFlux;  // Low flux flag
+  UShort_t iVz;      // Vertex bin from ESD
+  Double_t vZ;       // Z coordinate from ESD
+  Double_t cent;     // Centrality 
+  UShort_t iVzMc;    // Vertex bin from MC
+  Double_t vZMc;     // Z coordinate of IP vertex from MC
+  Double_t b;        // Impact parameter
+  Int_t    nPart;    // Number of participants 
+  Int_t    nBin;     // Number of binary collisions 
+  Double_t phiR;     // Reaction plane from MC
+  UShort_t nClusters;// Number of SPD clusters 
+  // Process the data 
+  UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, 
+                                    cent, nClusters);
+  /*UInt_t retMC  =*/ fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc, 
+                                          b, nPart, nBin, phiR);
+
+  // Bool_t isInelMC = true; // (triggers & AliAODForwardMult::kB);
+  // Bool_t isNSDMC  = (triggers & AliAODForwardMult::kMCNSD);
+  Bool_t isInel   = triggers & AliAODForwardMult::kInel;
+  // Bool_t isNSD    = triggers & AliAODForwardMult::kNSD;
+  Bool_t hasVtx   = retESD == AliFMDMCEventInspector::kOk;
   // Fill the event count histograms 
-  if (gotInel)              fHEventsTr->Fill(mcVtxZ);
-  if (gotInel && gotVertex) fHEventsTrVtx->Fill(mcVtxZ);
-  if (gotVertex)            fHEventsVtx->Fill(mcVtxZ);
-  fHEvents->Fill(mcVtxZ);
+  if (isInel)           fHEventsTr->Fill(vZMc);
+  if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
+  fHEvents->Fill(vZMc);
 
   // Cache of the hits 
   AliFMDFloatMap hitsByStrip;
@@ -539,10 +542,10 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*)
     Double_t phi = particle->Phi();
     
     if (isCharged && isPrimary) 
-      FillPrimary(gotInel, gotVertex, mcVtxZ, eta, phi);
+      FillPrimary(isInel, hasVtx, vZMc, eta, phi);
     
     // For the rest, ignore non-collisions, and events out of vtx range 
-    if (!gotInel || gotVertex) continue;
+    if (!isInel || !hasVtx) continue;
     
     Int_t nTrRef = particle->GetNumberOfTrackReferences();
     for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
@@ -558,7 +561,7 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*)
       // Get the detector coordinates 
       UShort_t d = 0, s = 0, t = 0;
       Char_t r = '\0';
-      // AliFMDStripIndex::Unpack(trackRef->UserId(), d, r, s, t);
+      AliFMDStripIndex::Unpack(trackRef->UserId(), d, r, s, t);
       
       // Check if mother (?) is charged and that we haven't dealt with it 
       // already
@@ -572,7 +575,7 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*)
       // Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t);
 
       // Fill strip information into histograms 
-      FillStrip(d, r, mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1);
+      FillStrip(d, r, vZMc, eta, phi, hitsByStrip(d,r,s,t) == 1);
 
       // Set the last processed track number - marking it as done for
       // this strip
@@ -588,20 +591,21 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*)
 
 //____________________________________________________________________
 void
-AliForwardMCCorrectionsTask::FillPrimary(Bool_t gotInel, Bool_t gotVtx, 
-                                   Double_t vz, Double_t eta, Double_t phi) 
+AliForwardMCCorrectionsTask::FillPrimary(Bool_t isInel, Bool_t hasVtx, 
+                                        Double_t vz, 
+                                        Double_t eta, Double_t phi) 
 {
   // 
   // Fill in primary information
   // 
   // Parameters:
-  //    gotInel   Got INEL trigger from ESD
+  //    isInel   Got INEL trigger from ESD
   //    gotVtx    Got vertex Z from ESD 
   //    vz        @f$z@f$ coordinate of interation point
   //    eta       Pseudo rapidity 
   //    phi       Azimuthal angle
   //
-  if (gotInel && gotVtx) {
+  if (isInel && hasVtx) {
     // This takes the place of hPrimary_FMD_<r>_vtx<v> and 
     // Analysed_FMD<r>_vtx<v>
     fPrimaryInnerTrVtx->Fill(vz,eta,phi);
@@ -698,32 +702,32 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*)
   //    option Not used 
   //
 
-  TList* list = dynamic_cast<TList*>(GetOutputData(1));
-  if (!list) {
+  fList = dynamic_cast<TList*>(GetOutputData(1));
+  if (!fList) {
     AliError("No output list defined");
     return;
   }
 
-  TList* primaries = static_cast<TList*>(list->FindObject("primaries"));
-  TList* hits      = static_cast<TList*>(list->FindObject("hits"));
-  TList* strips    = static_cast<TList*>(list->FindObject("strips"));
+  TList* primaries = static_cast<TList*>(fList->FindObject("primaries"));
+  TList* hits      = static_cast<TList*>(fList->FindObject("hits"));
+  TList* strips    = static_cast<TList*>(fList->FindObject("strips"));
   if (!primaries || !hits || !strips) { 
     AliError(Form("Could not find all sub-lists in %s (p=%p,h=%p,s=%p)",
-                 list->GetName(), primaries, hits, strips));
+                 fList->GetName(), primaries, hits, strips));
     return;
   }
 
   TH1I* eventsAll = 
-    static_cast<TH1I*>(list->FindObject(GetEventName(false,false)));
+    static_cast<TH1I*>(fList->FindObject(GetEventName(false,false)));
   TH1I* eventsTr = 
-    static_cast<TH1I*>(list->FindObject(GetEventName(true,false)));
+    static_cast<TH1I*>(fList->FindObject(GetEventName(true,false)));
   TH1I* eventsVtx = 
-    static_cast<TH1I*>(list->FindObject(GetEventName(false,true)));
+    static_cast<TH1I*>(fList->FindObject(GetEventName(false,true)));
   TH1I* eventsTrVtx = 
-    static_cast<TH1I*>(list->FindObject(GetEventName(true,true)));
+    static_cast<TH1I*>(fList->FindObject(GetEventName(true,true)));
   if (!eventsAll || !eventsTr || !eventsVtx || !eventsTrVtx) {
     AliError(Form("Could not find all event histograms in %s "
-                 "(a=%p,t=%p,v=%p,tv=%p)", list->GetName(), 
+                 "(a=%p,t=%p,v=%p,tv=%p)", fList->GetName(), 
                  eventsAll, eventsTr, eventsVtx, eventsTrVtx));
     return;
   }
@@ -738,7 +742,7 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*)
     static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('O',true)));
   if (!primIAll || !primOAll || !primITrVtx || !primOTrVtx) {
     AliError(Form("Could not find all primary particle histograms in %s "
-                 "(ai=%p,ao=%p,tvi=%p,tvo=%p)", list->GetName(), 
+                 "(ai=%p,ao=%p,tvi=%p,tvo=%p)", fList->GetName(), 
                  primIAll, primOAll, primITrVtx, primOTrVtx));
     return;
   }
@@ -767,6 +771,11 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*)
     return;
   }
 
+  // Output list 
+  TList* output = new TList;
+  output->SetOwner();
+  output->SetName(Form("%sResults", GetName()));
+
   // Calculate the over-all event selection efficiency 
   TH1D* selEff = new TH1D("selEff", "Event selection efficiency", 
                          fVtxAxis.GetNbins(), 
@@ -778,14 +787,15 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*)
   selEff->SetFillStyle(3001);
   selEff->Add(eventsAll);
   selEff->Divide(eventsTrVtx);
-  list->Add(selEff);
+  output->Add(selEff);
 
   // Loop over vertex bins and do vertex dependendt stuff 
   for (Int_t v = 1; v <= fVtxAxis.GetNbins(); v++) {
-    // Make a sub-list in the output 
+    // Make a sub-fList in the output 
     TList* vl = new TList;
+    vl->SetOwner();
     vl->SetName(Form("vtx%02d", v));
-    list->Add(vl);
+    output->Add(vl);
 
     // Get event counts 
     Int_t nEventsAll   = Int_t(eventsAll->GetBinContent(v));
@@ -858,10 +868,11 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*)
        doubleHit->SetFillStyle(3001);
        doubleHit->Sumw2();
        doubleHit->Divide(hStrips);
-       list->Add(doubleHit);
-      }
-    }
-  }
+       output->Add(doubleHit);
+      } // for q 
+    } // for d 
+  } // for v 
+  PostData(2, output);
 }
 
 //____________________________________________________________________
index 80b7f09b4d5b5aacf29660c960cac88e468ec844..7701145b7784a392fe62ed304afd9fd2c42fa8d3 100644 (file)
@@ -15,8 +15,8 @@
  */
 #include <AliAnalysisTaskSE.h>
 #include <AliESDFMD.h>
+#include "AliFMDMCEventInspector.h"
 #include <TH1I.h>
-class AliFMDAnaParameters;
 class AliESDEvent;
 class TH2D;
 class TH3D;
@@ -186,10 +186,11 @@ protected:
   void FillStrip(UShort_t d, Char_t r, 
                 Double_t vz, Double_t eta, Double_t phi,
                 Bool_t first);
+  AliFMDMCEventInspector fInspector; // Event inspector 
+  Bool_t fFirstEvent;        // First event flag 
   TH1I*  fHEvents;           // All Events
   TH1I*  fHEventsTr;         // Histogram of events w/trigger
   TH1I*  fHEventsTrVtx;      // Events w/trigger and vertex 
-  TH1I*  fHEventsVtx;        // Events w/vertex 
   TH1I*  fHTriggers;         // Triggers
   TH3D*  fPrimaryInnerAll;   // Distribution of primaries - all events
   TH3D*  fPrimaryOuterAll;   // Distribution of primaries - all events
index 4e29b526770eaa7750c923a1e510475e88bff4c0..5b5d77506b323499d145121643f418d9bd980a43 100644 (file)
@@ -38,6 +38,8 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
     fMCESDFMD(),
     fMCHistos(),
     fMCAODFMD(),
+    fRingSums(),
+    fMCRingSums(),
     fPrimary(0),
     fEventInspector(),
     fSharingFilter(),
@@ -61,6 +63,8 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
     fMCESDFMD(),
     fMCHistos(),
     fMCAODFMD(kTRUE),
+    fRingSums(),
+    fMCRingSums(),
     fPrimary(0),
     fEventInspector("event"),
     fSharingFilter("sharing"), 
@@ -88,6 +92,8 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMul
     fMCESDFMD(o.fMCESDFMD),
     fMCHistos(o.fMCHistos),
     fMCAODFMD(o.fMCAODFMD),
+    fRingSums(o.fRingSums),
+    fMCRingSums(o.fMCRingSums),
     fPrimary(o.fPrimary),
     fEventInspector(o.fEventInspector),
     fSharingFilter(o.fSharingFilter),
@@ -130,6 +136,8 @@ AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
   fAODFMD            = o.fAODFMD;
   fMCHistos          = o.fMCHistos;
   fMCAODFMD          = o.fMCAODFMD;
+  fRingSums          = o.fRingSums;
+  fMCRingSums        = o.fMCRingSums;
   fPrimary           = o.fPrimary;
   fList              = o.fList;
 
@@ -152,6 +160,14 @@ AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
   fCorrections.SetDebug(dbg);
   fHistCollector.SetDebug(dbg);
 }
+//____________________________________________________________________
+void
+AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
+{
+  fSharingFilter.SetOnlyPrimary(use);
+  fCorrections.SetSecondaryForMC(!use);
+}
+
 //____________________________________________________________________
 void
 AliForwardMCMultiplicityTask::InitializeSubs()
@@ -169,6 +185,8 @@ AliForwardMCMultiplicityTask::InitializeSubs()
   fAODFMD.Init(*pe);
   fMCHistos.Init(*pe);
   fMCAODFMD.Init(*pe);
+  fRingSums.Init(*pe);
+  fMCRingSums.Init(*pe);
 
   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
   fHData->SetStats(0);
@@ -176,7 +194,42 @@ AliForwardMCMultiplicityTask::InitializeSubs()
 
   fList->Add(fHData);
 
+  TList* rings = new TList;
+  rings->SetName("ringSums");
+  rings->SetOwner();
+  fList->Add(rings);
+
+  rings->Add(fRingSums.Get(1, 'I'));
+  rings->Add(fRingSums.Get(2, 'I'));
+  rings->Add(fRingSums.Get(2, 'O'));
+  rings->Add(fRingSums.Get(3, 'I'));
+  rings->Add(fRingSums.Get(3, 'O'));
+  fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
+  fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
+  fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
+  fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
+  fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
+
+  TList* mcRings = new TList;
+  mcRings->SetName("mcRingSums");
+  mcRings->SetOwner();
+  fList->Add(mcRings);
+
+  mcRings->Add(fMCRingSums.Get(1, 'I'));
+  mcRings->Add(fMCRingSums.Get(2, 'I'));
+  mcRings->Add(fMCRingSums.Get(2, 'O'));
+  mcRings->Add(fMCRingSums.Get(3, 'I'));
+  mcRings->Add(fMCRingSums.Get(3, 'O'));
+  fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
+  fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
+  fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
+  fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
+  fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
+
+
+
   fEventInspector.Init(*pv);
+  fSharingFilter.Init();
   fDensityCalculator.Init(*pe);
   fCorrections.Init(*pe);
   fHistCollector.Init(*pv,*pe);
@@ -249,13 +302,14 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   fMCAODFMD.Clear();
   fPrimary->Reset();
 
-  Bool_t   lowFlux  = kFALSE;
-  UInt_t   triggers = 0;
-  UShort_t ivz      = 0;
-  Double_t vz       = 0;
-  Double_t cent     = 0;
-  UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, 
-                                             ivz, vz, cent);
+  Bool_t   lowFlux   = kFALSE;
+  UInt_t   triggers  = 0;
+  UShort_t ivz       = 0;
+  Double_t vz        = 0;
+  Double_t cent      = 0;
+  UShort_t nClusters = 0;
+  UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
+                                              ivz, vz, cent, nClusters);
   UShort_t ivzMC    = 0;
   Double_t vzMC     = 0;
   Double_t phiR     = 0;
@@ -281,11 +335,13 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   fAODFMD.SetSNN(fEventInspector.GetEnergy());
   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
   fAODFMD.SetCentrality(cent);
+  fAODFMD.SetNClusters(nClusters);
 
   fMCAODFMD.SetTriggerBits(triggers);
   fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
   fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
   fMCAODFMD.SetCentrality(cent);
+  fMCAODFMD.SetNClusters(nClusters);
   
   //All events should be stored - HHD
   //if (isAccepted) MarkEventForStore();
@@ -346,11 +402,13 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   }
   fCorrections.CompareResults(fHistos, fMCHistos);
     
-  if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
+  if (!fHistCollector.Collect(fHistos, fRingSums, 
+                             ivz, fAODFMD.GetHistogram())) {
     AliWarning("Histogram collector failed");
     return;
   }
-  if (!fHistCollector.Collect(fMCHistos, ivz, fMCAODFMD.GetHistogram())) {
+  if (!fHistCollector.Collect(fMCHistos, fMCRingSums, 
+                             ivz, fMCAODFMD.GetHistogram())) {
     AliWarning("MC Histogram collector failed");
     return;
   }
@@ -411,6 +469,9 @@ AliForwardMCMultiplicityTask::Terminate(Option_t*)
   list->Add(dNdeta);
   list->Add(norm);
 
+  MakeRingdNdeta(list, "ringSums", list, "ringResults");
+  MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);
+
   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
   fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
index 96c9cd910cae5fdc50f4fe9c6836a0d81a068e7b..311f2d853942fdb4e8e2a2f247bbdf0481313459 100644 (file)
@@ -99,6 +99,12 @@ public:
   /** 
    * @} 
    */
+  /** 
+   * Process only primary MC tracks 
+   * 
+   * @param use 
+   */
+  void SetOnlyPrimary(Bool_t use);
   /** 
    * @{ 
    * @name Access to sub-algorithms 
@@ -186,6 +192,8 @@ protected:
   AliESDFMD              fMCESDFMD;     // MC 'Sharing corrected' ESD object
   AliForwardUtil::Histos fMCHistos;     // MC Cache histograms 
   AliAODForwardMult      fMCAODFMD;     // MC Output object
+  AliForwardUtil::Histos fRingSums;     // Cache histograms 
+  AliForwardUtil::Histos fMCRingSums;   // Cache histograms 
   TH2D*                  fPrimary;      // Per event primary particles 
 
   AliFMDMCEventInspector    fEventInspector;    // Algorithm
index e238ac832421f87839232076345db78fee8283be..7230dc5d588880a1e1ca70a156ff962798ba8cfa 100644 (file)
@@ -27,6 +27,7 @@
 #include "AliESDEvent.h"
 #include <TROOT.h>
 #include <TAxis.h>
+#include <THStack.h>
 #include <iostream>
 #include <iomanip>
 
@@ -207,6 +208,85 @@ AliForwardMultiplicityBase::MarkEventForStore() const
   ah->SetFillAOD(kTRUE);
 }
 
+//____________________________________________________________________
+void
+AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input, 
+                                          const char*  inName,
+                                          TList*       output,
+                                          const char*  outName,
+                                          Int_t        style) const
+{
+  // Make dN/deta for each ring found in the input list.  
+  // 
+  // A stack of all the dN/deta is also made for easy drawing. 
+  // 
+  // Note, that the distributions are normalised to the number of
+  // observed events only - they should be corrected for 
+  if (!input) return;
+  TList* list = static_cast<TList*>(input->FindObject(inName));
+  if (!list) { 
+    AliWarning(Form("No list %s found in %s", inName, input->GetName()));
+    return;
+  }
+  
+  TList* out = new TList;
+  out->SetName(outName);
+  out->SetOwner();
+  output->Add(out);
+
+  THStack*     dndetaRings = new THStack("all", "dN/d#eta per ring");
+  const char*  names[]     = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
+  const char** ptr         = names;
+  
+  while (*ptr) { 
+    TList* thisList = new TList;
+    thisList->SetOwner();
+    thisList->SetName(*ptr);
+    out->Add(thisList);
+
+    TH2D* h = static_cast<TH2D*>(list->FindObject(Form("%s_cache", *ptr)));
+    if (!h) { 
+      AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
+      ptr++;
+      continue;
+    }
+    TH2D* copy = static_cast<TH2D*>(h->Clone("sum"));
+    copy->SetDirectory(0);
+    thisList->Add(copy);
+    
+    TH1D* norm =static_cast<TH1D*>(h->ProjectionX("norm", 0, 0, ""));
+    for (Int_t i = 1; i <= copy->GetNbinsX(); i++) { 
+      for (Int_t j = 1; j <= copy->GetNbinsY(); j++) { 
+       Double_t c = copy->GetBinContent(i, j);
+       Double_t e = copy->GetBinError(i, j);
+       Double_t a = norm->GetBinContent(i);
+       copy->SetBinContent(i, j, a <= 0 ? 0 : c / a);
+       copy->SetBinError(i, j, a <= 0 ? 0 : e / a);
+      }
+    }
+
+    TH1D* res  =static_cast<TH1D*>(copy->ProjectionX("dndeta",1,
+                                                    h->GetNbinsY(),"e"));
+    TH1D* proj =static_cast<TH1D*>(h->ProjectionX("proj",1,h->GetNbinsY(),"e"));
+    res->SetTitle(*ptr);
+    res->Scale(1., "width");
+    copy->Scale(1., "width");
+    proj->Scale(1. / norm->GetMaximum(), "width");
+    norm->Scale(1. / norm->GetMaximum());
+
+    res->SetMarkerStyle(style);
+    norm->SetDirectory(0);
+    res->SetDirectory(0);
+    proj->SetDirectory(0);
+    thisList->Add(norm);
+    thisList->Add(res);
+    thisList->Add(proj);
+    dndetaRings->Add(res);
+    ptr++;
+  }
+  out->Add(dndetaRings);
+}
+
 //____________________________________________________________________
 void
 AliForwardMultiplicityBase::Print(Option_t* option) const
index 270b1369e4e80697e8a7670e612c6befece01889..558c6ca301487965dd0da3bb057f7352e3c42c0c 100644 (file)
@@ -303,7 +303,18 @@ protected:
    * 
    */
   virtual void MarkEventForStore() const;
-
+  /** 
+   * Make Ring @f$ dN/d\eta @f$ histogram and a stack 
+   * 
+   * @param input      List with summed signals 
+   * @param output     Output list 
+   * @param stackName  Stack name 
+   */
+  virtual void MakeRingdNdeta(const TList* input, 
+                             const char*  inName,
+                             TList*       output,
+                             const char*  outName,
+                             Int_t        style=20) const;
   Bool_t                 fEnableLowFlux;// Whether to use low-flux specific code
   Bool_t                 fFirstEvent;   // Whether the event is the first seen 
 private:
index ca348400aceb3f1d97e629c4a786f1cf62ae5327..9b445725f315bb905a557e423361b4b715f7a3d7 100644 (file)
@@ -34,6 +34,7 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask()
     fESDFMD(),
     fHistos(),
     fAODFMD(),
+    fRingSums(),
     fEventInspector(),
     fSharingFilter(),
     fDensityCalculator(),
@@ -53,6 +54,7 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
     fESDFMD(),
     fHistos(),
     fAODFMD(false),
+    fRingSums(),
     fEventInspector("event"),
     fSharingFilter("sharing"), 
     fDensityCalculator("density"),
@@ -76,6 +78,7 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplic
     fESDFMD(o.fESDFMD),
     fHistos(o.fHistos),
     fAODFMD(o.fAODFMD),
+    fRingSums(o.fRingSums),
     fEventInspector(o.fEventInspector),
     fSharingFilter(o.fSharingFilter),
     fDensityCalculator(o.fDensityCalculator),
@@ -115,6 +118,7 @@ AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
   fHistCollector     = o.fHistCollector;
   fHistos            = o.fHistos;
   fAODFMD            = o.fAODFMD;
+  fRingSums          = o.fRingSums;
   fList              = o.fList;
 
   return *this;
@@ -152,13 +156,31 @@ AliForwardMultiplicityTask::InitializeSubs()
 
   fHistos.Init(*pe);
   fAODFMD.Init(*pe);
+  fRingSums.Init(*pe);
 
   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
   fHData->SetStats(0);
   fHData->SetDirectory(0);
   fList->Add(fHData);
 
+  TList* rings = new TList;
+  rings->SetName("ringSums");
+  rings->SetOwner();
+  fList->Add(rings);
+
+  rings->Add(fRingSums.Get(1, 'I'));
+  rings->Add(fRingSums.Get(2, 'I'));
+  rings->Add(fRingSums.Get(2, 'O'));
+  rings->Add(fRingSums.Get(3, 'I'));
+  rings->Add(fRingSums.Get(3, 'O'));
+  fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
+  fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
+  fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
+  fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
+  fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
+
   fEventInspector.Init(*pv);
+  fSharingFilter.Init();
   fDensityCalculator.Init(*pe);
   fCorrections.Init(*pe);
   fHistCollector.Init(*pv,*pe);
@@ -215,13 +237,14 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
   fESDFMD.Clear();
   fAODFMD.Clear();
   
-  Bool_t   lowFlux  = kFALSE;
-  UInt_t   triggers = 0;
-  UShort_t ivz      = 0;
-  Double_t vz       = 0;
-  Double_t cent     = -1;
-  UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, 
-                                             ivz, vz, cent);
+  Bool_t   lowFlux   = kFALSE;
+  UInt_t   triggers  = 0;
+  UShort_t ivz       = 0;
+  Double_t vz        = 0;
+  Double_t cent      = -1;
+  UShort_t nClusters = 0;
+  UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
+                                              ivz, vz, cent, nClusters);
   
   if (found & AliFMDEventInspector::kNoEvent)    return;
   if (found & AliFMDEventInspector::kNoTriggers) return;
@@ -231,6 +254,7 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
   fAODFMD.SetSNN(fEventInspector.GetEnergy());
   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
   fAODFMD.SetCentrality(cent);
+  fAODFMD.SetNClusters(nClusters);
   MarkEventForStore();
   
   if (found & AliFMDEventInspector::kNoSPD)      return;
@@ -267,7 +291,8 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
     return;
   }
 
-  if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
+  if (!fHistCollector.Collect(fHistos, fRingSums, 
+                             ivz, fAODFMD.GetHistogram())) {
     AliWarning("Histogram collector failed");
     return;
   }
@@ -329,6 +354,8 @@ AliForwardMultiplicityTask::Terminate(Option_t*)
   list->Add(dNdeta);
   list->Add(norm);
 
+  MakeRingdNdeta(list, "ringSums", list, "ringResults");
+
   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
   fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
index 1d784b82c386457703186314292fe14ce5eedf58..1bbd5e3a265f806e11d1a63add3ca67f36f2bd7e 100644 (file)
@@ -179,6 +179,7 @@ protected:
   AliESDFMD              fESDFMD;       // Sharing corrected ESD object
   AliForwardUtil::Histos fHistos;       // Cache histograms 
   AliAODForwardMult      fAODFMD;       // Output object
+  AliForwardUtil::Histos fRingSums;     // Cache histograms 
 
   AliFMDEventInspector    fEventInspector;    // Algorithm
   AliFMDSharingFilter     fSharingFilter;     // Algorithm
index e3d29b89ca8f7af0534660d84b8e999ae3613192..c14619d8efb6163ad5092c420ede3c3cdd208860 100644 (file)
@@ -842,6 +842,7 @@ AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
   //
   if (!d) return 0;
   TList* list = new TList;
+  list->SetOwner();
   list->SetName(fName.Data());
   d->Add(list);
   return list;
index d7f103574f07d483bc26f434e1487347b08c587d..59d623d7cce5cefaf3e1aa91562866b6deeb8dd4 100644 (file)
@@ -43,7 +43,7 @@ public:
   static Color_t RingColor(UShort_t d, Char_t r)
   { 
     return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
-           + ((r == 'I' || r == 'i') ? 2 : -2));
+           + ((r == 'I' || r == 'i') ? 2 : -3));
   }
   //==================================================================
   /** 
@@ -606,6 +606,7 @@ public:
     { 
       return AliForwardUtil::RingColor(fDet, fRing);
     }
+    const char* GetName() const { return fName.Data(); } 
     UShort_t fDet;   // Detector
     Char_t   fRing;  // Ring
     TString  fName;  // Name
index 0b9484d4541cdd1312ba1bedd8200959b3a63874..580b7014ecda6fa88538104bbf46c66b17d4d38b 100644 (file)
@@ -5,6 +5,7 @@
 #include <TH1D.h>
 #include <THStack.h>
 #include <TList.h>
+#include <TFile.h>
 #include <AliAnalysisManager.h>
 #include <AliAODEvent.h>
 #include <AliAODHandler.h>
@@ -99,6 +100,57 @@ AliForwarddNdetaTask::UserExec(Option_t* option)
     bin->ProcessPrimary(forward, fTriggerMask, fVtxMin, fVtxMax, primary);
 }  
   
+//________________________________________________________________________
+void 
+AliForwarddNdetaTask::Terminate(Option_t *option) 
+{
+  // 
+  // Called at end of event processing.. 
+  //
+  // This is called once in the master 
+  // 
+  // Parameters:
+  //    option Not used 
+  AliBasedNdetaTask::Terminate(option);
+
+  THStack* truth      = new THStack("dndetaTruth",      "dN/d#eta MC Truth");
+  THStack* truthRebin = new THStack(Form("dndetaTruth_rebin%02d", fRebin), 
+                                   "dN/d#eta MC Truth");
+
+  TIter next(fListOfCentralities);
+  CentralityBin* bin = 0;
+  while ((bin = static_cast<CentralityBin*>(next()))) {
+    if (fCentAxis && bin->IsAllBin()) continue;
+
+    TList* results = bin->GetResults();
+    if (!results) continue; 
+
+    TH1* dndeta      = static_cast<TH1*>(results->FindObject("dndetaTruth"));
+    TH1* dndetaRebin = 
+      static_cast<TH1*>(results->FindObject(Form("dndetaTruth_rebin%02d",
+                                               fRebin)));
+    if (dndeta)      truth->Add(dndeta);
+    if (dndetaRebin) truthRebin->Add(dndetaRebin);
+  }
+  // If available output rebinned stack 
+  if (!truth->GetHists() || 
+      truth->GetHists()->GetEntries() <= 0) {
+    AliWarning("No MC truth histograms found");
+    delete truth;
+    truth = 0;
+  }
+  if (truth) fOutput->Add(truth);
+
+  // If available output rebinned stack 
+  if (!truthRebin->GetHists() || 
+      truthRebin->GetHists()->GetEntries() <= 0) {
+    AliWarning("No rebinned MC truth histograms found");
+    delete truthRebin;
+    truthRebin = 0;
+  }
+  if (truthRebin) fOutput->Add(truthRebin);
+  
+}
 
 //____________________________________________________________________
 TH2D*
@@ -132,8 +184,8 @@ void
 AliForwarddNdetaTask::CentralityBin::ProcessPrimary(const AliAODForwardMult* 
                                                    forward, 
                                                    Int_t triggerMask,
-                                                   Double_t /*vzMin*/
-                                                   Double_t /*vzMax*/
+                                                   Double_t vzMin
+                                                   Double_t vzMax
                                                    const TH2D* primary)
 { 
   // Check the centrality class unless this is the 'all' bin 
@@ -152,13 +204,18 @@ AliForwarddNdetaTask::CentralityBin::ProcessPrimary(const AliAODForwardMult*
   
   // translate real trigger mask to MC trigger mask
   Int_t mask = AliAODForwardMult::kB;
-  if (triggerMask == AliAODForwardMult::kNSD)
-    mask = AliAODForwardMult::kMCNSD;
+  if (triggerMask == AliAODForwardMult::kNSD) {
+    mask ^= AliAODForwardMult::kNSD;
+    mask =  AliAODForwardMult::kMCNSD;
+  }
 
   // Now use our normal check, but with the new mask, except 
-  if (!forward->CheckEvent(mask, -1000, -1000, 0, 0, 0)) return;
+  vzMin = vzMax = -10000; // ignore vertex 
+  if (!forward->CheckEvent(mask, vzMin, vzMax, 0, 0, 0)) return;
 
   fSumPrimary->Add(primary);
+  Int_t n = fSumPrimary->GetBinContent(0,0);
+  fSumPrimary->SetBinContent(0,0, ++n);
 }
 
 //________________________________________________________________________
@@ -174,7 +231,6 @@ AliForwarddNdetaTask::CentralityBin::End(TList*      sums,
                                         Bool_t      corrEmpty, 
                                         Bool_t      cutEdges,
                                         Int_t       triggerMask,
-                                        Int_t       color,
                                         Int_t       marker)
 {
   AliInfo(Form("In End of %s with corrEmpty=%d, cutEdges=%d, rootProj=%d", 
@@ -183,28 +239,73 @@ AliForwarddNdetaTask::CentralityBin::End(TList*      sums,
                                        shapeCorr, trigEff, 
                                        symmetrice, rebin, 
                                        rootProj, corrEmpty, cutEdges,
-                                       triggerMask, color, marker);
+                                       triggerMask, marker);
 
   fSumPrimary     = static_cast<TH2D*>(fSums->FindObject("truth"));
 
   if (!fSumPrimary) { 
     AliWarning("No MC truth histogram found");
-    fSums->ls();
-    return;
   }
-  Int_t n = (triggerMask == AliAODForwardMult::kNSD ? 
-            Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinMCNSD)) : 
-            Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinAll)));
-  AliInfo(Form("Normalising MC truth to %d", n));
-
-  TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1,
+  else {
+#if 0
+    Int_t n = fSumPrimary->GetBinContent(0,0);
+#else
+    Int_t n = (triggerMask == AliAODForwardMult::kNSD ? 
+              Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinMCNSD)) : 
+              Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinAll)));
+#endif
+    AliInfo(Form("Normalising MC truth to %d", n));
+    
+    TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1,
                                               fSumPrimary->GetNbinsY(),"e");
-  dndetaTruth->Scale(1./n, "width");
+    dndetaTruth->SetDirectory(0);
+    dndetaTruth->Scale(1./n, "width");
 
-  SetHistogramAttributes(dndetaTruth, color-2, marker+4, "Monte-Carlo truth");
     
-  fOutput->Add(dndetaTruth);
-  fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges));
+    SetHistogramAttributes(dndetaTruth, GetColor(), 30, "Monte-Carlo truth");
+    
+    fOutput->Add(dndetaTruth);
+    fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges));
+  }
+
+  if (!IsAllBin()) return;
+  TFile* file = TFile::Open("forward.root", "READ");
+  if (!file) return;
+  
+  TList* forward = static_cast<TList*>(file->Get("Forward"));
+  if (!forward) { 
+    AliError("List Forward not found in forward.root");
+    return;
+  }
+  TList* rings = static_cast<TList*>(forward->FindObject("ringResults"));
+  if (!rings) { 
+    AliError("List ringResults not found in forward.root");
+    return;
+  }
+  THStack* res = static_cast<THStack*>(rings->FindObject("all"));
+  if (!res) { 
+    AliError(Form("Stack all not found in %s", rings->GetName()));
+    return;
+  }
+  if (!fTriggers) { 
+    AliError("Triggers histogram not set");
+    return;
+  }
+  Double_t ntotal   = 0;
+  Double_t epsilonT = trigEff;
+  // TEMPORARY FIX
+  if (triggerMask == AliAODForwardMult::kNSD) {
+    // This is a local change 
+    epsilonT = 0.92; 
+    AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
+  }
+  AliInfo("Adding per-ring histograms to output");
+  Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal);
+  TIter next(res->GetHists());
+  TH1*  hist = 0;
+  while ((hist = static_cast<TH1*>(next()))) hist->Scale(scaler);
+  res->SetName("dndetaRings");
+  fOutput->Add(res);
 }
 
 //________________________________________________________________________
index 0859946f21dc7b1dc815ca20ade341153f1b244a..309e4b0abba6acd42caf6ff02c137a301fed8587 100644 (file)
@@ -51,6 +51,14 @@ public:
    * @param option Not used 
    */
   virtual void UserExec(Option_t* option);
+  /** 
+   * Called at end of event processing.
+   *
+   * This is called once in the master 
+   * 
+   * @param option Not used 
+   */
+  virtual void Terminate(Option_t* option);
 protected:
   /** 
    * Copy constructor 
@@ -162,7 +170,6 @@ protected:
                     Bool_t      corrEmpty, 
                     Bool_t      cutEdges, 
                     Int_t       triggerMask,
-                    Int_t       color, 
                     Int_t       marker);
   protected: 
     TH2D*           fSumPrimary;    //  Sum of primary histograms
index c6afac18bffa7f155cc3488ac97c70cf68c21558..33086534ac670b46bb2be934354eb3bc298d6ff0 100644 (file)
@@ -3,11 +3,13 @@
  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
  * @date   Wed Mar 23 14:07:10 2011
  * 
- * @brief  Script to visualise the dN/deta 
+ * @brief  Script to visualise the dN/deta for pp and PbPb
  *
  * This script is independent of any AliROOT code - and should stay
  * that way.
  * 
+ * The script is <emph>very</emph> long - sigh - the joy of drawing
+ * things nicely in ROOT
  * 
  * @ingroup pwg2_forward_dndeta
  */
@@ -377,10 +379,12 @@ struct dNdetaDrawer
                                             sys,snn,trg,
                                             centLow,centHigh,onlya));
     if (!ret) { 
+#if 0
       Warning("FetchResults", "No other data found for sys=%d, sNN=%d, "
              "trigger=%d %d%%-%d%% central %s", 
              sys, snn, trg, centLow, centHigh, 
              onlya ? " (ALICE results)" : "all");
+#endif
       return 0;
     }
     thisOther = reinterpret_cast<TMultiGraph*>(ret);
@@ -410,22 +414,16 @@ struct dNdetaDrawer
        Error("FetchResults", "Couldn't find list 'all' in %s", 
              list->GetName());
       else 
-       FetchResults(all, name, FetchOthers(0,0), -1, 0, max, rmax, amax);
+       FetchResults(all, name, FetchOthers(0,0), -1000, 0, max, rmax, amax);
       return;
     }
     
-    Int_t   nCol = gStyle->GetNumberOfColors();
     for (UShort_t i = 0; i < n; i++) { 
-      UShort_t centLow  = fCentAxis->GetXbins()->At(i);
-      UShort_t centHigh = fCentAxis->GetXbins()->At(i+1);
+      UShort_t centLow  = fCentAxis->GetBinLowEdge(i+1);
+      UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
       TString  lname    = Form("cent%03d_%03d", centLow, centHigh);
       TList*   thisCent = static_cast<TList*>(list->FindObject(lname));
-
-      Float_t fc   = (centLow+double(centHigh-centLow)/2) / 100;
-      Int_t   icol = TMath::Min(nCol-1,int(fc * nCol + .5));
-      Int_t   col  = gStyle->GetColorPalette(icol);
-      Info("FetchResults","Centrality %d-%d color index %d (=%d*%f) -> %d", 
-          centLow, centHigh, icol, nCol, fc, col);
+      Int_t    col      = GetCentralityColor(i+1);
 
       TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
       if (!thisCent)
@@ -437,6 +435,18 @@ struct dNdetaDrawer
     }
   } 
   //__________________________________________________________________
+  Int_t GetCentralityColor(Int_t bin) const
+  {
+    UShort_t centLow  = fCentAxis->GetBinLowEdge(bin);
+    UShort_t centHigh = fCentAxis->GetBinUpEdge(bin);
+    Float_t  fc       = (centLow+double(centHigh-centLow)/2) / 100;
+    Int_t    nCol     = gStyle->GetNumberOfColors();
+    Int_t    icol     = TMath::Min(nCol-1,int(fc * nCol + .5));
+    Int_t    col      = gStyle->GetColorPalette(icol);
+    //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
+    return col;
+  }
+  //__________________________________________________________________
   void SetAttributes(TH1* h, Int_t color)
   {
     if (!h) return;
@@ -457,8 +467,9 @@ struct dNdetaDrawer
   //__________________________________________________________________
   void ModifyTitle(TNamed* h, const char* centTxt)
   {
-    if (!centTxt || !h) return;
-    h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
+    return;
+    // if (!centTxt || !h) return;
+    // h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
   }
 
   //__________________________________________________________________
@@ -491,11 +502,17 @@ struct dNdetaDrawer
     TH1* tracks      = FetchResult(list, "tracks");
     if (tracks) tracks->SetTitle("ALICE Tracks");
     SetAttributes(dndeta,     color);
-    SetAttributes(dndetaMC,   color+2);
+    SetAttributes(dndetaMC,   fCentAxis ? color : color+2);
     SetAttributes(dndetaTruth,color);
     SetAttributes(dndetaSym,  color);
-    SetAttributes(dndetaMCSym,color+2);
-    SetAttributes(tracks,     color+3);
+    SetAttributes(dndetaMCSym,fCentAxis ? color : color+2);
+    SetAttributes(tracks,     fCentAxis ? color : color+3);
+    if (dndetaMC && fCentAxis) 
+      dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
+    if (dndetaMCSym && fCentAxis) 
+      dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2);
+    if (dndetaTruth && fCentAxis) 
+      dndetaTruth->SetMarkerStyle(34);
     ModifyTitle(dndeta,     centTxt);
     ModifyTitle(dndetaMC,   centTxt);
     ModifyTitle(dndetaTruth,centTxt);
@@ -508,7 +525,14 @@ struct dNdetaDrawer
     max = TMath::Max(max, AddHistogram(fResults, dndetaMC,    dndetaMCSym));
     max = TMath::Max(max, AddHistogram(fResults, dndeta,      dndetaSym));
     max = TMath::Max(max, AddHistogram(fResults, tracks));
-    
+
+    THStack* rings = static_cast<THStack*>(list->FindObject("dndetaRings"));
+    if (rings) { 
+      TIter next(rings->GetHists());
+      TH1*  hist = 0;
+      while ((hist = static_cast<TH1*>(next()))) 
+       max = TMath::Max(max, AddHistogram(fResults, hist));
+    }
     // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f", 
     //      dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
 
@@ -540,11 +564,13 @@ struct dNdetaDrawer
        fRatios->Add(Ratio(dndeta, tracks, rmax));
     }
 
+    if (dndetaMC) { 
+      fRatios->Add(Ratio(dndeta,    dndetaMC,    rmax));
+      fRatios->Add(Ratio(dndetaSym, dndetaMCSym, rmax));
+    }
     if (dndetaTruth) { 
       fRatios->Add(Ratio(dndeta,      dndetaTruth, rmax));
       fRatios->Add(Ratio(dndetaSym,   dndetaTruth, rmax));
-      fRatios->Add(Ratio(dndetaMC,    dndetaTruth, rmax));
-      fRatios->Add(Ratio(dndetaMCSym, dndetaTruth, rmax));
     }
   }
   //__________________________________________________________________
@@ -580,12 +606,6 @@ struct dNdetaDrawer
     c->SetBorderSize(0);
     c->SetBorderMode(0);
 
-#if 1
-    Info("Plot", "y1=%f, y2=%f, y3=%f extra: %s %s", 
-        y1, y2, y2, (fShowRatios ? "ratios" : ""), 
-        (fShowLeftRight ? "right/left" : ""));
-    Info("Plot", "Maximum is %f", max);
-#endif
     PlotResults(max, y1);
     c->cd();
 
@@ -628,27 +648,54 @@ struct dNdetaDrawer
                   Double_t x1, Double_t y1, Double_t x2, Double_t y2)
   {
     TLegend* l = new TLegend(x1,y1,x2,y2);
-    l->SetNColumns(2);
+    l->SetNColumns(fCentAxis ? 1 : 2);
     l->SetFillColor(0);
     l->SetFillStyle(0);
     l->SetBorderSize(0);
     l->SetTextFont(132);
 
     TIter    next(stack->GetHists());
-    TObject* hist = 0;
-    while ((hist = next())) { 
+    TH1*     hist = 0;
+    TObjArray unique;
+    unique.SetOwner();
+    while ((hist = static_cast<TH1*>(next()))) { 
       TString n(hist->GetTitle());
       if (n.Contains("mirrored")) continue;
-      l->AddEntry(hist, hist->GetTitle(), "pl");
+      if (unique.FindObject(n.Data())) continue;
+      TObjString* s = new TObjString(n);
+      s->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
+                    ((hist->GetMarkerColor() & 0xFFFF) <<  0));
+      unique.Add(s);
+      // l->AddEntry(hist, hist->GetTitle(), "pl");
     }
     if (mg) {
       TIter nexto(mg->GetListOfGraphs());
-      while ((hist = nexto())) { 
-       TString n(hist->GetTitle());
+      TGraph* g = 0;
+      while ((g = static_cast<TGraph*>(nexto()))) { 
+       TString n(g->GetTitle());
        if (n.Contains("mirrored")) continue;
-       l->AddEntry(hist, hist->GetTitle(), "pl");
+       if (unique.FindObject(n.Data())) continue;
+       TObjString* s = new TObjString(n);
+       s->SetUniqueID(((g->GetMarkerStyle() & 0xFFFF) << 16) |
+                      ((g->GetMarkerColor() & 0xFFFF) <<  0));
+       unique.Add(s);
+       // l->AddEntry(hist, hist->GetTitle(), "pl");
       }
     }
+    unique.ls();
+    TIter nextu(&unique);
+    TObject* s = 0;
+    Int_t i = 0;
+    while ((s = nextu())) { 
+      TLegendEntry* dd = l->AddEntry(Form("data%2d", i++), 
+                                    s->GetName(), "lp");
+      Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF;
+      Int_t color = (s->GetUniqueID() >>  0) & 0xFFFF;
+      dd->SetLineColor(kBlack);
+      if (fCentAxis) dd->SetMarkerColor(kBlack);
+      else           dd->SetMarkerColor(color);
+      dd->SetMarkerStyle(style);
+    }
     TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
     d1->SetLineColor(kBlack);
     d1->SetMarkerColor(kBlack);
@@ -661,6 +708,38 @@ struct dNdetaDrawer
     l->Draw();
   }
   //__________________________________________________________________
+  /** 
+   * Build centrality legend 
+   * 
+   * @param axis    Stack to include 
+   * @param x1      Lower X coordinate in the range [0,1]
+   * @param y1      Lower Y coordinate in the range [0,1]
+   * @param x2      Upper X coordinate in the range [0,1]
+   * @param y2             Upper Y coordinate in the range [0,1]
+   */
+  void BuildCentLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
+  {
+    if (!fCentAxis) return;
+
+    TLegend* l = new TLegend(x1,y1,x2,y2);
+    l->SetNColumns(1);
+    l->SetFillColor(0);
+    l->SetFillStyle(0);
+    l->SetBorderSize(0);
+    l->SetTextFont(132);
+
+    Int_t n = fCentAxis->GetNbins();
+    for (Int_t i = 1; i <= n; i++) { 
+      Double_t low = fCentAxis->GetBinLowEdge(i);
+      Double_t upp = fCentAxis->GetBinUpEdge(i);
+      TLegendEntry* e = l->AddEntry(Form("dummy%02d", i),
+                                   Form("%3d%% - %3d%%", 
+                                        int(low), int(upp)), "pl");
+      e->SetMarkerColor(GetCentralityColor(i));
+    }
+    l->Draw();
+  }
+  //__________________________________________________________________
   /** 
    * Plot the results
    *    
@@ -703,15 +782,14 @@ struct dNdetaDrawer
     }
 
     // Make a legend in the main result pad
-    BuildLegend(fResults, fOthers, .15,p1->GetBottomMargin()+.01,.90,.35);
-#if 0
-    TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
-    l->SetNColumns(2);
-    l->SetFillColor(0);
-    l->SetFillStyle(0);
-    l->SetBorderSize(0);
-    l->SetTextFont(132);
-#endif
+    BuildCentLegend(.12, 1-p1->GetTopMargin()-.01-.5,  
+                   .35, 1-p1->GetTopMargin()-.01-.1);
+    Double_t x1 = (fCentAxis ? .7 : .15); 
+    Double_t x2 = (fCentAxis ? 1-p1->GetRightMargin()-.01: .90);
+    Double_t y1 = (fCentAxis ? .5: p1->GetBottomMargin()+.01); 
+    Double_t y2 = (fCentAxis ? 1-p1->GetTopMargin()-.01-.15 : .35);
+                  
+    BuildLegend(fResults, fOthers, x1, y1, x2, y2);
 
     // Put a title on top
     fTitle.ReplaceAll("@", " ");
@@ -725,11 +803,14 @@ struct dNdetaDrawer
     TString     eS;
     UShort_t    snn = fSNNString->GetUniqueID();
     const char* sys = fSysString->GetTitle();
+    if (snn == 2750) snn = 2760;
     if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
     else            eS = Form("%3dGeV", snn);
-    TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s", 
+    TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s%s}=%s, %s", 
                                           sys, 
+                                          (fCentAxis ? "_{NN}" : ""),
                                           eS.Data(), 
+                                          fCentAxis ? "by centrality" : 
                                           fTrigString->GetTitle()));
     tt->SetNDC();
     tt->SetTextFont(132);
@@ -886,12 +967,17 @@ struct dNdetaDrawer
 
     
     // Make a legend
-    TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
-    l2->SetNColumns(2);
-    l2->SetFillColor(0);
-    l2->SetFillStyle(0);
-    l2->SetBorderSize(0);
-    l2->SetTextFont(132);
+    Double_t xx1 = (fCentAxis ? .7                           : .15); 
+    Double_t xx2 = (fCentAxis ? 1-p3->GetRightMargin()-.01   : .90);
+    Double_t yy1 = p3->GetBottomMargin()+.01;
+    Double_t yy2 = (fCentAxis ? 1-p3->GetTopMargin()-.01-.15 : .5);
+    BuildLegend(fLeftRight, 0, xx1, yy1, xx2, yy2);
+    // TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
+    // l2->SetNColumns(2);
+    // l2->SetFillColor(0);
+    // l2->SetFillStyle(0);
+    // l2->SetBorderSize(0);
+    // l2->SetTextFont(132);
 
     // Make a nice band from 0.9 to 1.1
     TGraphErrors* band = new TGraphErrors(2);
@@ -936,10 +1022,12 @@ struct dNdetaDrawer
     if (!list) return 0;
     
     TH1* ret = static_cast<TH1*>(list->FindObject(name));
+#if 0
     if (!ret) {
       // all->ls();
       Warning("GetResult", "Histogram %s not found", name);
     }
+#endif
 
     return ret;
   }
@@ -961,10 +1049,7 @@ struct dNdetaDrawer
     // Rebin if needed 
     Rebin(hist);
 
-    // Info("AddHistogram", "Adding %s to %s", 
-    //      hist->GetName(), stack->GetName());
     stack->Add(hist, option);
-    // stack->ls();
     return hist->GetMaximum();
   }
   //__________________________________________________________________
@@ -992,8 +1077,6 @@ struct dNdetaDrawer
     sym = Symmetrice(hist);
     stack->Add(sym, option);
 
-    // Info("AddHistogram", "Adding %s and %s to %s", 
-    //      hist->GetName(), sym->GetName(), stack->GetName());
     return hist->GetMaximum();
   }
 
@@ -1239,32 +1322,34 @@ struct dNdetaDrawer
   {
     if (!o1 || !o2) return 0;
 
-    TH1* r = 0;
+    TH1*        r  = 0;
+    const TAttMarker* m1 = 0;
+    const TAttMarker* m2 = 0;
     const TH1* h1 = dynamic_cast<const TH1*>(o1); 
     if (h1) { 
-      // Info("Ratio", "First is a TH1");
+      m1 = h1;
       const TH1* h2 = dynamic_cast<const TH1*>(o2); 
       if (h2) { 
-       // Info("Ratio", "Second is a TH1");
-       r = RatioHH(h1,h2,max);
+       m2 = h2;
+       r = RatioHH(h1,h2);
       }
       else {
        const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
        if (g2) { 
-         // Info("Ratio", "Second os a TGraph");
-         r = RatioHG(h1,g2,max);      
+         m2 = g2;
+         r = RatioHG(h1,g2);      
        }
       }
     }
     else {
       const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
       if (g1) { 
-       // Info("Ratio", "First is a TGraphAsymmErrors");
+       m1 = g1;
        const TGraphAsymmErrors* g2 = 
          dynamic_cast<const TGraphAsymmErrors*>(o2);
        if (g2) {
-         // Info("Ratio", "Second is a TGraphAsymmErrors");
-         r = RatioGG(g1, g2, max);
+         m2 = g2;
+         r = RatioGG(g1, g2);
        }
       }
     }
@@ -1278,9 +1363,16 @@ struct dNdetaDrawer
       delete r; 
       r = 0; 
     }
-    // for (Int_t bin = 1; bin <= r->GetNbinsX(); bin++) 
-    //   if (r->GetBinContent(bin) != 0) return r;
-      
+    if (r) {
+      r->SetMarkerStyle(m1->GetMarkerStyle());
+      r->SetMarkerColor(m2->GetMarkerColor());
+      r->SetMarkerSize(0.9*m1->GetMarkerSize());
+      r->SetName(Form("%s_over_%s", o1->GetName(), o2->GetName()));
+      r->SetTitle(Form("%s / %s", o1->GetTitle(), o2->GetTitle()));
+      r->SetDirectory(0);
+      max = TMath::Max(RatioMax(r), max);
+    }
+
     return r;
   }
   //__________________________________________________________________
@@ -1290,25 +1382,15 @@ struct dNdetaDrawer
    * 
    * @param h  Numerator 
    * @param g  Divisor 
-   * @param max Maximum diviation from 1 
    * 
    * @return h/g 
    */
-  TH1* RatioHG(const TH1* h, const TGraph* g, Double_t& max) const 
+  TH1* RatioHG(const TH1* h, const TGraph* g) const 
   {
     if (!h || !g) return 0;
 
     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
-    ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
-    ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
     ret->Reset();
-    ret->SetMarkerStyle(h->GetMarkerStyle());
-    ret->SetMarkerColor(g->GetMarkerColor());
-    ret->SetMarkerSize(0.9*h->GetMarkerSize());
-    // ret->SetMarkerStyle(g->GetMarkerStyle());
-    // ret->SetMarkerColor(h->GetMarkerColor());
-    // ret->SetMarkerSize(0.9*g->GetMarkerSize());
-    ret->SetDirectory(0);
     Double_t xlow  = g->GetX()[0];
     Double_t xhigh = g->GetX()[g->GetN()-1];
     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
@@ -1326,8 +1408,6 @@ struct dNdetaDrawer
       ret->SetBinContent(i, c / f);
       ret->SetBinError(i, h->GetBinError(i) / f);
     }
-    if (ret->GetEntries() > 0) 
-      max = TMath::Max(RatioMax(ret), max);
     return ret;
   }
   //__________________________________________________________________
@@ -1336,25 +1416,14 @@ struct dNdetaDrawer
    * 
    * @param h1 First histogram (numerator) 
    * @param h2 Second histogram (denominator)
-   * @param max Maximum diviation from 1 
    * 
    * @return h1 / h2
    */
-  TH1* RatioHH(const TH1* h1, const TH1* h2, Double_t& max) const
+  TH1* RatioHH(const TH1* h1, const TH1* h2) const
   {
     if (!h1 || !h2) return 0;
-    TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s", 
-                                              h1->GetName(), 
-                                              h2->GetName())));
-    t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle()));
+    TH1* t1 = static_cast<TH1*>(h1->Clone("tmp"));
     t1->Divide(h2);
-    // t1->SetMarkerColor(h1->GetMarkerColor());
-    // t1->SetMarkerStyle(h2->GetMarkerStyle());
-    t1->SetMarkerColor(h2->GetMarkerColor());
-    t1->SetMarkerStyle(h1->GetMarkerStyle());
-    t1->SetMarkerSize(0.9*h1->GetMarkerSize());
-    t1->SetDirectory(0);
-    max = TMath::Max(RatioMax(t1), max);
     return t1;
   }
   //__________________________________________________________________
@@ -1363,12 +1432,11 @@ struct dNdetaDrawer
    * 
    * @param g1 Numerator 
    * @param g2 Denominator
-   * @param max Maximum diviation from 1 
    * 
    * @return g1 / g2 in a histogram 
    */
   TH1* RatioGG(const TGraphAsymmErrors* g1, 
-              const TGraphAsymmErrors* g2, Double_t& max) const
+              const TGraphAsymmErrors* g2) const
   {
     Int_t    nBins = g1->GetN();
     TArrayF  bins(nBins+1);
@@ -1383,22 +1451,13 @@ struct dNdetaDrawer
       if (i == 0) dx  = dxi;
       else if (dxi != dx) dx = 0;
     }
-    TString name(Form("%s_%s", g1->GetName(), g2->GetName()));
-    TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle()));
     TH1* h = 0;
     if (dx != 0) {
-      h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
+      h = new TH1F("tmp", "tmp", nBins, bins[0], bins[nBins]);
     }
     else {
-      h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray);
+      h = new TH1F("tmp", "tmp", nBins, bins.fArray);
     }
-    h->SetMarkerStyle(g1->GetMarkerStyle());
-    h->SetMarkerColor(g2->GetMarkerColor());
-    h->SetMarkerSize(0.9*g1->GetMarkerSize());
-    // h->SetMarkerStyle(g2->GetMarkerStyle());
-    // h->SetMarkerColor(g1->GetMarkerColor());
-    // h->SetMarkerSize(0.9*g2->GetMarkerSize());
-    h->SetDirectory(0);
 
     Double_t low  = g2->GetX()[0];
     Double_t high = g2->GetX()[g2->GetN()-1];
@@ -1413,7 +1472,6 @@ struct dNdetaDrawer
       h->SetBinContent(i+1, c1 / c2);
       h->SetBinError(i+1, e1 / c2);
     }
-    max = TMath::Max(RatioMax(h), max);
     return h;
   }  
   /* @} */
@@ -1543,7 +1601,6 @@ void UpdateRange(dNdetaDrawer::RangeParam* p)
   Int_t    last  = p->fMasterAxis->GetLast();
   Double_t x1    = p->fMasterAxis->GetBinCenter(first);
   Double_t x2    = p->fMasterAxis->GetBinCenter(last);
-  //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2);
 
   if (p->fSlave1Axis) { 
     Int_t i1 = p->fSlave1Axis->FindBin(x1);
index 89528561dc1b8076944aeb17c795f3ce6d6a1885..50a8e15fe92e873269345b0e6d987a0e742a7326 100644 (file)
@@ -28,6 +28,22 @@ ForwardAODConfig(AliForwardMultiplicityBase* task)
   // Whether to enable low flux specific code 
   task->SetEnableLowFlux(kFALSE);
 
+  // Would like to use dynamic cast but CINT interprets that as a 
+  // static cast - sigh!
+  Bool_t mc = false;
+  if (task->IsA() == AliForwardMCMultiplicityTask::Class()) 
+    mc = true;
+
+#if 0 
+  if (mc) {
+    AliForwardMCMultiplicityTask* mcTask = 
+      static_cast<AliForwardMCMultiplicityTask*>(task);
+    mcTask->SetOnlyPrimary(true);
+  }
+#endif
+  Double_t nXi = mc ? 1 : .5;
+  Bool_t   includeSigma = true;
+
   // --- Event inspector ---------------------------------------------
   // Set the number of SPD tracklets for which we consider the event a
   // low flux event
@@ -37,17 +53,29 @@ ForwardAODConfig(AliForwardMultiplicityBase* task)
 
   // --- Sharing filter ----------------------------------------------
   // Set the low cut used for sharing - overrides settings in eloss fits
-  task->GetSharingFilter().SetLowCut(0.3);
+  task->GetSharingFilter().SetLowCut(0.15);
   // Set the number of xi's (width of landau peak) to stop at 
-  task->GetSharingFilter().SetNXi(1);
+  task->GetSharingFilter().SetNXi(nXi);
+  // Set whether or not to include sigma in cut
+  task->GetSharingFilter().SetIncludeSigma(includeSigma);
+  // Enable use of angle corrected signals in the algorithm 
+  task->GetSharingFilter().SetUseAngleCorrectedSignals(true);
 
   // --- Density calculator 
   // Set the maximum number of particle to try to reconstruct 
-  task->GetDensityCalculator().SetMaxParticles(3);
+  task->GetDensityCalculator().SetMaxParticles(10);
   // Wet whether to use poisson statistics to estimate N_ch
   task->GetDensityCalculator().SetUsePoisson(false);
   // Set the lower multiplicity cut.  Overrides setting in energy loss fits.
-  task->GetDensityCalculator().SetMultCut(0.3); //was 0.3
+  task->GetDensityCalculator().SetMultCut(-1); //was 0.3
+  // Set the lower per-ring multiplicity cuts 
+  task->GetDensityCalculator().SetMultCuts(-1,-1,-1,-1,-1);
+  // USe this many times xi+sigma below MPV 
+  task->GetDensityCalculator().SetNXi(nXi);
+  // Set whether or not to include sigma in cut
+  task->GetDensityCalculator().SetIncludeSigma(includeSigma);
+  // Set whether or not to use the phi acceptance 
+  task->GetDensityCalculator().SetUsePhiAcceptance(true);
 
   // --- Corrector ---------------------------------------------------
   // Whether to use the secondary map correction
@@ -82,7 +110,7 @@ ForwardAODConfig(AliForwardMultiplicityBase* task)
   // Set the overall debug level (1: some output, 3: a lot of output)
   task->SetDebug(0);
   // Set the debug level of a single algorithm 
-  // task->GetEventInspector().SetDebug(4);
+  task->GetSharingFilter().SetDebug(0);
 
   // --- Set limits on fits the energy -------------------------------
   // Maximum relative error on parameters 
diff --git a/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C b/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C
new file mode 100644 (file)
index 0000000..4e5481c
--- /dev/null
@@ -0,0 +1,851 @@
+#ifdef BUILD
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliMCEvent.h"
+#include "AliESDEvent.h"
+#include "AliStack.h"
+#include "AliMultiplicity.h"
+#include "AliFMDMCEventInspector.h"
+#include "AliAODForwardMult.h"
+#include "AliLog.h"
+#include <TH1I.h>
+#include <TH2D.h>
+#include <TAxis.h>
+#include <TList.h>
+#include <TObjArray.h>
+#include <TParameter.h>
+#include <TStopwatch.h>
+#include <TROOT.h>
+#include <THStack.h>
+#include <TStyle.h>
+
+//====================================================================
+/**
+ * Task to evaluate trigger bias in pp 
+ * 
+ */
+class EvaluateTrigger : public AliAnalysisTaskSE
+{
+public:
+  enum { 
+    kNone = 0x0, 
+    kESD  = 0x1, 
+    kMC   = 0x2
+  };
+  /** 
+   * Constructor
+   */
+  EvaluateTrigger() 
+    : AliAnalysisTaskSE(),
+      fInel(),
+      fInelGt0(),
+      fNSD(),
+      fInspector(), 
+      fFirstEvent(true),
+      fData(0), 
+      fTriggers(0), 
+      fTrackletRequirement(kESD),
+      fVertexRequirement(kESD), 
+      fVertexAxis(0, 0, 0), 
+      fVertexESD(0),
+      fVertexMC(0), 
+      fM(0)
+  {}
+  /** 
+   * Constructor 
+   */
+  EvaluateTrigger(const char* /*name*/) 
+    : AliAnalysisTaskSE("evaluateTrigger"),
+      fInel(AliAODForwardMult::kInel),
+      fInelGt0(AliAODForwardMult::kInelGt0),
+      fNSD(AliAODForwardMult::kNSD),
+      fInspector("eventInspector"), 
+      fFirstEvent(true), 
+      fData(0), 
+      fTriggers(0),
+      fTrackletRequirement(kESD),
+      fVertexRequirement(kESD), 
+      fVertexAxis(10, -10, 10), 
+      fVertexESD(0),
+      fVertexMC(0), 
+      fM(0)
+  {
+    DefineOutput(1, TList::Class());
+    DefineOutput(2, TList::Class());
+  }
+  void SetVertexRequirement(UShort_t m) { fVertexRequirement = m; }
+  void SetTrackletRequirement(UShort_t m) { fTrackletRequirement = m; }
+  void SetVertexAxis(Int_t nBins, Double_t low, Double_t high) 
+  {
+    fVertexAxis.Set(nBins, low, high);
+  }
+  /** 
+   * Intialize 
+   */
+  void Init() {}
+  /** 
+   * Create output objects 
+   */
+  void UserCreateOutputObjects()
+  {
+    fList = new TList;
+    fList->SetOwner();
+    fList->SetName(GetName());
+
+    Double_t mb[] = { 0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 1000 };
+    Int_t    nM   = 10;
+    TAxis mAxis(nM, mb);
+    TAxis eAxis(200, -4, 6);
+    TAxis pAxis(40, 0, 2*TMath::Pi());
+
+    fData = new TH2D("data", "Cache", 
+                    eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(), 
+                    pAxis.GetNbins(), pAxis.GetXmin(), pAxis.GetXmax());
+    fData->SetDirectory(0);
+    fData->SetXTitle("#eta");
+    fData->SetYTitle("#varphi [radians]");
+    fData->SetZTitle("N_{ch}(#eta,#varphi)");
+    fData->Sumw2();
+    
+    fM = new TH1D("m", "Distribution of N_{ch}|_{|#eta|<1}",nM+1, -0.5, nM+.5); 
+    fM->SetXTitle("N_{ch}|_{|#eta|<1}");
+    fM->SetYTitle("Events");
+    fM->SetFillColor(kRed+1);
+    fM->SetFillStyle(3001);
+    fM->SetDirectory(0);
+    fList->Add(fM);
+
+    for (Int_t i = 0; i <= nM; i++) { 
+      TString lbl;
+      if (i == 0)       lbl = "all";
+      else if (i == nM) lbl = Form("%d+",int(mAxis.GetBinLowEdge(i)+.5));
+      else              lbl = Form("<%d",int(mAxis.GetBinUpEdge(i)+.5));
+      fM->GetXaxis()->SetBinLabel(i+1, lbl);
+    }
+
+    fTriggers = new TH1I("triggers", "Triggers", 6, -.5, 5.5);
+    fTriggers->SetDirectory(0);
+    fTriggers->GetXaxis()->SetBinLabel(1, "INEL (MC)");
+    fTriggers->GetXaxis()->SetBinLabel(2, "INEL (ESD)");
+    fTriggers->GetXaxis()->SetBinLabel(3, "INEL>0 (MC)");
+    fTriggers->GetXaxis()->SetBinLabel(4, "INEL>0 (ESD)");
+    fTriggers->GetXaxis()->SetBinLabel(5, "NSD (MC)");
+    fTriggers->GetXaxis()->SetBinLabel(6, "NSD (ESD)");
+    fTriggers->SetFillColor(kYellow+1);
+    fTriggers->SetFillStyle(3001);
+    fList->Add(fTriggers);
+
+    fVertexESD = new TH1D("vertexESD", "ESD vertex distribution", 
+                         fVertexAxis.GetNbins(), 
+                         fVertexAxis.GetXmin(), 
+                         fVertexAxis.GetXmax());
+    fVertexESD->SetDirectory(0);
+    fVertexESD->SetFillColor(kRed+1);
+    fVertexESD->SetFillStyle(3001);
+    fVertexESD->SetXTitle("v_{z} [cm]");
+    fVertexESD->SetYTitle("P(v_{z})");
+    fList->Add(fVertexESD);
+
+    fVertexMC = new TH1D("vertexMC", "MC vertex distribution", 
+                         fVertexAxis.GetNbins(), 
+                         fVertexAxis.GetXmin(), 
+                         fVertexAxis.GetXmax());
+    fVertexMC->SetDirectory(0);
+    fVertexMC->SetFillColor(kBlue+1);
+    fVertexMC->SetFillStyle(3001);
+    fVertexMC->SetXTitle("v_{z} [cm]");
+    fVertexMC->SetYTitle("P(v_{z})");
+    fList->Add(fVertexMC);
+
+    fInel.CreateObjects(fList, fM, fData);
+    fInelGt0.CreateObjects(fList, fM, fData);
+    fNSD.CreateObjects(fList, fM, fData);
+
+
+    fInspector.DefineOutput(fList);
+    fInspector.Init(fVertexAxis);
+
+    PostData(1, fList);
+  }
+  /** 
+   * Event processing 
+   */
+  void UserExec(Option_t*) 
+  {
+    // Get the input data - MC event
+    AliMCEvent*  mcEvent = MCEvent();
+    if (!mcEvent) { 
+      AliWarning("No MC event found");
+      return;
+    }
+    
+    // Get the input data - ESD event
+    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
+    if (!esd) { 
+      AliWarning("No ESD event found for input event");
+      return;
+    }
+
+    if (fFirstEvent && esd->GetESDRun()) {
+      fInspector.ReadRunDetails(esd);
+
+      AliInfo(Form("Initializing with parameters from the ESD:\n"
+                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
+                  "         AliESDEvent::GetBeamType()     ->%s\n"
+                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
+                  "         AliESDEvent::GetMagneticField()->%f\n"
+                  "         AliESDEvent::GetRunNumber()    ->%d\n",
+                  esd->GetBeamEnergy(),
+                  esd->GetBeamType(),
+                  esd->GetCurrentL3(),
+                  esd->GetMagneticField(),
+                  esd->GetRunNumber()));
+      
+      fFirstEvent = false;
+    }
+
+    // Get the particle stack 
+    AliStack* stack = mcEvent->Stack();
+
+    // Some variables 
+    UInt_t   triggers; // Trigger bits
+    Bool_t   lowFlux;  // Low flux flag
+    UShort_t iVz;      // Vertex bin from ESD
+    Double_t vZ;       // Z coordinate from ESD
+    Double_t cent;     // Centrality 
+    UShort_t iVzMc;    // Vertex bin from MC
+    Double_t vZMc;     // Z coordinate of IP vertex from MC
+    Double_t b;        // Impact parameter
+    Int_t    nPart;    // Number of participants 
+    Int_t    nBin;     // Number of binary collisions 
+    Double_t phiR;     // Reaction plane from MC
+    
+    // Process the data 
+    Int_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, cent);
+    Int_t retMC  = fInspector.ProcessMC(mcEvent, triggers, iVzMc, 
+                                       vZMc, b, nPart, nBin, phiR);
+
+    Bool_t hasESDVtx = retESD == AliFMDEventInspector::kOk;
+    Bool_t hasMCVtx  = retMC  == AliFMDEventInspector::kOk;
+    if (hasESDVtx) fVertexESD->Fill(vZ);
+    if (hasMCVtx)  fVertexMC->Fill(vZMc);
+
+    Bool_t isMcInel = true; // (triggers & AliAODForwardMult::kB);
+    Bool_t isMcNSD  = (triggers & AliAODForwardMult::kMCNSD);
+
+    Int_t mESD = 0;
+    const AliMultiplicity* spdmult = esd->GetMultiplicity();
+    if (!spdmult) {
+      AliWarning("No SPD multiplicity");
+    }
+    else { 
+      // Check if we have one or more tracklets 
+      // in the range -1 < eta < 1 to set the INEL>0 
+      // trigger flag. 
+      Int_t n = spdmult->GetNumberOfTracklets();
+      for (Int_t j = 0; j < n; j++) 
+       if(TMath::Abs(spdmult->GetEta(j)) < 1) mESD++;
+    }
+
+    // Reset cache 
+    fData->Reset();
+    Int_t mMC = 0; // Number of particles in |eta|<1
+
+    // Loop over all tracks 
+    Int_t nTracks = mcEvent->GetNumberOfTracks();
+    for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
+      AliMCParticle* particle = 
+       static_cast<AliMCParticle*>(mcEvent->GetTrack(iTr));
+    
+      // Check the returned particle 
+      if (!particle) continue;
+
+      // Check if this charged and a primary 
+      Bool_t isCharged = particle->Charge() != 0;
+      Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
+
+      if (!isCharged || !isPrimary) continue;
+
+      
+      // Fill (eta,phi) of the particle into histograsm for b
+      Double_t eta = particle->Eta();
+      Double_t phi = particle->Phi();
+      
+      fData->Fill(eta, phi);
+      if (TMath::Abs(eta) <= 1) mMC++;
+    }
+    Int_t m = mESD;
+    if (fTrackletRequirement == kMC) m = mMC;
+    fM->Fill(m);
+
+    bool isMcInelGt0 = isMcInel && (mMC > 0);
+    
+    bool hasVertex   = true;
+    if (fVertexRequirement & kMC)  hasVertex = hasVertex && hasMCVtx;
+    if (fVertexRequirement & kESD) hasVertex = hasVertex && hasESDVtx;
+
+    if (isMcInel) {
+      fTriggers->Fill(0);
+      bool triggered = (triggers & AliAODForwardMult::kInel);
+      if (triggered) fTriggers->Fill(1);
+      fInel.AddEvent(triggered, hasVertex, m, fData);
+    }
+    if (isMcInelGt0) {
+      fTriggers->Fill(2);
+      bool triggered = (triggers & AliAODForwardMult::kInelGt0);
+      if (triggered) fTriggers->Fill(3);
+      fInelGt0.AddEvent(triggered, hasVertex, m, fData);
+    }
+    if (isMcNSD) {
+      fTriggers->Fill(4);
+      bool triggered = (triggers & AliAODForwardMult::kNSD);
+      if (triggered) fTriggers->Fill(5);
+      fNSD.AddEvent(triggered, hasVertex, m, fData);
+    }
+    PostData(1, fList);
+  }
+  /** 
+   * End of job processing 
+   */
+  void Terminate(Option_t*)
+  {
+    fList = dynamic_cast<TList*>(GetOutputData(1));
+    if (!fList) {
+      AliError(Form("No output list defined (%p)", GetOutputData(1)));
+      if (GetOutputData(1)) GetOutputData(1)->Print();
+      return;
+    }
+
+
+    TList* output = new TList;
+    output->SetName(GetName());
+    output->SetOwner();
+
+    fVertexMC = static_cast<TH1D*>(fList->FindObject("vertexMC"));
+    fVertexESD = static_cast<TH1D*>(fList->FindObject("vertexESD"));
+    fM         = static_cast<TH1D*>(fList->FindObject("m"));
+    if (fVertexMC) { 
+      TH1D* vtxMC = static_cast<TH1D*>(fVertexMC->Clone("vertexMC"));
+      vtxMC->SetDirectory(0);
+      vtxMC->Scale(1. / vtxMC->GetEntries());
+      output->Add(vtxMC);
+    }
+    if (fVertexESD) { 
+      TH1D* vtxESD = static_cast<TH1D*>(fVertexESD->Clone("vertexESD"));
+      vtxESD->SetDirectory(0);
+      vtxESD->Scale(1. / vtxESD->GetEntries());
+      output->Add(vtxESD);
+    }
+    if (fM) { 
+      TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
+      m->SetDirectory(0);
+      m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
+      m->Scale(1. / m->GetBinContent(1));
+      output->Add(m);
+    }      
+
+    fInel.Finish(fList, output);
+    fInelGt0.Finish(fList, output);
+    fNSD.Finish(fList, output);
+
+    PostData(2, output);
+  }
+    
+protected:
+  //__________________________________________________________________
+  /** 
+   * Structure to hold per trigger type information 
+   */
+  struct TriggerType : public TNamed
+  {
+    //________________________________________________________________
+    /** 
+     * Structure to hold per multiplicity bin information 
+     */
+    struct MBin : public TNamed
+    {
+      TH2D*     fTruth;
+      TH2D*     fTriggered; 
+      TH2D*     fAccepted;
+      TH1I*     fCounts;
+      UShort_t  fLow;
+      UShort_t  fHigh;
+      /** 
+       * Constructor 
+       */
+      MBin() 
+       : fTruth(0), fTriggered(0), fAccepted(0), 
+         fCounts(0), fLow(0), fHigh(1000) {}
+      /** 
+       * Constructor 
+       * 
+       * @param p      Parent list 
+       * @param low    Low cut 
+       * @param high   High cut 
+       * @param eAxis  Eta axis 
+       * @param pAxis  Phi axis 
+       */
+      MBin(TList* p, UShort_t low, UShort_t high, const TH2D* dHist) 
+       : fTruth(0), 
+         fTriggered(0), 
+         fAccepted(0),
+         fCounts(0), 
+         fLow(low), 
+         fHigh(high)
+      {
+       if (low >= high) SetName("all");
+       else             SetName(Form("m%03d_%03d", fLow, fHigh));
+       TList* l = new TList;
+       l->SetOwner();
+       l->SetName(GetName());
+       p->Add(l);
+
+       fTruth = static_cast<TH2D*>(dHist->Clone(("truth")));
+       fTruth->SetTitle("MC truth");
+       fTruth->SetDirectory(0);
+       fTruth->SetZTitle("#sum_i^{N_X} N_{ch}(#eta,#varphi)");
+       fTruth->Sumw2();
+       fTruth->Reset();
+       l->Add(fTruth);
+
+       fTriggered = static_cast<TH2D*>(fTruth->Clone(("triggered")));
+       fTriggered->SetTitle("Triggered");
+       fTriggered->SetDirectory(0);
+       fTriggered->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
+       fTriggered->Sumw2();
+       fTriggered->Reset();
+       l->Add(fTriggered);
+
+       fAccepted = static_cast<TH2D*>(fTruth->Clone(("accepted")));
+       fAccepted->SetTitle("Accepted");
+       fAccepted->SetDirectory(0);
+       fAccepted->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
+       fAccepted->Sumw2();
+       fAccepted->Reset();
+       l->Add(fAccepted);
+       
+       fCounts = new TH1I("counts", "Event counts", 3, -.5, 2.5);
+       fCounts->SetDirectory(0);
+       fCounts->GetXaxis()->SetBinLabel(1, "Truth");
+       fCounts->GetXaxis()->SetBinLabel(2, "Triggered");
+       fCounts->GetXaxis()->SetBinLabel(3, "Accepted");
+       fCounts->SetYTitle("# events");
+       l->Add(fCounts);
+      }
+      /** 
+       * Add event observation
+       * 
+       * @param triggered Whether the event was triggered
+       * @param event     Data for this event 
+       */
+      void AddEvent(Bool_t triggered, Bool_t hasVtx, const TH2D* event) 
+      {
+       fCounts->Fill(0);
+       fTruth->Add(event);
+       if (triggered) { 
+         fCounts->Fill(1);
+         fTriggered->Add(event);
+         if (hasVtx) {
+           fCounts->Fill(2);
+           fAccepted->Add(event);
+         }
+       }
+      }
+      /** 
+       * End of job processing 
+       * 
+       * @param p      Parent list
+       * @param o      Output parent list
+       * @param stack  Stack of histograms
+       * 
+       * @return Trigger efficiency
+       */ 
+     Double_t Finish(const TList* p, TList* o, THStack* stack) 
+      {
+       TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
+       if (!l) { 
+         Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
+         return 0;
+       }
+       fTruth     = static_cast<TH2D*>(l->FindObject("truth"));
+       fTriggered = static_cast<TH2D*>(l->FindObject("triggered"));
+       fAccepted  = static_cast<TH2D*>(l->FindObject("accepted"));
+       fCounts    = static_cast<TH1I*>(l->FindObject("counts"));
+       
+       Int_t    nTruth     = fCounts->GetBinContent(1);
+       Int_t    nTriggered = fCounts->GetBinContent(2);
+       Int_t    nAccepted  = fCounts->GetBinContent(3);
+       Double_t eff        = 0;
+       if (nTruth > 0) eff = double(nTriggered) / nTruth;
+       else if (nTriggered == nTruth) eff = 1;
+
+       if (nTruth > 0)     fTruth->Scale(1. / nTruth);
+       if (nTriggered > 0) fTriggered->Scale(1. / nTriggered);
+       if (nAccepted > 0)  fAccepted->Scale(1. / nAccepted);
+
+       if (fLow >= fHigh) 
+         Info("Finish", "%-6s  [all]  E_X=N_T/N_X=%9d/%-9d=%f "
+              "E_V=N_A/N_T=%9d/%-9d=%f", 
+              p->GetName(), nTriggered, nTruth, eff, nAccepted, nTriggered, 
+              (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
+       else
+         Info("Finish", "%-6s [%2d-%2d] E_X=N_T/N_X=%9d/%-9d=%f "
+              "E_V=N_A/N_T=%9d/%-9d=%f", 
+              p->GetName(), fLow, fHigh, nTriggered, nTruth, eff, 
+              nAccepted, nTriggered, 
+              (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
+       
+       TList* out = new TList;
+       out->SetName(GetName());
+       out->SetOwner();
+       o->Add(out);
+       
+       out->Add(fTruth);
+       out->Add(fTriggered);
+       out->Add(new TParameter<double>("eff", eff));
+       
+       TH2D* bias = static_cast<TH2D*>(fAccepted->Clone("bias"));
+       bias->Divide(fTruth);
+       bias->SetDirectory(0);
+       bias->SetZTitle("Trigger bias (accepted/truth)");
+       out->Add(bias);
+
+       TH1D* truth_px = static_cast<TH1D*>(fTruth->ProjectionX("truth_eta"));
+       truth_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", 
+                              fLow, fHigh));
+       truth_px->Scale(1. / fTruth->GetNbinsY());
+       truth_px->SetDirectory(0);
+       truth_px->SetLineColor(kRed+1);
+       truth_px->SetMarkerColor(kRed+1);
+       truth_px->SetFillColor(kRed+1);
+       truth_px->SetFillStyle(0);
+       out->Add(truth_px);
+
+       TH1D* triggered_px = 
+         static_cast<TH1D*>(fTriggered->ProjectionX("triggered_eta"));
+       triggered_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", 
+                              fLow, fHigh));
+       triggered_px->Scale(1. / fTriggered->GetNbinsY());
+       triggered_px->SetDirectory(0);
+       triggered_px->SetLineColor(kGreen+1);
+       triggered_px->SetMarkerColor(kGreen+1);
+       triggered_px->SetFillColor(kGreen+1);
+       triggered_px->SetFillStyle(0);
+       out->Add(triggered_px);
+
+       TH1D* accepted_px = 
+         static_cast<TH1D*>(fAccepted->ProjectionX("accepted_eta"));
+       accepted_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", 
+                                  fLow, fHigh));
+       accepted_px->Scale(1. / fAccepted->GetNbinsY());
+       accepted_px->SetLineColor(kBlue+1);
+       accepted_px->SetMarkerColor(kBlue+1);
+       accepted_px->SetFillColor(kBlue+1);
+       accepted_px->SetDirectory(0);
+       out->Add(accepted_px);
+
+       THStack* data = new THStack("data", "Data distributions");
+       data->Add(truth_px);
+       data->Add(triggered_px);
+       data->Add(accepted_px);
+       out->Add(data);
+
+#if 0
+       TH1D* bias_px = static_cast<TH1D*>(bias->ProjectionX("bias_eta"));
+       bias_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", 
+                              fLow, fHigh));
+       bias_px->Scale(1. / bias->GetNbinsY());
+#else
+       TH1D* bias_px = static_cast<TH1D*>(accepted_px->Clone("bias_px"));
+       bias_px->Divide(truth_px);
+       bias_px->SetYTitle("Trigger bias (triggered/truth)");
+#endif
+       bias_px->SetDirectory(0);
+       bias_px->SetMarkerStyle(20);
+       bias_px->SetFillStyle(0);
+       bias_px->SetMinimum(0);
+       out->Add(bias_px);
+
+       stack->Add(bias_px);
+
+       return eff;
+      }        
+      ClassDef(MBin,1);
+    };
+
+    /** 
+     * Constructor 
+     */
+    TriggerType() 
+      : TNamed(), 
+       fMask(0),
+       fM(0),
+       fBins(0)
+    {}
+    //--- Constructor ------------------------------------------------
+    /** 
+     * Constructor 
+     * 
+     * @param mask  Trigger mask
+     */
+    TriggerType(UShort_t mask) 
+      : TNamed(AliAODForwardMult::GetTriggerString(mask), ""),
+       fMask(mask), 
+       fM(0), 
+       fBins(0)
+    {
+    }
+    /** 
+     * Create objects 
+     * 
+     * @param list   PArent list
+     * @param mAxis  Multiplicity axis 
+     * @param eAxis  Eta axis
+     * @param pAxis  Phi axis
+     */
+    void CreateObjects(TList* list, 
+                      const TH1D* mHist, 
+                      const TH2D* dHist)
+    {
+      TList* ours  = new TList;
+      ours->SetName(GetName());
+      ours->SetOwner();
+      list->Add(ours);
+
+      fM = static_cast<TH1D*>(mHist->Clone("m"));
+      fM->SetDirectory(0);
+      ours->Add(fM);
+      
+      fBins = new TObjArray;
+      fBins->AddAt(new MBin(ours, 0, 0, dHist), 0);
+      for (Int_t i = 1; i <= fM->GetNbinsX(); i++) { 
+       Double_t low  = fM->GetXaxis()->GetBinLowEdge(i);
+       Double_t high = fM->GetXaxis()->GetBinUpEdge(i);
+
+       fBins->AddAt(new MBin(ours, low, high, dHist), i);
+      }
+    }
+    /** 
+     * Find bin corresponding to m
+     * 
+     * @param m Multiplicity 
+     * 
+     * @return Bin. 
+     */
+    MBin* FindBin(UShort_t m) 
+    { 
+      Int_t b = fM->GetXaxis()->FindBin(m);
+      if (b <= 0) return 0;
+      if (b >= fM->GetNbinsX()+1) b = fM->GetNbinsX();
+      return static_cast<MBin*>(fBins->At(b));
+    }
+    /** 
+     * Add event observation 
+     * 
+     * @param triggered  IF this is triggered
+     * @param m          Multiplicity 
+     * @param data       Observation 
+     */
+    void AddEvent(Bool_t triggered, Bool_t hasVtx, UShort_t m, const TH2D* data)
+    {
+      fM->AddBinContent(1);
+      fM->AddBinContent(TMath::Min(fM->GetNbinsX(), m+2));
+
+      MBin* all = static_cast<MBin*>(fBins->At(0));
+      all->AddEvent(triggered, hasVtx, data);
+      
+      MBin* bin = FindBin(m);
+      bin->AddEvent(triggered, hasVtx, data);      
+    }      
+    /** 
+     * End of job processing 
+     * 
+     * @param p Parent list 
+     * @param o Parent output list 
+     */
+    void Finish(const TList* p, TList* o)
+    {
+      TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
+      if (!l) { 
+       Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
+       return;
+      }
+      
+      TList* ours  = new TList;
+      ours->SetName(GetName());
+      ours->SetOwner();
+      o->Add(ours);
+
+      fM = static_cast<TH1D*>(l->FindObject("m"));
+      if (!fM) { 
+       Warning("Finish", "Didn't find object 'm' in %s", l->GetName());
+       return;
+      }
+      TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
+      m->SetDirectory(0);
+      m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
+      if (m->GetBinContent(1) > 0) 
+       m->Scale(1. / m->GetBinContent(1));
+      ours->Add(m);
+
+      Int_t nBin = fM->GetNbinsX();
+      TH1D* effs = static_cast<TH1D*>(fM->Clone("effs"));
+      effs->SetYTitle("#epsilon_{X}");
+      effs->SetFillColor(kRed+1);
+      effs->SetDirectory(0);
+      effs->SetMinimum(0);
+
+      gStyle->SetPalette(1);
+      Int_t   nCol = gStyle->GetNumberOfColors();
+      THStack* stack = new THStack("biases", "Trigger biases");
+      for (Int_t i = 0; i <= nBin; i++) { 
+       MBin* bin = static_cast<MBin*>(fBins->At(i));
+       effs->SetBinContent(i+1, bin->Finish(l, ours, stack));
+       TH1* h = static_cast<TH1*>(stack->GetHists()->At(i));
+       Int_t col = kBlack;
+       if (i != 0) { 
+         Int_t icol = TMath::Min(nCol-1,int(double(i)/nBin * nCol + .5));
+         col        = gStyle->GetColorPalette(icol);
+       }
+       h->SetMarkerColor(col);
+       h->SetFillColor(col);
+       h->SetLineColor(col);
+      }
+
+      ours->Add(stack);
+      ours->Add(effs);
+    } 
+    UShort_t   fMask;
+    TH1D*      fM;
+    TObjArray* fBins;
+    ClassDef(TriggerType,1);
+  };
+  TriggerType            fInel;
+  TriggerType            fInelGt0;
+  TriggerType            fNSD;
+  AliFMDMCEventInspector fInspector;
+  TList*                 fList;
+  Bool_t                 fFirstEvent;
+  TH2D*                  fData;
+  TH1I*                  fTriggers;
+  UInt_t                 fTrackletRequirement;
+  UInt_t                 fVertexRequirement;
+  TAxis                  fVertexAxis;
+  TH1D*                  fVertexESD;
+  TH1D*                  fVertexMC;
+  TH1D*                  fM;
+  ClassDef(EvaluateTrigger,1);
+};
+#else 
+//====================================================================
+void MakeEvaluateTriggers(const char* esddir, 
+                         Int_t       nEvents    = -1, 
+                         UInt_t      vtx        = 0x1,
+                         UInt_t      trk        = 0x1,
+                         UInt_t      vz         = 10,
+                         Int_t       proof      = 0)
+{
+  // --- Libraries to load -------------------------------------------
+  gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+
+  // --- Check for proof mode, and possibly upload pars --------------
+  if (proof> 0) { 
+    gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C");
+    if (!LoadPars(proof)) { 
+      Error("MakeAOD", "Failed to load PARs");
+      return;
+    }
+  }
+  
+  // --- Our data chain ----------------------------------------------
+  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
+  TChain* chain = MakeChain("ESD", esddir,true);
+  // If 0 or less events is select, choose all 
+  if (nEvents <= 0) nEvents = chain->GetEntries();
+  
+  // --- Set the macro path ------------------------------------------
+  gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
+                          "$ALICE_ROOT/ANALYSIS/macros",
+                          gROOT->GetMacroPath()));
+
+  // --- Creating the manager and handlers ---------------------------
+  AliAnalysisManager *mgr  = new AliAnalysisManager("Triggers", 
+                                                   "Forward multiplicity");
+  AliAnalysisManager::SetCommonFileName("triggers.root");
+
+  // --- ESD input handler -------------------------------------------
+  AliESDInputHandler *esdHandler = new AliESDInputHandler();
+  mgr->SetInputEventHandler(esdHandler);      
+       
+  // --- Monte Carlo handler -----------------------------------------
+  AliMCEventHandler* mcHandler = new AliMCEventHandler();
+  mgr->SetMCtruthEventHandler(mcHandler);
+  mcHandler->SetReadTR(true);    
+
+  // --- Add tasks ---------------------------------------------------
+  // Physics selection 
+  gROOT->Macro(Form("AddTaskPhysicsSelection.C(1,1,0)"));
+
+#if 0
+  // --- Fix up physics selection to give proper A,C, and E triggers -
+  AliInputEventHandler* ih =
+    static_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
+  AliPhysicsSelection* ps = 
+    static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
+  // Ignore trigger class when selecting events.  This mean that we
+  // get offline+(A,C,E) events too
+  ps->SetSkipTriggerClassSelection(true);
+#endif
+
+  // --- compile our code --------------------------------------------
+  gSystem->AddIncludePath("-I${ALICE_ROOT}/PWG2/FORWARD/analysis2 "
+                         "-I${ALICE_ROOT}/ANALYSIS "
+                         "-I${ALICE_ROOT}/include -DBUILD=1");
+  gROOT->LoadMacro("${ALICE_ROOT}/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C++g");
+  
+  // --- Make our object ---------------------------------------------
+  EvaluateTrigger* task = new EvaluateTrigger("triggers");
+  mgr->AddTask(task);
+  task->SetVertexRequirement(vtx);
+  task->SetTrackletRequirement(trk);
+  task->SetVertexAxis(10, -vz, vz);
+
+  // --- create containers for input/output --------------------------
+  AliAnalysisDataContainer *sums = 
+    mgr->CreateContainer("triggerSums", TList::Class(), 
+                        AliAnalysisManager::kOutputContainer, 
+                        AliAnalysisManager::GetCommonFileName());
+  AliAnalysisDataContainer *output = 
+    mgr->CreateContainer("triggerResults", TList::Class(), 
+                        AliAnalysisManager::kParamContainer, 
+                        AliAnalysisManager::GetCommonFileName());
+  
+  // --- connect input/output ----------------------------------------
+  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task, 1, sums);
+  mgr->ConnectOutput(task, 2, output);
+  
+  // --- Run the analysis --------------------------------------------
+  TStopwatch t;
+  if (!mgr->InitAnalysis()) {
+    Error("MakeAOD", "Failed to initialize analysis train!");
+    return;
+  }
+  // Skip terminate if we're so requested and not in Proof or full mode
+  mgr->SetSkipTerminate(false);
+  // Some informative output 
+  mgr->PrintStatus();
+  if (proof) mgr->SetDebugLevel(3);
+  if (mgr->GetDebugLevel() < 1 && !proof) 
+    mgr->SetUseProgressBar(kTRUE,100);
+
+  // Run the train 
+  t.Start();
+  Printf("=== RUNNING ANALYSIS on %9d events ==========================",
+        nEvents);
+  mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
+  t.Stop();
+  t.Print();
+}
+#endif
index 6bffe59d061cde9628bee75cfa59cab65fcc0a1a..47164ffd5dbf9825781c32f7399dc36ca79844ca 100644 (file)
@@ -1150,9 +1150,11 @@ GetData(UShort_t sys,
        mp->Add(AliceCentralInelGt7000());
       }
     }
+#if 0
     else 
       Warning("GetData", "No other results for sys=%d, energy=%d",
              sys, energy);
+#endif
   }
   else if (sys == 2) { 
     // Nothing for PbPb so far 
@@ -1162,7 +1164,7 @@ GetData(UShort_t sys,
       en = Form(", #sqrt{s_{NN}}=%dGeV", energy);
     else 
       en = Form(", #sqrt{s_{NN}}=%f4.2TeV", float(energy)/1000);
-    Warning("GetData", "No other data for PbP b yet");
+    // Warning("GetData", "No other data for PbPb yet");
   }
   else 
     Warning("GetData", "Unknown system %d", sys);