]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
AliMCAuxHandler: A class to handle auxillary input (such as hits or digits)
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 24 Apr 2013 09:31:46 +0000 (09:31 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 24 Apr 2013 09:31:46 +0000 (09:31 +0000)
in MC analysis passes.

AliFMDMCHitHandler: A class to import FMD hits from simulations.

QAPlotter can now do variance rather than min/max in trends.

RunQAMT.sh can now use existing trending_<X>.root files via option -L

Added script RunGridQA.sh to run custom QA passes using MakeQATrain.C

Other minor updates

PWGLF/FORWARD/analysis2/AliFMDMCHitHandler.cxx [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliFMDMCHitHandler.h [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliForwardUtil.cxx
PWGLF/FORWARD/analysis2/AliForwardUtil.h
PWGLF/FORWARD/analysis2/AliMCAuxHandler.cxx [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliMCAuxHandler.h [new file with mode: 0644]
PWGLF/FORWARD/analysis2/qa/QAPlotter.C
PWGLF/FORWARD/analysis2/qa/RunFinalQA.C
PWGLF/FORWARD/analysis2/qa/RunGridQA.sh [new file with mode: 0755]
PWGLF/FORWARD/analysis2/qa/RunQAMT.sh

diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCHitHandler.cxx b/PWGLF/FORWARD/analysis2/AliFMDMCHitHandler.cxx
new file mode 100644 (file)
index 0000000..e1b842c
--- /dev/null
@@ -0,0 +1,350 @@
+#include "AliFMDMCHitHandler.h"
+#include "AliAnalysisManager.h"
+#include <TError.h>
+#include <AliLog.h>
+#include <TFile.h>
+#include <TClonesArray.h>
+#include <TROOT.h>
+#include <AliStack.h>
+#include <AliMCEvent.h>
+
+ClassImp(AliFMDMCHitHandler)
+#if 0 // For Emacs - do not remove
+;
+#endif
+
+//____________________________________________________________________
+AliFMDMCHitHandler::AliFMDMCHitHandler(const char* name,
+                                      const char* what,
+                                      AliMCEventHandler* parent)
+  : AliMCEventHandler(name, what),
+    fParent(parent), 
+    fFile(0),
+    fTree(0),
+    fDir(0),
+    fArray(0),
+    fNEvents(0), 
+    fNEventsPerFile(0),
+    fNEventsInContainer(0),
+    fEvent(0), 
+    fFileNumber(0),
+    fTreeName(""),
+    fFileBase("")
+{
+  // Constructor 
+  // 
+  // Parameters: 
+  //   name The name 
+}
+
+//____________________________________________________________________
+TString*
+AliFMDMCHitHandler::GetParentPath() const
+{
+  if (!fParent) { 
+    AliWarning("No parent");
+    return 0;
+  }
+  return fParent->GetInputPath();
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitHandler::Init(Option_t* opt)
+{
+  // Initialize the input
+  // 
+  // @param opt Options 
+  // 
+  // @return true on success 
+  AliDebugF(10,"AliFMDMCHitHandler::Init(\"%s\")", opt);
+
+  TString option(opt);
+  if (option.EqualTo("proof") || option.EqualTo("local")) return true;
+
+  TString t = "Tree";
+  TString b = "";
+  TClass* cl = gROOT->GetClass(GetTitle());
+  if (cl) { 
+    if (cl->InheritsFrom("AliHit")) { 
+      t += "H";
+      b =  "Hits";
+    }
+    else if (cl->InheritsFrom("AliSDigit")) {
+      t += "S";
+      b =  "SDigits";
+    }
+    else if (cl->InheritsFrom("AliDigit")) {
+      t += "D";
+      b =  "Digits";
+    }
+    else 
+      t = "";
+  }
+  if (!t.IsNull()) fTreeName = t;
+  if (!b.IsNull()) fFileBase = b;
+
+
+  fArray = new TClonesArray(GetTitle());
+  
+  TTree* treeE = fParent->GetTree();
+  if (!treeE) { 
+    AliError("Parent does not have an events tree");
+    return false;
+  }
+
+  // Get number of events in this directory 
+  fNEventsPerFile = -1;
+  fNEvents        = treeE->GetEntries();
+  fEvent          = 0;
+  fFileNumber     = 0;
+
+  if (!OpenFile(fFileNumber)) return false;
+
+  return true;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitHandler::BeginEvent(Long64_t entry)
+{
+  // Called at the beginning of an event 
+  // 
+  // @param entry Entry in tree 
+  // 
+  // @return true on success
+  AliDebugF(10,"AliFMDMCHitHandler::BeginEvent(%lld)", entry);
+
+  if (entry == -1) 
+    fEvent++;
+  else 
+    fEvent = entry;
+
+  if (fEvent >= fNEvents) { 
+    AliWarningF("Event number out of range %d/%d", fEvent, fNEvents);
+    return false;
+  }
+
+  if (fNEventsPerFile < 0) {
+    TTree* treeK = fParent->TreeK();
+    if (!treeK) { 
+      AliError("Parent does not have a kinematics tree");
+      return false;
+    }
+  
+    TFile* fileK = treeK->GetCurrentFile();
+    if (!fileK) { 
+      AliError("Kinematics tree has no associated file");
+      return false;
+    }
+    // Get the number of events per file 
+    fNEventsPerFile = fileK->GetNkeys() - fileK->GetNProcessIDs();
+  }
+  return LoadEvent(fEvent);
+}
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitHandler::Notify(const char* path)
+{
+  // Called when the input file is changed 
+  // 
+  // @param path New path 
+  //
+  // @return true on success
+  AliDebugF(10,"AliFMDMCHitHandler::Notify(\"%s\")", path);
+  return true;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitHandler::FinishEvent()
+{
+  // Called at the end of an event 
+  // 
+  // @return true on success
+  AliDebug(10,"AliFMDMCHitHandler::FinishEvent()");
+  return true;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitHandler::Terminate()
+{
+  // Called at the end of a job 
+  // 
+  // @return true on success 
+  AliDebug(10,"AliFMDMCHitHandler::Terminate()");
+  return true;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitHandler::TerminateIO()
+{
+  // Called at the end of a sub-job
+  // 
+  // @return true on success
+  AliDebug(10,"AliFMDMCHitHandler::TerminateIO()");
+  return true;
+}
+
+//____________________________________________________________________
+void
+AliFMDMCHitHandler::ResetIO()
+{
+  // Reset the I/O
+  // 
+  //
+  AliDebug(10,"AliFMDMCHitHandler::ResetIO()");
+
+  TString* path = GetParentPath();
+  AliDebugF(10,"Got parent path %s", path ? path->Data() : "null");
+
+  if (fFile) { 
+    delete fFile;
+    fFile = 0;
+  }
+}
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitHandler::OpenFile(Int_t fileNo)
+{
+  TString* path = GetParentPath();
+  AliDebugF(10,"Got parent path %s", path ? path->Data() : "null");
+  if (!path) return false;
+
+  TString ext("");
+  if (fileNo > 0) ext = TString::Format("%d", fileNo);
+
+  TString w(GetTitle());
+  if (w.EndsWith("s")) w.Chop();
+    
+  TString fn = TString::Format("%s%s.%s%s.root", 
+                              path->Data(), GetName(), 
+                              fFileBase.Data(), ext.Data());
+  Info("Init", "Opening %s", fn.Data());
+  fFile = TFile::Open(fn, "READ");
+  if (!fFile) { 
+    AliErrorF("Failed to open %s", fn.Data());
+    return false;
+  }
+
+  return true;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitHandler::LoadEvent(Int_t iev)
+{
+  // Load an event 
+  // 
+  // @param iev Event number 
+  // 
+  // @return true on success 
+  AliDebugF(10,"AliFMDMCHitHandler::LoadEvent(%d)", iev);
+
+  Int_t iNew = iev / fNEventsPerFile;
+  if (iNew != fFileNumber) { 
+    fFileNumber = iNew;
+    if (!OpenFile(fFileNumber)) return false;
+  }
+  if (!fFile) return false;
+
+  TString folder = TString::Format("Event%d", iev);
+  fFile->GetObject(folder, fDir);
+  if (!fDir) { 
+    AliWarningF("Folder %s not found in file", folder.Data());
+    return false;
+  }
+
+  fDir->GetObject(fTreeName, fTree);
+  if (!fTree) { 
+    AliWarningF("Folder %s does not contain the %s tree %s", 
+               folder.Data(), GetTitle(), fTreeName.Data());
+    return false;
+  }
+
+  fTree->SetBranchAddress(GetName(), &fArray);
+  return true;
+}
+
+
+//____________________________________________________________________
+AliFMDMCHitHandler*
+AliFMDMCHitHandler::Create(const char* name, const char* what)
+{
+  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) { 
+    ::Error("AliFMDMCHitHandler::Create", "No analysis manager");
+    return 0;
+  }
+  
+  AliVEventHandler* vmc = mgr->GetMCtruthEventHandler();
+  if (!vmc) { 
+    ::Error("AliFMDMCHitHandler::Create", "No MC truth handler");
+    return 0;
+  }
+  
+  AliMCEventHandler* mc = dynamic_cast<AliMCEventHandler*>(vmc);
+  if (!mc) { 
+    ::Error("AliFMDMCHitHandler::Create", 
+           "MC truth handler not a AliMCEventHandler, but %s", 
+           vmc->ClassName());
+    return 0;
+  }
+
+  AliFMDMCHitHandler* ret = new AliFMDMCHitHandler(name, what, mc);
+  mc->AddSubsidiaryHandler(ret);
+  
+  return ret;
+}
+
+//____________________________________________________________________
+TClonesArray*
+AliFMDMCHitHandler::GetParticleArray(AliFMDMCHitHandler* handler, 
+                                    Int_t particle)
+{
+  if (!handler) { 
+    ::Error("AliFMDMCHitHandler::GetArray", "No handler passed");
+    return 0;
+  }
+
+  AliMCEventHandler* mc = handler->GetParent();
+  if (!mc) {
+    ::Error("AliFMDMCHitHandler::GetArray", "Handler has no parent");
+    return 0;
+  }
+  
+  AliMCEvent* event = mc->MCEvent();
+  if (!event) { 
+    ::Error("AliFMDMCHitHandler::GetArray", "No MC event");
+    return 0;
+  }
+  
+  AliStack* stack = event->Stack();
+  if (!event) { 
+    ::Error("AliFMDMCHitHandler::GetArray", "Event has no stack");
+    return 0;
+  }
+  
+  TTree* tree = handler->GetTree();
+  if (!tree) {
+    ::Error("AliFMDMCHitHandler::GetArray", "Handler has no tree");
+    return 0;
+  }
+    
+  Int_t treeIdx = stack->TreeKEntry(particle);
+  if (treeIdx < 0 || treeIdx >= tree->GetEntries()) { 
+    ::Error("AliFMDMCHitHandler::GetArray", 
+           "Index %d of %d out of bounds [0,%lld]", treeIdx, particle, 
+           tree->GetEntries()-1);
+    return 0;
+  }
+  
+  tree->GetEntry(treeIdx);
+  
+  return handler->GetArray();
+}
+
+  
+    
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCHitHandler.h b/PWGLF/FORWARD/analysis2/AliFMDMCHitHandler.h
new file mode 100644 (file)
index 0000000..2aae177
--- /dev/null
@@ -0,0 +1,185 @@
+/**
+ * @file   AliFMDMCHitHandler.h
+ * @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
+ * @date   Thu Feb  7 12:04:02 2013
+ * 
+ * @brief  
+ * 
+ * 
+ */
+#ifndef ALIFMDMCHITHANDLER_H
+#define ALIFMDMCHITHANDLER_H
+#include <AliMCEventHandler.h>
+class TFile;
+class TTree;
+
+/**
+ * This class defines an input handler for simulated data which will
+ * connect the FMD Hit tree.  It is intended to be added to the global
+ * MC input handler using AliMCEventHandler::AddSubsidiaryHandler
+ * 
+ */
+class AliFMDMCHitHandler : public AliMCEventHandler
+{
+public:
+  /** 
+   * Constructor 
+   * 
+   * @param name Name 
+   */
+  AliFMDMCHitHandler(const char* name="FMD",
+                    const char* clsName="AliFMDHit",
+                    AliMCEventHandler* parent=0);
+
+  virtual ~AliFMDMCHitHandler() {}
+  /** 
+   * @{
+   * @name Interface member functions 
+   */
+  /** 
+   * Intialize 
+   * 
+   * @param t Not used
+   * @param o Not used
+   * 
+   * @return always true 
+   */  
+  virtual Bool_t Init(TTree* t,Option_t* o) { return AliMCEventHandler::Init(t,o); }
+  /** 
+   * Initialize the input
+   * 
+   * @param opt Options 
+   * 
+   * @return true on success 
+   */
+  virtual Bool_t Init(Option_t* opt);
+  /**
+   * Called at the beginning of an event 
+   * 
+   * @param entry Entry in tree 
+   * 
+   * @return true on success
+   */
+  virtual Bool_t BeginEvent(Long64_t entry);
+  /** 
+   * Called when the input file is changed 
+   * 
+   * @param path New path 
+   *
+   * @return true on success
+   */
+  virtual Bool_t Notify() { return AliMCEventHandler::Notify(); }
+  /** 
+   * Called when the input file is changed 
+   * 
+   * @param path New path 
+   *
+   * @return true on success
+   */
+  virtual Bool_t Notify(const char* path);
+  /** 
+   * Called at the end of an event 
+   * 
+   * @return true on success
+   */
+  virtual Bool_t FinishEvent();
+  /** 
+   * Called at the end of a job 
+   * 
+   * @return true on success 
+   */
+  virtual Bool_t Terminate();
+  /** 
+   * Called at the end of a sub-job
+   * 
+   * @return true on success
+   */  
+  virtual Bool_t TerminateIO();
+  /** 
+   * Reset the I/O
+   * 
+   */
+  virtual void ResetIO();
+  /** 
+   * Load an event 
+   * 
+   * @param iev Event number 
+   * 
+   * @return true on success 
+   */
+  virtual Bool_t LoadEvent(Int_t iev);
+  /** 
+   * Set the number of events in the container 
+   * 
+   * @param nev Number of events 
+   */
+  virtual void  SetNumberOfEventsInContainer(Int_t nev) {
+    this->fNEventsInContainer = nev;}
+  virtual Bool_t OpenFile(Int_t ev);
+  /* @} */
+
+  AliMCEventHandler* GetParent() { return fParent; }
+  /** 
+   * Get the tree 
+   * 
+   * @return The connected hits tree
+   */
+  virtual TTree*  GetTree() const { return fTree;}  
+  TClonesArray*      GetArray() { return fArray; }
+  /** 
+   * Static member function to create and attach this handler 
+   * 
+   * @param name Name of the handler 
+   * 
+   * @return Newly allocated handler or null
+   */
+  static AliFMDMCHitHandler* Create(const char* name="FMD",
+                                   const char* what="Hits");
+  static TClonesArray* GetParticleArray(AliFMDMCHitHandler* handler, 
+                                       Int_t particle);
+protected:
+  AliFMDMCHitHandler(const AliFMDMCHitHandler& o)
+    : AliMCEventHandler(), 
+      fParent(o.fParent),
+      fFile(0), 
+      fTree(0), 
+      fDir(0),
+      fArray(0),
+      fNEvents(0), 
+      fNEventsPerFile(0),
+      fNEventsInContainer(0),
+      fEvent(0), 
+      fFileNumber(0), 
+      fTreeName(""),
+      fFileBase("")
+  {}
+  AliFMDMCHitHandler& operator=(const AliFMDMCHitHandler& o)
+  {
+    if (&o == this) return *this;
+    // AliMCEventHandler::operator=(o);
+    fParent = o.fParent;
+    fFile   = o.fFile;
+    fTree   = o.fTree;
+    return *this;
+  }
+  TString* GetParentPath() const;
+  AliMCEventHandler* fParent; // Parent MC handler 
+  TFile*             fFile;                //!
+  TTree*             fTree;                //!
+  TDirectory*        fDir;                 //!
+  TClonesArray*      fArray;               //!
+  Int_t              fNEvents;             //!
+  Int_t              fNEventsPerFile;      //!
+  Int_t              fNEventsInContainer;  //!
+  Int_t              fEvent;               //!
+  Int_t              fFileNumber;          //!
+  TString            fTreeName;            //!
+  TString            fFileBase;            //! 
+  ClassDef(AliFMDMCHitHandler,1); // Connect FMD hits tree
+};
+
+#endif
+// Local Variables:
+//  mode: C++
+// End:
+
index 0e6e306d9f2ffbd83f3f596d453494ea3987c698..9836dc5954a481827e0b01e3378a94e614c59299 100644 (file)
@@ -469,6 +469,22 @@ namespace {
     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigmaN);
   }
 
+  Double_t landauGausComposite(Double_t* xp, Double_t* pp)
+  {
+    Double_t x           = xp[0];
+    Double_t cP          = pp[AliForwardUtil::ELossFitter::kC];
+    Double_t deltaP      = pp[AliForwardUtil::ELossFitter::kDelta];
+    Double_t xiP         = pp[AliForwardUtil::ELossFitter::kXi];
+    Double_t sigmaP      = pp[AliForwardUtil::ELossFitter::kSigma];
+    Double_t cS          = pp[AliForwardUtil::ELossFitter::kSigma+1];
+    Double_t deltaS      = pp[AliForwardUtil::ELossFitter::kSigma+2];
+    Double_t xiS         = pp[AliForwardUtil::ELossFitter::kSigma+3];
+    Double_t sigmaS      = pp[AliForwardUtil::ELossFitter::kSigma+4];
+
+    return (cP * AliForwardUtil::LandauGaus(x,deltaP,xiP,sigmaP,0) + 
+           cS * AliForwardUtil::LandauGaus(x,deltaS,xiS,sigmaS,0));
+  }
+    
   // 
   // Utility function to use in TF1 defintition 
   //
@@ -926,34 +942,45 @@ AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
   dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
   
   // Get the bin with maximum 
-  Int_t    maxBin = dist->GetMaximumBin();
-  Double_t maxE   = dist->GetBinLowEdge(maxBin);
+  Int_t    peakBin = dist->GetMaximumBin();
+  Double_t peakE   = dist->GetBinLowEdge(peakBin);
   
   // Get the low edge 
-  dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
-  Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
+  dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
+  Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
-  Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
+  Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
 
+  Int_t    minEb = dist->GetXaxis()->FindBin(minE);
+  Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
+  Double_t intg  = dist->Integral(minEb, maxEb);
+  if (intg <= 0) {
+    ::Warning("Fit1Particle", 
+             "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0", 
+             dist->GetName(), minE, maxE, minEb, maxEb, intg);
+    return 0;
+  }
+    
   // Restore the range 
   dist->GetXaxis()->SetRangeUser(0, fMaxRange);
   
   // Define the function to fit 
-  TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
+  TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxE,kSigmaN+1);
 
   // Set initial guesses, parameter names, and limits  
-  landau1->SetParameters(1,0.5,0.07,0.1,sigman);
+  landau1->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
   landau1->SetNpx(500);
   landau1->SetParLimits(kDelta, minE, fMaxRange);
   landau1->SetParLimits(kXi,    0.00, fMaxRange);
-  landau1->SetParLimits(kSigma, 0.01, 0.1);
+  landau1->SetParLimits(kSigma, 1e-5, fMaxRange);
   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
   else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
 
   // Do the fit, getting the result object 
-  TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
-  landau1->SetRange(minE, fMaxRange);
+  ::Info("Fit1Particle", "Fitting in the range %f,%f", minE, maxE);
+  TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxE);
+  // landau1->SetRange(minE, fMaxRange);
   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
   fFunctions.AddAtAndExpand(landau1, 0);
 
@@ -998,20 +1025,30 @@ AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
   Double_t minE    = f->GetXmin();
 
+  Int_t    minEb = dist->GetXaxis()->FindBin(minE);
+  Int_t    maxEb = dist->GetXaxis()->FindBin(maxEi);
+  Double_t intg  = dist->Integral(minEb, maxEb);
+  if (intg <= 0) {
+    ::Warning("FitNParticle",
+             "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0", 
+             dist->GetName(), minE, maxEi, minEb, maxEb, intg);
+    return 0;
+  }
+
   // Array of weights 
   TArrayD a(n-1);
   for (UShort_t i = 2; i <= n; i++) 
     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
   // Make the fit function 
-  TF1* landaun     = MakeNLandauGaus(r->Parameter(kC),
-                                    r->Parameter(kDelta),
-                                    r->Parameter(kXi),
-                                    r->Parameter(kSigma),
-                                    r->Parameter(kSigmaN),
-                                    n,a.fArray,minE,maxEi);
+  TF1* landaun = MakeNLandauGaus(r->Parameter(kC),
+                                r->Parameter(kDelta),
+                                r->Parameter(kXi),
+                                r->Parameter(kSigma),
+                                r->Parameter(kSigmaN),
+                                n, a.fArray, minE, maxEi);
   landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
   landaun->SetParLimits(kXi,     0.00, fMaxRange);       // xi
-  landaun->SetParLimits(kSigma,  0.01, 1);            // sigma
+  landaun->SetParLimits(kSigma,  1e-5, fMaxRange);       // sigma
   // Check if we're using the noise sigma 
   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
   else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
@@ -1022,14 +1059,133 @@ AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
   }
 
   // Do the fit 
+  ::Info("Fit1Particle", "Fitting in the range %f,%f", minE, maxEi);
   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
   
-  landaun->SetRange(minE, fMaxRange);
+  // landaun->SetRange(minE, fMaxRange);
   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
   fFunctions.AddAtAndExpand(landaun, n-1);
   
   return landaun;
 }  
+//____________________________________________________________________
+TF1*
+AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
+{
+  // 
+  // Fit a composite particle signal to the passed energy loss
+  // distribution
+  // 
+  // Parameters:
+  //    dist    Data to fit the function to 
+  //    sigman If larger than zero, the initial guess of the
+  //               detector induced noise. If zero or less, then this 
+  //               parameter is ignored in the fit (fixed at 0)
+  // 
+  // Return:
+  //    The function fitted to the data 
+  //
+
+  // Find the fit range 
+  dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
+  
+  // Get the bin with maximum 
+  Int_t    peakBin = dist->GetMaximumBin();
+  Double_t peakE   = dist->GetBinLowEdge(peakBin);
+  
+  // Get the low edge 
+  dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
+  Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
+  Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
+  Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
+
+  // Get the range in bins and the integral of that range 
+  Int_t    minEb = dist->GetXaxis()->FindBin(minE);
+  Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
+  Double_t intg  = dist->Integral(minEb, maxEb);
+  if (intg <= 0) {
+    ::Warning("Fit1Particle", 
+             "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0", 
+             dist->GetName(), minE, maxE, minEb, maxEb, intg);
+    return 0;
+  }
+    
+  // Restore the range 
+  dist->GetXaxis()->SetRangeUser(0, fMaxRange);
+  
+  // Define the function to fit 
+  TF1* seed = new TF1("landauSeed", landauGaus1, minE,maxE,kSigmaN+1);
+
+  // Set initial guesses, parameter names, and limits  
+  seed->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
+  seed->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
+  seed->SetNpx(500);
+  seed->SetParLimits(kDelta, minE, fMaxRange);
+  seed->SetParLimits(kXi,    0.00, fMaxRange);
+  seed->SetParLimits(kSigma, 1e-5, fMaxRange);
+  if (sigman <= 0)  seed->FixParameter(kSigmaN, 0);
+  else              seed->SetParLimits(kSigmaN, 0, fMaxRange);
+
+  // Do the fit, getting the result object 
+  ::Info("FitComposite", "Fitting seed in the range %f,%f", minE, maxE);
+  /* TFitResultPtr r = */ dist->Fit(seed, "RNQS", "", minE, maxE);
+
+  maxE = dist->GetXaxis()->GetXmax();
+  TF1* comp = new TF1("composite", landauGausComposite, 
+                     minE, maxE, kSigma+1+4);
+  comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
+                   "C#prime", "#Delta_{p}#prime", "#xi#prime", "#sigma#prim");
+  comp->SetParameters(0.8 * seed->GetParameter(kC),  // 0 Primary weight 
+                     seed->GetParameter(kDelta),    // 1 Primary Delta
+                     seed->GetParameter(kDelta)/10, // 2 primary Xi
+                     seed->GetParameter(kDelta)/5,  // 3 primary sigma
+                     1.20 * seed->GetParameter(kC), // 5 Secondary weight
+                     seed->GetParameter(kDelta),    // 6 secondary Delta
+                     seed->GetParameter(kXi),       // 7 secondary Xi
+                     seed->GetParameter(kSigma));   // 8 secondary sigma
+                     
+  // comp->SetParLimits(kC,       minE, fMaxRange); // C
+  comp->SetParLimits(kDelta,      minE, fMaxRange); // Delta
+  comp->SetParLimits(kXi,         0.00, fMaxRange); // Xi 
+  comp->SetParLimits(kSigma,      1e-5, fMaxRange); // Sigma
+  // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
+  comp->SetParLimits(kSigma+2,    minE/10, fMaxRange); // Delta
+  comp->SetParLimits(kSigma+3,    0.00,    fMaxRange); // Xi 
+  comp->SetParLimits(kSigma+4,    1e-6,    fMaxRange); // Sigma
+  comp->SetLineColor(kRed+1);
+  comp->SetLineWidth(3);
+  
+  // Do the fit, getting the result object 
+  ::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
+  /* TFitResultPtr r = */ dist->Fit(comp, "RNQS", "", minE, maxE);
+
+#if 0
+  TF1* part1 = static_cast<TF1*>(seed->Clone("part1"));
+  part1->SetLineColor(kGreen+1);
+  part1->SetLineWidth(4);
+  part1->SetRange(minE, maxE);
+  part1->SetParameters(comp->GetParameter(0), // C 
+                      comp->GetParameter(1), // Delta
+                      comp->GetParameter(2), // Xi
+                      comp->GetParameter(3), // sigma
+                      0);
+  part1->Save(minE,maxE,0,0,0,0);
+  dist->GetListOfFunctions()->Add(part1);
+
+  TF1* part2 = static_cast<TF1*>(seed->Clone("part2"));
+  part2->SetLineColor(kBlue+1);
+  part2->SetLineWidth(4);
+  part2->SetRange(minE, maxE);
+  part2->SetParameters(comp->GetParameter(4), // C 
+                      comp->GetParameter(5), // Delta
+                      comp->GetParameter(6), // Xi
+                      comp->GetParameter(7), // sigma
+                      0);
+  part2->Save(minE,maxE,0,0,0,0);
+  dist->GetListOfFunctions()->Add(part2);
+#endif
+  return comp;
+}
 
 //====================================================================
 AliForwardUtil::Histos::~Histos()
index d286fe41c7f82d6b9097547c29134f42da8e8249..d780ac430451ecc6b1753e9b26b6bffa84bbcaba 100644 (file)
@@ -176,6 +176,21 @@ public:
    * @return Short integer value of magnetic field in kG 
    */
   static Short_t ParseMagneticField(Float_t field);
+  /** 
+   * Get a string representation of the magnetic field
+   * 
+   * @param field Magnetic field in kG
+   * 
+   * @return String representation of the magnetic field
+   */
+  static const char* MagneticFieldString(Short_t field);
+  /* @} */
+
+  //==================================================================
+  /** 
+   * @{ 
+   * @name Recalculate @f$\eta@f$, @f$\phi@f$, etc. 
+   */
   /** 
    * Get the radius of a strip. 
    * 
@@ -207,17 +222,13 @@ public:
    */  
   static Double_t GetPhiFromStrip(Char_t ring, UShort_t strip, 
                                  Double_t phi, Double_t xvtx, Double_t yvtx);
-  /** 
-   * Get a string representation of the magnetic field
-   * 
-   * @param field Magnetic field in kG
-   * 
-   * @return String representation of the magnetic field
-   */
-   static const char* MagneticFieldString(Short_t field);
   /* @} */
 
-  //__________________________________________________________________
+  //==================================================================
+  /** 
+   * @{ 
+   * @name Manager related tasks 
+   */
   /** 
    * Get the AOD event - either from the input (AOD analysis) or the
    * output (ESD analysis)
@@ -245,8 +256,9 @@ public:
    * @return true if the needed task was found 
    */
   static Bool_t CheckForTask(const char* clsOrName, Bool_t cls=true);
+  /* @} */
 
-  //__________________________________________________________________
+  //==================================================================
   /** 
    * @{ 
    * @name Member functions to store and retrieve analysis parameters 
@@ -263,6 +275,7 @@ public:
   static void GetParameter(TObject* o, ULong_t& value);
   /* @} */
 
+  //==================================================================
   /** 
    * @{ 
    * @name Energy stragling functions 
@@ -531,6 +544,18 @@ public:
      * @return The function fitted to the data 
      */
     TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
+    /** 
+     * Fit a composite distribution of energy loss from both primaries
+     * and secondaries
+     * 
+     * @param dist   Distribution 
+     * @param sigman If larger than zero, the initial guess of the
+     *                detector included noise.  If zero or less this
+     *                parameter is fixed to 0.
+     * 
+     * @return Function fitted to the data 
+     */
+    TF1* FitComposite(TH1* dist, Double_t sigman);
     /**
      * Get Lower cut on data 
      *
@@ -748,7 +773,7 @@ public:
      */
     TH1* GetOutputHist(const TList* d, const char* name) const;
     /** 
-     * 
+     * Get the colour of this ring 
      * 
      * 
      * @return 
diff --git a/PWGLF/FORWARD/analysis2/AliMCAuxHandler.cxx b/PWGLF/FORWARD/analysis2/AliMCAuxHandler.cxx
new file mode 100644 (file)
index 0000000..e6ba89a
--- /dev/null
@@ -0,0 +1,377 @@
+#include "AliMCAuxHandler.h"
+#include "AliAnalysisManager.h"
+#include <TError.h>
+#include <AliLog.h>
+#include <TFile.h>
+#include <TClonesArray.h>
+#include <TROOT.h>
+#include <AliStack.h>
+#include <AliMCEvent.h>
+
+ClassImp(AliMCAuxHandler)
+#if 0 // For Emacs - do not remove
+;
+#endif
+
+//____________________________________________________________________
+AliMCAuxHandler::AliMCAuxHandler(const char* name,
+                                      const char* what,
+                                      AliMCEventHandler* parent)
+  : AliMCEventHandler(name, what),
+    fParent(parent), 
+    fFile(0),
+    fTree(0),
+    fDir(0),
+    fArray(0),
+    fNEvents(0), 
+    fNEventsPerFile(0),
+    fNEventsInContainer(0),
+    fEvent(0), 
+    fFileNumber(0),
+    fTreeName(""),
+    fFileBase("")
+{
+  // Constructor 
+  // 
+  // Parameters: 
+  //   name The name 
+}
+
+//____________________________________________________________________
+TString*
+AliMCAuxHandler::GetParentPath() const
+{
+  if (!fParent) { 
+    AliWarning("No parent");
+    return 0;
+  }
+  return fParent->GetInputPath();
+}
+
+//____________________________________________________________________
+Bool_t
+AliMCAuxHandler::Init(Option_t* opt)
+{
+  // Initialize the input
+  // 
+  // @param opt Options 
+  // 
+  // @return true on success 
+  AliDebugF(10,"AliMCAuxHandler::Init(\"%s\")", opt);
+
+  TString option(opt);
+  if (option.EqualTo("proof") || option.EqualTo("local")) return true;
+
+  TString t = "Tree";
+  TString b = "";
+  TClass* cl = gROOT->GetClass(GetTitle());
+  if (cl) { 
+    if (cl->InheritsFrom("AliHit")) { 
+      t += "H";
+      b =  "Hits";
+    }
+    else if (cl->InheritsFrom("AliSDigit")) {
+      t += "S";
+      b =  "SDigits";
+    }
+    else if (cl->InheritsFrom("AliDigit")) {
+      t += "D";
+      b =  "Digits";
+    }
+    else 
+      t = "";
+  }
+  if (!t.IsNull()) fTreeName = t;
+  if (!b.IsNull()) fFileBase = b;
+
+
+  fArray = new TClonesArray(GetTitle());
+  
+  TTree* treeE = fParent->GetTree();
+  if (!treeE) { 
+    AliError("Parent does not have an events tree");
+    return false;
+  }
+
+  // Get number of events in this directory 
+  fNEventsPerFile = -1;
+  fNEvents        = treeE->GetEntries();
+  fEvent          = 0;
+  fFileNumber     = 0;
+
+  if (!OpenFile(fFileNumber)) return false;
+
+  return true;
+}
+//____________________________________________________________________
+Bool_t
+AliMCAuxHandler::BeginEvent(Long64_t entry)
+{
+  // Called at the beginning of an event 
+  // 
+  // @param entry Entry in tree 
+  // 
+  // @return true on success
+  AliDebugF(10,"AliMCAuxHandler::BeginEvent(%lld)", entry);
+
+  if (entry == -1) 
+    fEvent++;
+  else 
+    fEvent = entry;
+
+  if (fEvent >= fNEvents) { 
+    AliWarningF("Event number out of range %d/%d", fEvent, fNEvents);
+    return false;
+  }
+
+  if (fNEventsPerFile < 0) {
+    TTree* treeK = fParent->TreeK();
+    if (!treeK) { 
+      AliError("Parent does not have a kinematics tree");
+      return false;
+    }
+  
+    TFile* fileK = treeK->GetCurrentFile();
+    if (!fileK) { 
+      AliError("Kinematics tree has no associated file");
+      return false;
+    }
+    // Get the number of events per file 
+    fNEventsPerFile = fileK->GetNkeys() - fileK->GetNProcessIDs();
+  }
+  return LoadEvent(fEvent);
+}
+//____________________________________________________________________
+Bool_t
+AliMCAuxHandler::Notify(const char* path)
+{
+  // Called when the input file is changed 
+  // 
+  // @param path New path 
+  //
+  // @return true on success
+  AliDebugF(10,"AliMCAuxHandler::Notify(\"%s\")", path);
+  return true;
+}
+//____________________________________________________________________
+Bool_t
+AliMCAuxHandler::FinishEvent()
+{
+  // Called at the end of an event 
+  // 
+  // @return true on success
+  AliDebug(10,"AliMCAuxHandler::FinishEvent()");
+  return true;
+}
+//____________________________________________________________________
+Bool_t
+AliMCAuxHandler::Terminate()
+{
+  // Called at the end of a job 
+  // 
+  // @return true on success 
+  AliDebug(10,"AliMCAuxHandler::Terminate()");
+  return true;
+}
+//____________________________________________________________________
+Bool_t
+AliMCAuxHandler::TerminateIO()
+{
+  // Called at the end of a sub-job
+  // 
+  // @return true on success
+  AliDebug(10,"AliMCAuxHandler::TerminateIO()");
+  return true;
+}
+
+//____________________________________________________________________
+void
+AliMCAuxHandler::ResetIO()
+{
+  // Reset the I/O
+  // 
+  //
+  AliDebug(10,"AliMCAuxHandler::ResetIO()");
+
+  TString* path = GetParentPath();
+  AliDebugF(10,"Got parent path %s", path ? path->Data() : "null");
+
+  if (fFile) { 
+    delete fFile;
+    fFile = 0;
+  }
+}
+//____________________________________________________________________
+Bool_t
+AliMCAuxHandler::OpenFile(Int_t fileNo)
+{
+  TString* path = GetParentPath();
+  AliDebugF(10,"Got parent path %s", path ? path->Data() : "null");
+  if (!path) return false;
+
+  TString ext("");
+  if (fileNo > 0) ext = TString::Format("%d", fileNo);
+
+  TString w(GetTitle());
+  if (w.EndsWith("s")) w.Chop();
+    
+  TString fn = TString::Format("%s%s.%s%s.root", 
+                              path->Data(), GetName(), 
+                              fFileBase.Data(), ext.Data());
+  Info("Init", "Opening %s", fn.Data());
+  fFile = TFile::Open(fn, "READ");
+  if (!fFile) { 
+    AliErrorF("Failed to open %s", fn.Data());
+    return false;
+  }
+
+  return true;
+}
+
+//____________________________________________________________________
+Bool_t
+AliMCAuxHandler::LoadEvent(Int_t iev)
+{
+  // Load an event 
+  // 
+  // @param iev Event number 
+  // 
+  // @return true on success 
+  AliDebugF(10,"AliMCAuxHandler::LoadEvent(%d)", iev);
+
+  Int_t iNew = iev / fNEventsPerFile;
+  if (iNew != fFileNumber) { 
+    fFileNumber = iNew;
+    if (!OpenFile(fFileNumber)) return false;
+  }
+  if (!fFile) return false;
+
+  TString folder = TString::Format("Event%d", iev);
+  fFile->GetObject(folder, fDir);
+  if (!fDir) { 
+    AliWarningF("Folder %s not found in file", folder.Data());
+    return false;
+  }
+
+  fDir->GetObject(fTreeName, fTree);
+  if (!fTree) { 
+    AliWarningF("Folder %s does not contain the %s tree %s", 
+               folder.Data(), GetTitle(), fTreeName.Data());
+    return false;
+  }
+
+  fTree->SetBranchAddress(GetName(), &fArray);
+  return true;
+}
+
+//____________________________________________________________________
+Int_t
+AliMCAuxHandler::GetNEntry() const
+{
+  if (!fTree) return 0;
+  return fTree->GetEntries();
+}
+
+//____________________________________________________________________
+TClonesArray*
+AliMCAuxHandler::GetEntryArray(Int_t entry)
+{
+  if (!fTree) return 0;
+  if (!fArray) return 0;
+  if (entry < 0 || entry >= fTree->GetEntries()) { 
+    AliErrorF("Entry # %d out of bounds [0,%lld]", 
+             entry, fTree->GetEntries());
+    return 0;
+  }
+  fArray->Clear();
+  
+  if (fTree->GetEntry(entry) <= 0) return 0;
+
+  return fArray;
+}
+  
+
+//____________________________________________________________________
+AliMCAuxHandler*
+AliMCAuxHandler::Create(const char* name, const char* what)
+{
+  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) { 
+    ::Error("AliMCAuxHandler::Create", "No analysis manager");
+    return 0;
+  }
+  
+  AliVEventHandler* vmc = mgr->GetMCtruthEventHandler();
+  if (!vmc) { 
+    ::Error("AliMCAuxHandler::Create", "No MC truth handler");
+    return 0;
+  }
+  
+  AliMCEventHandler* mc = dynamic_cast<AliMCEventHandler*>(vmc);
+  if (!mc) { 
+    ::Error("AliMCAuxHandler::Create", 
+           "MC truth handler not a AliMCEventHandler, but %s", 
+           vmc->ClassName());
+    return 0;
+  }
+
+  AliMCAuxHandler* ret = new AliMCAuxHandler(name, what, mc);
+  mc->AddSubsidiaryHandler(ret);
+  
+  return ret;
+}
+
+//____________________________________________________________________
+TClonesArray*
+AliMCAuxHandler::GetParticleArray(AliMCAuxHandler* handler, 
+                                    Int_t particle)
+{
+  if (!handler) { 
+    ::Error("AliMCAuxHandler::GetArray", "No handler passed");
+    return 0;
+  }
+
+  AliMCEventHandler* mc = handler->GetParent();
+  if (!mc) {
+    ::Error("AliMCAuxHandler::GetArray", "Handler has no parent");
+    return 0;
+  }
+  
+  AliMCEvent* event = mc->MCEvent();
+  if (!event) { 
+    ::Error("AliMCAuxHandler::GetArray", "No MC event");
+    return 0;
+  }
+  
+  AliStack* stack = event->Stack();
+  if (!event) { 
+    ::Error("AliMCAuxHandler::GetArray", "Event has no stack");
+    return 0;
+  }
+
+  handler->GetArray()->Clear();
+  TTree* tree = handler->GetTree();
+  if (!tree) {
+    ::Error("AliMCAuxHandler::GetArray", "Handler has no tree");
+    return 0;
+  }
+    
+  Int_t treeIdx = stack->TreeKEntry(particle);
+  if (treeIdx < 0 || treeIdx >= tree->GetEntries()) { 
+    ::Error("AliMCAuxHandler::GetArray", 
+           "Index %d of %d out of bounds [0,%lld]", treeIdx, particle, 
+           tree->GetEntries()-1);
+    return 0;
+  }
+  
+  tree->GetEntry(treeIdx);
+  
+  return handler->GetArray();
+}
+
+  
+    
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/PWGLF/FORWARD/analysis2/AliMCAuxHandler.h b/PWGLF/FORWARD/analysis2/AliMCAuxHandler.h
new file mode 100644 (file)
index 0000000..36a4402
--- /dev/null
@@ -0,0 +1,187 @@
+/**
+ * @file   AliMCAuxHandler.h
+ * @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
+ * @date   Thu Feb  7 12:04:02 2013
+ * 
+ * @brief  
+ * 
+ * 
+ */
+#ifndef ALIFMDMCHITHANDLER_H
+#define ALIFMDMCHITHANDLER_H
+#include <AliMCEventHandler.h>
+class TFile;
+class TTree;
+
+/**
+ * This class defines an input handler for simulated data which will
+ * connect the FMD Hit tree.  It is intended to be added to the global
+ * MC input handler using AliMCEventHandler::AddSubsidiaryHandler
+ * 
+ */
+class AliMCAuxHandler : public AliMCEventHandler
+{
+public:
+  /** 
+   * Constructor 
+   * 
+   * @param name Name 
+   */
+  AliMCAuxHandler(const char* name="FMD",
+                    const char* clsName="AliFMDHit",
+                    AliMCEventHandler* parent=0);
+
+  virtual ~AliMCAuxHandler() {}
+  /** 
+   * @{
+   * @name Interface member functions 
+   */
+  /** 
+   * Intialize 
+   * 
+   * @param t Not used
+   * @param o Not used
+   * 
+   * @return always true 
+   */  
+  virtual Bool_t Init(TTree* t,Option_t* o) { return AliMCEventHandler::Init(t,o); }
+  /** 
+   * Initialize the input
+   * 
+   * @param opt Options 
+   * 
+   * @return true on success 
+   */
+  virtual Bool_t Init(Option_t* opt);
+  /**
+   * Called at the beginning of an event 
+   * 
+   * @param entry Entry in tree 
+   * 
+   * @return true on success
+   */
+  virtual Bool_t BeginEvent(Long64_t entry);
+  /** 
+   * Called when the input file is changed 
+   * 
+   * @param path New path 
+   *
+   * @return true on success
+   */
+  virtual Bool_t Notify() { return AliMCEventHandler::Notify(); }
+  /** 
+   * Called when the input file is changed 
+   * 
+   * @param path New path 
+   *
+   * @return true on success
+   */
+  virtual Bool_t Notify(const char* path);
+  /** 
+   * Called at the end of an event 
+   * 
+   * @return true on success
+   */
+  virtual Bool_t FinishEvent();
+  /** 
+   * Called at the end of a job 
+   * 
+   * @return true on success 
+   */
+  virtual Bool_t Terminate();
+  /** 
+   * Called at the end of a sub-job
+   * 
+   * @return true on success
+   */  
+  virtual Bool_t TerminateIO();
+  /** 
+   * Reset the I/O
+   * 
+   */
+  virtual void ResetIO();
+  /** 
+   * Load an event 
+   * 
+   * @param iev Event number 
+   * 
+   * @return true on success 
+   */
+  virtual Bool_t LoadEvent(Int_t iev);
+  /** 
+   * Set the number of events in the container 
+   * 
+   * @param nev Number of events 
+   */
+  virtual void  SetNumberOfEventsInContainer(Int_t nev) {
+    this->fNEventsInContainer = nev;}
+  virtual Bool_t OpenFile(Int_t ev);
+  /* @} */
+
+  AliMCEventHandler* GetParent() { return fParent; }
+  /** 
+   * Get the tree 
+   * 
+   * @return The connected hits tree
+   */
+  virtual TTree*  GetTree() const { return fTree;}  
+  TClonesArray*  GetArray() const { return fArray; }
+  Int_t         GetNEntry() const;
+  TClonesArray* GetEntryArray(Int_t entry);
+  /** 
+   * Static member function to create and attach this handler 
+   * 
+   * @param name Name of the handler 
+   * 
+   * @return Newly allocated handler or null
+   */
+  static AliMCAuxHandler* Create(const char* name="FMD",
+                                   const char* what="Hits");
+  static TClonesArray* GetParticleArray(AliMCAuxHandler* handler, 
+                                       Int_t particle);
+protected:
+  AliMCAuxHandler(const AliMCAuxHandler& o)
+    : AliMCEventHandler(), 
+      fParent(o.fParent),
+      fFile(0), 
+      fTree(0), 
+      fDir(0),
+      fArray(0),
+      fNEvents(0), 
+      fNEventsPerFile(0),
+      fNEventsInContainer(0),
+      fEvent(0), 
+      fFileNumber(0), 
+      fTreeName(""),
+      fFileBase("")
+  {}
+  AliMCAuxHandler& operator=(const AliMCAuxHandler& o)
+  {
+    if (&o == this) return *this;
+    // AliMCEventHandler::operator=(o);
+    fParent = o.fParent;
+    fFile   = o.fFile;
+    fTree   = o.fTree;
+    return *this;
+  }
+  TString* GetParentPath() const;
+  AliMCEventHandler* fParent; // Parent MC handler 
+  TFile*             fFile;                //!
+  TTree*             fTree;                //!
+  TDirectory*        fDir;                 //!
+  TClonesArray*      fArray;               //!
+  Int_t              fNEvents;             //!
+  Int_t              fNEventsPerFile;      //!
+  Int_t              fNEventsInContainer;  //!
+  Int_t              fEvent;               //!
+  Int_t              fFileNumber;          //!
+  TString            fTreeName;            //!
+  TString            fFileBase;            //! 
+  ClassDef(AliMCAuxHandler,1); // Connect FMD hits tree
+};
+
+#endif
+// Local Variables:
+//  mode: C++
+// End:
+
index df275ba6d1fffb7d1184d78f1af3453dfa173b55..6fc2bc0ed2d16dacd8ebe507700deaa11228832d 100644 (file)
@@ -55,7 +55,7 @@ struct QAPlotter : public QABase
      * @param d Detector 
      * @param r Ring 
      */
-    Ring(UShort_t d, Char_t r)
+    Ring(UShort_t d, Char_t r, Bool_t useVar=false)
       : QARing(d, r),
         fGChi2(0),
         fGC(0),
@@ -66,7 +66,8 @@ struct QAPlotter : public QABase
         fGSingles(0),
         fGLoss(0),
         fGBeta(0),
-        fGOccupancy(0)
+        fGOccupancy(0),
+       fUseVar(useVar)
     {
       fGChi2           = new TGraphAsymmErrors;
       fGC              = new TGraphAsymmErrors;
@@ -143,6 +144,12 @@ struct QAPlotter : public QABase
       Double_t y  = q.mean;
       Double_t el = y-q.min;
       Double_t eh = q.max-y;
+      if (fUseVar) { 
+       Info("UpdateGraph", "Setting errors on %s to variance",
+            g->GetName());
+       el = q.var; 
+       eh = q.var; 
+      }
       if (TMath::Abs(y) < 1e-6) return;
       Int_t    i  = g->GetN();
       g->SetPoint(i, runNo, y);
@@ -185,20 +192,23 @@ struct QAPlotter : public QABase
     TGraph*            fGLoss;     // Graph of 'lost' data 
     TGraph*            fGBeta;     // Graph of Poisson vs ELoss correlation
     TGraphAsymmErrors* fGOccupancy;// Graph of mean occupancy              
+    Bool_t             fUseVar;    // Use variance 
   };
   /** 
    * Constructor 
    */
-  QAPlotter(Int_t prodYear=0, Char_t prodLetter='\0') 
+  QAPlotter(Int_t prodYear=0, Char_t prodLetter='\0', Bool_t useVar=false
     : QABase(false, prodYear, prodLetter),
       fNAccepted(0),
-      fVz(0)
+      fVz(0),
+      fUseVar(useVar)
   {
-    fFMD1i = new Ring(1, 'I'); 
-    fFMD2i = new Ring(2, 'I'); 
-    fFMD2o = new Ring(2, 'O'); 
-    fFMD3i = new Ring(3, 'I'); 
-    fFMD3o = new Ring(3, 'O'); 
+    Info("QAPlotter", "Do we use variance? %s", fUseVar ? "yes" : "no");
+    fFMD1i = new Ring(1, 'I', useVar); 
+    fFMD2i = new Ring(2, 'I', useVar); 
+    fFMD2o = new Ring(2, 'O', useVar); 
+    fFMD3i = new Ring(3, 'I', useVar); 
+    fFMD3o = new Ring(3, 'O', useVar); 
     fNAccepted = new TGraph;
     fNAccepted->SetName("nAccepted");
     fNAccepted->SetMarkerStyle(20);
@@ -254,6 +264,7 @@ struct QAPlotter : public QABase
     UInt_t nEntries = fTree->GetEntries();
     UInt_t j = 0;
     fRuns.Set(nEntries);
+    Info("Run", "Got %d runs", nEntries);
     for (UInt_t i = 0; i < nEntries; i++) {
       fTree->GetEntry(i);
 
@@ -503,13 +514,18 @@ struct QAPlotter : public QABase
   {
     h->GetXaxis()->SetNoExponent();
     // h->GetXaxis()->SetTitleOffset(1);
-    h->SetYTitle(title);
+    TString ytitle(title);
+    if (fUseVar) ytitle.Append(" (errors: variance)");
+    else         ytitle.Append(" (errors: min/max)");
+    h->SetYTitle(ytitle.Data());
     h->SetXTitle("Run #");
+    Info("AddRuns", "%s: %s vs %s", h->GetName(), 
+        h->GetXaxis()->GetTitle(), h->GetYaxis()->GetTitle());
 
     Int_t    r1  = h->GetXaxis()->GetXmin();
     Int_t    r2  = h->GetXaxis()->GetXmax();
     Double_t lx  = 0;
-    Double_t tx  = .045; // (r2 - r1) / 18;
+    Double_t tx  = .025; // (r2 - r1) / 18;
     Double_t wx  = 1 - fCanvas->GetLeftMargin() - fCanvas->GetRightMargin();
     Double_t dy  = .025;
     Double_t y   = fCanvas->GetBottomMargin()+dy;
@@ -620,6 +636,7 @@ struct QAPlotter : public QABase
   UInt_t        fLast;      // Last run
   TArrayI       fRuns;      // Seen runs 
   TObjArray     fFiles;
+  Bool_t        fUseVar;    // Use variance rather than min/max 
 };
  
 //
index 7582eabaf49e9e5481816c0d4133975707b05dfe..7f617df9f3c39931a9ac9b484362b34b06370eed 100644 (file)
@@ -139,7 +139,8 @@ GetListOfFiles(const char* input=".")
  * @ingroup pwglf_forward_qa_scripts
  */
 void 
-RunFinalQA(const char* dir, Int_t prodYear=0, const char* prodLetter="")
+RunFinalQA(const char* dir, Int_t prodYear=0, const char* prodLetter="",
+          Bool_t useVar=false)
 {
    int ret = 0;
    gROOT->SetMacroPath(Form(".:%s",gROOT->GetMacroPath()));
@@ -149,7 +150,9 @@ RunFinalQA(const char* dir, Int_t prodYear=0, const char* prodLetter="")
    gROOT->LoadMacro("QABase.h+g");
    gROOT->LoadMacro("QAPlotter.C+g");
 
-   QAPlotter p(prodYear, prodLetter[0]);
+   Info("RunFinalQA", "Final QA: %d%c (variance: %s)", 
+       prodYear, prodLetter[0], (useVar ? "true" : "false"));
+   QAPlotter p(prodYear, prodLetter[0], useVar);
   
    TList* l = GetListOfFiles(dir);
    TIter next(l);
diff --git a/PWGLF/FORWARD/analysis2/qa/RunGridQA.sh b/PWGLF/FORWARD/analysis2/qa/RunGridQA.sh
new file mode 100755 (executable)
index 0000000..e354143
--- /dev/null
@@ -0,0 +1,238 @@
+#!/bin/bash
+
+# --- Usage ----------------------------------------------------------
+usage()
+{
+    cat <<EOF
+Usage: $0 [OPTIONS]
+
+Options:
+       -p,--production PERIOD      LHC Period identifier [$prod]
+       -P,--pass       PASS        Production pass identifier [$prodpass]
+       -n,--part       PART        Part identifier [$part]
+       -r,--runs       RUNS        List of runs or file [$runs]
+       -c,--clean                  Clean previous job output [$clean]
+       -v,--verbose                Increase verbosity [$verb]
+EOF
+
+}
+
+# --- Check AliEn token ----------------------------------------------
+check_token()
+{
+    uid=`id -u`
+    genv_file=/tmp/gclient_env_${uid}
+    
+    if test ! -f ${genv_file} ; then 
+       echo "No such file: ${genv_file}, please do alien-token-init" \
+           >/dev/stderr
+       exit 1
+    fi
+    . ${genv_file}
+    alien-token-info | grep -q "Token is still valid"
+    if test $? -ne 0 ; then 
+       echo "Token not valid, please re-new" > /dev/stderr 
+       exit 1
+    fi
+}
+
+# --- Diagnostics output ---------------------------------------------
+verb=0
+mess()
+{
+    if test $1 -gt $verb ; then return ; fi 
+    shift
+    echo $*
+}
+
+# --- Parse production information -----------------------------------
+parse_prod()
+{
+    prodyear=`echo $prodfull | sed 's/LHC\(..\).*/\1/'` 
+    prodletter=`echo $prodfull | sed "s/LHC${prodyear}\(.\).*/\1/"` 
+    prodpost=`echo $prodfull | sed "s/LHC${prodyear}${prodletter}//"` 
+}
+
+parse_pass()
+{
+    passno=`echo $passfull | sed 's/.*pass\([0-9]*\).*/\1/'`  
+    passpost=`echo $passfull | sed -n "s/.*pass${passno}_//p"` 
+    passpre=`echo $passfull | sed -n "s/pass.*//p"` 
+}
+
+# --- Append path element --------------------------------------------
+append_to_path()
+{
+    local tmp=$1 ; shift 
+    local add=$1
+    if test "x$tmp" != "x" ; then tmp="${tmp}/" ; fi 
+    echo ${tmp}${add}
+}
+
+# --- Get the path and search pattern --------------------------------
+path=
+search=
+setup_input()
+{
+    local datd=data
+    local esdd=ESDs/
+    case x$prodpost in 
+       x_*) ;; 
+       x) ;; 
+       *)  mess 3 "Assuming simulation output"
+           datd=sim 
+           esdd= 
+           ;; 
+    esac
+    
+    local paid=
+    if echo "$passno" | grep -q -E '^[0-9]*[.]?[0-9]*$' ; then 
+       if test "x$passfull" != "x" && test $passno -gt 0 ; then 
+           paid=pass${passno}
+       fi
+    else
+       paid=${passfull}
+       passpre=
+       post=
+    fi
+    local post=${passpost}
+    case x$post in 
+       x_*) ;; 
+       x) ;; 
+       *) post="_${post}" ;; 
+    esac
+
+    # Assume official productions 
+    path=/alice/${datd}/${year}/${prodfull}/
+    search="${esdd}${passpre}${paid}${post}"
+    search=`append_to_path "$search" "*/AliESDs.root"` 
+}
+
+# --- Setup the runs -------------------------------------------------
+unique=
+setup_runs()
+{
+    if test -f $runs ; then 
+       unique=`basename $runs .txt` 
+       unique=`basename $unique .list` 
+       unique=`basename $unique .runs` 
+       return
+    fi
+    local l=`echo "$runs" | tr ',+:.' ' '` 
+    local first=
+    local last=
+    for i in $l ; do 
+       if test "x$first" = "x" ; then first=$i ; fi 
+       last=$l 
+    done
+    unique=${first}_${last} 
+    runs=`echo "$l" | tr ' ' ','` 
+}
+
+# --- Clean previous attempt -----------------------------------------
+clean_previous()
+{
+    local alienu=`alien_whoami | tr -d ' '` 
+    local alienl=`echo ${alienu} | cut -c1` 
+    local aliend="/alice/cern.ch/user/${alienl}/${alienu}/${name}"
+    mess 1 "Removing local directory ${name}"
+    if test $noact -lt 1 ; then rm -rf ${name} ; fi
+    mess 1 "Removing alien directory ${aliend}"
+    if test $noact -lt 1 ; then alien_rmdir ${aliend} ; fi
+}
+
+# --- Parse command line options -------------------------------------
+redir=/dev/null
+part=
+clean=0
+noact=0
+while test $# -gt 0 ; do 
+    case $1 in 
+       -h|--help)        usage ; exit 0 ;; 
+       -v|--verbose)     let verb=$verb+1   ;; 
+       -l|--log-file)    redir=             ;; 
+       -p|--production)  prodfull=$2; shift; parse_prod ; year=20${prodyear} ;;
+       -P|--pass)        passfull=$2; shift; parse_pass ;;
+       -n|--part)        part=$2 ; shift ;; 
+       -r|--runs)        runs=$2               ; shift ;;
+       -c|--clean)       clean=1 ; shift ;; 
+       -s|--noact)       noact=1 ; shift ;;
+       *) echo "$0: Unknown argument: $1" > /dev/stderr ; exit 1 ;; 
+    esac
+    shift
+done
+
+# --- Initial checks -------------------------------------------------
+if test "x$prodfull" = "x" ; then 
+    echo "No production specified" > /dev/stderr 
+    exit 1
+fi
+if test "x$passfull" = "x" ; then 
+    echo "No pass specified" > /dev/stderr 
+    exit 1
+fi
+if test "x$runs" = "x" ; then 
+    echo "No runs specified" > /dev/stderr 
+    exit 1
+fi
+
+check_token 
+setup_input
+setup_runs
+name="QA_${prodfull}_${passfull}"
+if test "x$part" != "x" ; then 
+    if echo "$part" | grep -q -E '^[0-9]*[.]?[0-9]*$' ; then 
+       part="part${part}"
+    fi
+    name="${name}_${part}"
+fi
+
+if test $clean -gt 0 ; then clean_previous ; fi 
+    
+
+# --- Some friendly information --------------------------------------
+cat <<EOF
+       Year:                   ${year}
+       Production:             ${prodfull} 
+         Year:                 ${prodyear}
+         Letter:               ${prodletter}
+         Suffix:               ${prodpost}
+       Pass:                   ${passfull}
+         Number:               ${passno}
+         Prefix:               ${passpre}
+         Postfix:              ${passpost}
+       Lock file:              ${lock}
+       Log:                    ${redir}
+       Runs:                   ${runs}
+       Path:                   ${path}
+       Search:                 ${search}
+       Part:                   ${part}
+       Name:                   ${name}
+EOF
+
+# --- Do the actual running ------------------------------------------
+url="alien:///alice/data/${year}/${prodfull}\?run=${runs}\&pattern=${search}#esdTree"
+train=MakeQATrain
+opts=(--class=${train} \
+    --name=${name} \
+    --cent \
+    --url="${url}")
+mess 1 "Will do runTrain \"${opts[@]}\" $@" 
+if test $noact -gt  0 ; then 
+    sleep 1 
+else 
+    runTrain "${opts[@]}" $@
+fi
+
+mess 1 "Will monitor in directory $name" 
+if test $noact -gt 0 ; then 
+    sleep 1 
+else 
+    (cd ${name} && aliroot -b -x -q -l Watch.C)
+fi
+
+#
+# EOF
+#
+
+
index b97e02050d7eda8866498448831ec3c5d7aca90a..12567bc9b025616f741854656bd7d1132a170b2d 100755 (executable)
@@ -31,13 +31,18 @@ Options:
        -l,--log-file              Log file output [$redir]
        -b,--barrel     MODE       Fetch barrel data  [$barrel]
        -f,--force                 Force re-run analysis [$force]
+       -V,--variance              Errors=variance (not min/max) [$variance]
+        -L,--local                 Local trending_<X>.root files [$from_local]
+       -d,--directory  DIR        Search custom AliEn directory [$path]
 
 Note the option -j and the options -p and -P are mutually exclusive,
 The option -Q is only used if the options -p and -P are given.
 Production identifiers are of the form LHC11h, LHC11h3, or LHC11h_2. 
-Pass identifers are of the form pass2, pass1_HLT, or cpass1.  If barrel mode
-is >0, then do not assume ESD directory.  If mode>1, then get 
-trending_barrel.root and QAresults_barrel.root
+Pass identifers are of the form pass2, pass1_HLT, or cpass1.  
+If barrel mode>0, then do not assume ESD directory.  
+If barrel mode>1, then get trending_barrel.root and QAresults_barrel.root
+Option -d is for hand-made QA passes. 
+If optiond -d is not specified then official QA passes are assumed.
 EOF
 }
 
@@ -214,10 +219,11 @@ other=QAresults.root
 files=
 path=
 numf=0
+from_local=0
 get_filelist()
 {
     mess 3 "Getting file list" 
-    
+
     local datd=data
     local esdd=ESDs/
     if test ${barrel} -gt 0 ; then
@@ -253,8 +259,14 @@ get_filelist()
        *) post="_${post}" ;; 
     esac
 
-    path=/alice/${datd}/${year}/${prodfull}/
-    local search="${esdd}${passpre}${paid}${post}"
+    local search=
+    if test "x$path" = "x" ; then 
+       # Assume official productions 
+       path=/alice/${datd}/${year}/${prodfull}/
+       search="${esdd}${passpre}${paid}${post}"
+    else
+       search="*"
+    fi
 
     if test $qanumber -gt 0 ; then 
        qapost=`printf "QA%02d" $qanumber` 
@@ -267,9 +279,14 @@ get_filelist()
        Path:                   $path
        Search:                 $search
 EOF
-    mess 1 "Getting list of files from AliEn - can take minutes - be patient"
-    mess 2 "alien_find ${path} ${search}"
-    files=`alien_find ${path} ${search} | grep -v "\(files found\|AND THE\)" 2>> ${redir}` 
+    if test $from_local -lt 1 ; then 
+
+       mess 1 "Getting list of files from AliEn - can take minutes - be patient"
+       mess 2 "alien_find ${path} ${search}"
+       files=`alien_find ${path} ${search} | grep -v "\(files found\|AND THE\)" 2>> ${redir}` 
+    else 
+       files=`ls ${top}/${store}/trending_*.root | sed 's,${top}/${store}/,,g'`
+    fi
     for i in $files ; do 
        let numf=$numf+1
     done 
@@ -285,8 +302,11 @@ EOF
 fix_perm()
 {
     # if test ! -f $1 ; then return ; fi 
-    chmod g+rwX $1 >> ${redir} 2>&1 
-    chmod o+rX $1 >> ${redir} 2>&1 
+    chmod g+rwX $1 >> /dev/null 2>&1 
+    chmod o+rX $1 >> /dev/null 2>&1 
+    # 
+    # chmod g+rwX $1 >> ${redir} 2>&1 
+    # chmod o+rX $1 >> ${redir} 2>&1 
 }
 
 # --- Check if a file is OK ------------------------------------------
@@ -399,8 +419,13 @@ submit_runs()
     for i in $@ ; do 
        let cur=$sta+$counter
 
-       local b=`echo $i | sed -e "s,${path},,"` 
-       local r=`echo $b | sed -e "s,/.*,,"` 
+       local r
+       if test $from_local -lt 1 ; then 
+           local b=`echo $i | sed -e "s,${path},,"` 
+           r=`echo $b | sed -e "s,/.*,,"` 
+       else
+           r=`basename $i .root | sed 's/trending_//'` 
+       fi
 
        printf "%3d/%3d: %s\n" $cur $max $r 
        runs[$counter]=$r
@@ -476,11 +501,12 @@ copy_style()
 }      
 
 # --- Run the final trending -----------------------------------------
+variance=0
 make_trend()
 {
     local dir=$1 
     local ret=0
-    mess 1 -n "Analysing $dir ... "
+    mess 1 -n "Analysing for trend $dir ... "
     (cd $dir 
        rm -f trend_*_*.html 
        rm -f trend_*_*.pdf
@@ -488,7 +514,7 @@ make_trend()
 
        root -l -b <<EOF 
 .L RunFinalQA.C
-RunFinalQA(".", $prodyear, "$prodletter");
+RunFinalQA(".", $prodyear, "$prodletter", $variance);
 .q
 EOF
        mess 1 -n " ... "
@@ -589,12 +615,15 @@ while test $# -gt 0 ; do
        -R|--also-results) also_results=1     ;; 
        -Q|--qa-number)    qanumber=$2        ; shift ;;
        -l|--log-file)     redir=             ;; 
+       -L|--local)        from_local=0       ;;
+       -V|--variance)     variance=1         ;;
        -b|--barrel)       barrel=$2          ; shift ;;
        -f|--force)        let force=$force+1 ;; 
        -p|--production) 
            prodfull=$2; shift; parse_prod ; year=20${prodyear} ;;
        -P|--pass) 
            passfull=$2; shift; parse_pass ;;
+       -d|--directory)    path=$2              ; shift ;;
        *) echo "$0: Unknown argument: $1" > /dev/stderr ; exit 1 ;; 
     esac
     shift
@@ -674,6 +703,8 @@ cat <<EOF
        Lock file:              ${lock}
        Log:                    ${redir}
        Force:                  ${force}
+       Use variance:           ${variance}
+       Use pre-downloaded:     ${from_local}
 EOF
 # --- Do a search to find our files ----------------------------------
 get_filelist