Minor fixes
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Mar 2006 23:51:58 +0000 (23:51 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Mar 2006 23:51:58 +0000 (23:51 +0000)
14 files changed:
FMD/AliFMDIndex.cxx [new file with mode: 0644]
FMD/AliFMDIndex.h [new file with mode: 0644]
FMD/AliFMDRawReader.cxx
FMD/AliFMDRawStream.cxx
FMD/Doxyfile
FMD/FMDbaseLinkDef.h
FMD/Reconstruct.C
FMD/libFMDbase.pkg
FMD/scripts/GeoGeometry.C [new file with mode: 0644]
FMD/scripts/GetMedia.C
FMD/scripts/NodeGeometry.C [new file with mode: 0644]
FMD/scripts/SimpleGeometry.C [deleted file]
FMD/scripts/TestIndex.C [new file with mode: 0644]
FMD/scripts/WriteMedArrays.C [new file with mode: 0644]

diff --git a/FMD/AliFMDIndex.cxx b/FMD/AliFMDIndex.cxx
new file mode 100644 (file)
index 0000000..799a3d3
--- /dev/null
@@ -0,0 +1,158 @@
+/**************************************************************************
+ * Copyright(c) 2004, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+/* $Id$ */
+/** @file    AliFMDDigit.cxx
+    @author  Christian Holm Christensen <cholm@nbi.dk>
+    @date    Mon Mar 27 12:37:41 2006
+    @brief   Digits for the FMD 
+*/
+//////////////////////////////////////////////////////////////////////
+//
+//  Class that holds an FMD index.  That is, it holds the detector
+//  coordinates for a given strip:
+//
+//     Variable | Type     | Range   | Description
+//     ---------+----------+---------+------------------
+//     detector | UShort_t | 1-3     | Detector number 
+//     ring     | Char_t   | 'I'/'O' | Ring identifier 
+//     sector   | UShort_t | 0-39    | Sector number
+//     strip    | UShort_t | 0-511   | Strip number
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "AliFMDIndex.h"       // ALIFMDINDEX_H
+#include "Riostream.h"         // ROOT_Riostream
+#include <TString.h>            // ROOT_TString
+#include <AliFMDMap.h>
+
+//====================================================================
+ClassImp(AliFMDIndex)
+#if 0
+  ; // This is here to keep Emacs from indenting the next line
+#endif
+
+//____________________________________________________________________
+AliFMDIndex::AliFMDIndex()
+  : fDetector(0), 
+    fRing('\0'), 
+    fSector(0), 
+    fStrip(0), 
+    fHash(-1)
+{}
+
+//____________________________________________________________________
+AliFMDIndex::AliFMDIndex(const AliFMDIndex& o)
+  : fDetector(o.fDetector), 
+    fRing(o.fRing), 
+    fSector(o.fSector), 
+    fStrip(o.fStrip), 
+    fHash(o.fHash)
+{
+  // Copy constructor 
+}
+
+//____________________________________________________________________
+AliFMDIndex::AliFMDIndex(UShort_t detector, 
+                        Char_t   ring, 
+                        UShort_t sector, 
+                        UShort_t strip)
+  : fDetector(detector), 
+    fRing(ring), 
+    fSector(sector), 
+    fStrip(strip), 
+    fHash(-1)
+{
+  //
+  // Creates a base data digit object
+  //
+  // Parameters 
+  //
+  //    detector  Detector # (1, 2, or 3)                      
+  //    ring     Ring ID ('I' or 'O')
+  //    sector   Sector # (For inner/outer rings: 0-19/0-39)
+  //    strip    Strip # (For inner/outer rings: 0-511/0-255)
+}
+
+//____________________________________________________________________
+AliFMDIndex& 
+AliFMDIndex::operator=(const AliFMDIndex& o)
+{
+  // Assignment operator 
+  fDetector = o.fDetector;
+  fRing     = o.fRing;
+  fSector   = o.fSector;
+  fStrip    = o.fStrip;
+  fHash     = o.fHash;
+  return *this;
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDIndex::Hash() const 
+{
+  if (fHash < 0) {
+    size_t ringi = (fRing == 'I' ||  fRing == 'i' ? 0 : 1);
+    fHash = (fStrip + 
+            AliFMDMap::kMaxStrips * 
+            (fSector + AliFMDMap::kMaxSectors * 
+             (ringi + AliFMDMap::kMaxRings * (fDetector-1))));
+  }
+  return fHash;
+}
+
+
+//____________________________________________________________________
+void
+AliFMDIndex::Print(Option_t* /* option*/) const 
+{
+  // Print digit to standard out 
+  cout << Name() << flush;
+}
+
+//____________________________________________________________________
+const char*
+AliFMDIndex::Name() const 
+{ 
+  if (fName.IsNull()) 
+    fName = Form("FMD%d%c[%2d,%3d]", fDetector, fRing, fSector, fStrip);
+  return fName.Data();
+}
+
+//====================================================================
+ClassImp(AliFMDObjIndex)
+#if 0
+  ; // This is here to keep Emacs from indenting the next line
+#endif
+
+//____________________________________________________________________
+Int_t 
+AliFMDObjIndex::Compare(const TObject* o) const
+{
+  const AliFMDObjIndex* a = dynamic_cast<const AliFMDObjIndex*>(o);
+  if (!a) {
+    Fatal("Compare", 
+         "trying to compare to something not a AliFMDObjIndex object, "
+         "but a %s object", o->ClassName());
+    return 0;
+  }
+  if (this->operator<(*a)) return -1;
+  if (this->operator==(*a)) return 0;
+  return 1;
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/AliFMDIndex.h b/FMD/AliFMDIndex.h
new file mode 100644 (file)
index 0000000..a16f0bf
--- /dev/null
@@ -0,0 +1,173 @@
+#ifndef ALIFMDINDEX_H
+#define ALIFMDINDEX_H
+/** @file    AliFMDIndex.h
+    @author  Christian Holm Christensen <cholm@nbi.dk>
+    @date    Mon Mar 27 12:37:41 2006
+    @brief   FMD detector coordinates
+*/
+//___________________________________________________________________
+//
+//  Class that holds an FMD index.  That is, it holds the detector
+//  coordinates for a given strip:
+//
+//     Variable | Type     | Range   | Description
+//     ---------+----------+---------+------------------
+//     detector | UShort_t | 1-3     | Detector number 
+//     ring     | Char_t   | 'I'/'O' | Ring identifier 
+//     sector   | UShort_t | 0-39    | Sector number
+//     strip    | UShort_t | 0-511   | Strip number
+//
+#ifndef ROOT_Rtypes
+# include <Rtypes.h>
+#endif
+#ifndef ROOT_TObject
+# include <TObject.h>
+#endif
+#ifndef ROOT_TString
+# include <TString.h>
+#endif
+#include <iosfwd>
+
+//____________________________________________________________________
+/** @class AliFMDIndex AliFMDIndex.h <FMD/AliFMDIndex.h>
+    @brief FMD detector coordinates 
+    @ingroup FMD_base
+ */
+class AliFMDIndex 
+{
+public: 
+  /** CTOR */
+  AliFMDIndex();
+  /** Copy CTOR 
+      @param o Object to copy from */
+  AliFMDIndex(const AliFMDIndex& o);
+  /** Constrctor 
+      @param detector Detector 
+      @param ring     Ring
+      @param sector   Sector
+      @param strip    Strip */
+  AliFMDIndex(UShort_t detector, 
+             Char_t   ring='\0', 
+             UShort_t sector=0, 
+             UShort_t strip=0);
+  /** Assignment operator 
+      @param o Object to assign from 
+      @return Reference to this object  */
+  AliFMDIndex& operator=(const AliFMDIndex& o);
+  /** Comparison operator 
+      @param o Object to compare to 
+      @return @c true if these refer to the same index  */
+  bool operator==(const AliFMDIndex& o) const;
+  /** Comparison operator 
+      @param o Object to compare to 
+      @return @c true if this is smaller than @a o */
+  bool operator<(const AliFMDIndex& o) const;
+  /** DTOR */
+  virtual ~AliFMDIndex() {}
+  /** @return Detector # */
+  UShort_t     Detector()         const { return fDetector; }
+  /** @return Ring ID */
+  Char_t       Ring()             const { return fRing;     }
+  /** @return sector # */
+  UShort_t     Sector()                   const { return fSector;   }
+  /** @return strip # */
+  UShort_t     Strip()            const { return fStrip;    }
+  /** @param x Detector # */
+  void SetDetector(UShort_t x)    { fHash = -1; fDetector = x; }
+  /** @param x Ring ID */
+  void SetRing(Char_t x)          { fHash = -1; fRing = x; }
+  /** @param x sector # */
+  void SetSector(UShort_t x)      { fHash = -1; fSector = x; }
+  /** @param x strip # */
+  void SetStrip(UShort_t x)       { fHash = -1; fStrip = x; }
+  /** Print information 
+      @param opt Not used */
+  virtual void Print(Option_t* opt="") const;
+  /** @return Name */
+  const char*  Name() const;
+protected:
+  Int_t Hash() const;
+  UShort_t fDetector;      // (Sub) Detector # (1,2, or 3)
+  Char_t   fRing;          // Ring ID ('I' or 'O')
+  UShort_t fSector;        // Sector # (phi division)
+  UShort_t fStrip;         // Strip # (radial division)
+  mutable TString  fName;  //! Cached name
+  mutable Int_t    fHash;  //! Cached hash value
+  ClassDef(AliFMDIndex, 1) // Base class for FMD digits 
+};
+
+//____________________________________________________________________
+class AliFMDObjIndex : public TObject, public AliFMDIndex
+{
+public:
+  /** CTOR */
+  AliFMDObjIndex() {}
+  /** Copy CTOR 
+      @param o Object to copy from */
+  AliFMDObjIndex(const AliFMDObjIndex& o) : TObject(o), AliFMDIndex(o) {}
+  /** Construct from a pure index
+      @param o Object to copy from */
+  explicit AliFMDObjIndex(const AliFMDIndex& o) : AliFMDIndex(o) {}
+  /** Constrctor 
+      @param detector Detector 
+      @param ring     Ring
+      @param sector   Sector
+      @param strip    Strip */
+  AliFMDObjIndex(UShort_t detector, 
+                Char_t   ring='\0', 
+                UShort_t sector=0, 
+                UShort_t strip=0) 
+    : AliFMDIndex(detector, ring, sector, strip)
+  {}
+  /** DTOR */
+  virtual ~AliFMDObjIndex() {}
+  /** @return name */
+  virtual const char* GetName() const { return AliFMDIndex::Name(); }
+  /** sort compare for TCollection's
+      @param o Object to compare to
+      @return  -1 if this is @e smaller than @a o, 0 if @e equal to 
+      @a o, and 1 if this is @e larger than @a o */
+  virtual Int_t Compare(const TObject* o) const;
+  /** @return always true */
+  Bool_t IsSortable() const { return kTRUE; }
+  ClassDef(AliFMDObjIndex, 1) // Base class for FMD digits 
+};
+//____________________________________________________________________
+inline
+bool
+AliFMDIndex::operator==(const AliFMDIndex& o) const
+{
+  return (o.Hash() == Hash());
+}
+
+//____________________________________________________________________
+inline
+bool
+AliFMDIndex::operator<(const AliFMDIndex& rhs) const
+{
+  return (Hash() < rhs.Hash());
+}
+
+#if 0
+//____________________________________________________________________
+inline
+bool
+operator<(const AliFMDIndex& lhs, const AliFMDIndex& rhs)
+{
+  return (lhs.Detector() < rhs.Detector() ? true : 
+         (lhs.Ring()    < rhs.Ring()     ? true :
+          (lhs.Sector() < rhs.Sector()   ? true :
+           (lhs.Strip() < rhs.Strip()    ? true : false))));
+}
+#endif 
+#endif
+//____________________________________________________________________
+//
+// Local Variables:
+//   mode: C++
+// End:
+//
+//
+// EOF
+//
index a730238..5e8c7ff 100644 (file)
@@ -112,8 +112,6 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
   // Get sample rate 
   AliFMDParameters* pars = AliFMDParameters::Instance();
   AliFMDRawStream input(fReader);
-  // Select FMD DDL's 
-  fReader->Select(AliFMDParameters::kBaseDDL>>8);
 
   UShort_t stripMin = 0;
   UShort_t stripMax = 127;
index 6cc5a33..7c47e95 100644 (file)
@@ -43,6 +43,8 @@ AliFMDRawStream::AliFMDRawStream(AliRawReader* reader)
   : AliAltroRawStream(reader)
 {
   fNoAltroMapping = kFALSE;
+  // Select FMD DDL's 
+  SelectRawData(AliFMDParameters::kBaseDDL>>8);
 }
 
 //_____________________________________________________________________________
@@ -50,7 +52,6 @@ Bool_t
 AliFMDRawStream::ReadChannel(UInt_t& ddl, UInt_t& addr, 
                             UInt_t& len, UShort_t* data)
 {
-  UInt_t       prevddl   = 0;
   Int_t        l         = 0;
   static Int_t last      = 0xFFFF; // 0xFFFF means signal is used
   Bool_t       next      = kTRUE;
@@ -65,13 +66,12 @@ AliFMDRawStream::ReadChannel(UInt_t& ddl, UInt_t& addr,
        AliDebug(15, Form("New hardware address, was 0x%x, now 0x%x", 
                          GetPrevHWAddress(), GetHWAddress()));
        addr = GetPrevHWAddress();
-       ddl  = AliFMDParameters::kBaseDDL + prevddl;
-       len  = l+1;
+       ddl  = AliFMDParameters::kBaseDDL + GetPrevDDLNumber();
+       len  = l+1; // Need to add one - l points to last valid index
        last = signal;
        break;
       }
     }
-    prevddl  = GetPrevDDLNumber();
     Int_t t  = GetTime();
     l        = TMath::Max(l, t);
     data[t]  = signal;
index 5028e20..0ef205f 100644 (file)
@@ -207,18 +207,18 @@ PERL_PATH              = /usr/bin/perl
 # Configuration options related to the dot tool   
 #---------------------------------------------------------------------------
 CLASS_DIAGRAMS         = YES
-HIDE_UNDOC_RELATIONS   = YES
+HIDE_UNDOC_RELATIONS   = NO
 HAVE_DOT               = YES
 CLASS_GRAPH            = YES
 COLLABORATION_GRAPH    = YES
-GROUP_GRAPHS           = YES
+GROUP_GRAPHS           = NO
 UML_LOOK               = YES
 TEMPLATE_RELATIONS     = NO
-INCLUDE_GRAPH          = YES
-INCLUDED_BY_GRAPH      = YES
+INCLUDE_GRAPH          = NO
+INCLUDED_BY_GRAPH      = NO
 CALL_GRAPH             = NO
 GRAPHICAL_HIERARCHY    = YES
-DIRECTORY_GRAPH        = YES
+DIRECTORY_GRAPH        = NO
 DOT_IMAGE_FORMAT       = png
 DOT_PATH               = 
 DOTFILE_DIRS           = 
index 75fdfa6..daff1c5 100644 (file)
@@ -13,6 +13,8 @@
 #pragma link off all classes;
 #pragma link off all functions;
  
+#pragma link C++ class  AliFMDIndex+;
+#pragma link C++ class  AliFMDObjIndex+;
 #pragma link C++ class  AliFMDBaseDigit+;
 #pragma link C++ class  AliFMDDigit+;
 #pragma link C++ class  AliFMDSDigit+;
index 51570bb..e89c638 100644 (file)
@@ -29,6 +29,20 @@ Reconstruct()
   cdb->SetDefaultStorage("local://$ALICE_ROOT");
   AliLog::SetModuleDebugLevel("FMD", 2);
   AliReconstruction rec;   
+  AliCDBEntry* align = cdb->Get("FMD/Align/Data");
+  if (align) {
+    TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
+    if (array) {
+      Int_t nAlign = array->GetEntries();
+      for (Int_t i = 0; i < nAlign; i++) {
+       AliAlignObjAngles* a = static_cast<AliAlignObjAngles*>(array->At(i));
+       if (!a->ApplyToGeometry()) {
+         Warning("ApplyAlignement", "Failed to apply alignment to %s", 
+                 a->GetVolPath());
+       }
+      }
+    }
+  }
   rec.SetRunLocalReconstruction("FMD");
   rec.SetRunVertexFinder(kFALSE);
   rec.SetRunTracking(""); 
index 1bc646c..80559af 100644 (file)
@@ -2,7 +2,8 @@
 #
 # $Id$
 
-SRCS           =  AliFMDDigit.cxx              \
+SRCS           =  AliFMDIndex.cxx              \
+                  AliFMDDigit.cxx              \
                   AliFMDBoolMap.cxx            \
                   AliFMDUShortMap.cxx          \
                   AliFMDCalibPedestal.cxx      \
diff --git a/FMD/scripts/GeoGeometry.C b/FMD/scripts/GeoGeometry.C
new file mode 100644 (file)
index 0000000..e30e4bc
--- /dev/null
@@ -0,0 +1,642 @@
+//____________________________________________________________________
+//
+//
+// $Id$
+//
+// Script I used for rapid prototyping of the FMD3 geometry - in
+// particular the support cone 
+//
+/** @defgroup geo_geom Simple geometry
+    @ingroup FMD_script
+*/
+#ifndef __CINT___
+# ifndef ROOT_TGeoVolume
+#  include <TGeoVolume.h>
+# endif
+# ifndef ROOT_TGeoMaterial
+#  include <TGeoMaterial.h>
+# endif
+# ifndef ROOT_TGeoMedium
+#  include <TGeoMedium.h>
+# endif
+# ifndef ROOT_TGeoMatrix
+#  include <TGeoMatrix.h>
+# endif
+# ifndef ROOT_TGeoXtru
+#  include <TGeoXtru.h>
+# endif
+# ifndef ROOT_TGeoPcon
+#  include <TGeoPcon.h>
+# endif
+# ifndef ROOT_TGeoTube
+#  include <TGeoTube.h>
+# endif
+# ifndef ROOT_TGeoManager
+#  include <TGeoManager.h>
+# endif
+# ifndef ROOT_TGeoPhysicalNode
+#  include <TGeoPhysicalNode.h>
+# endif
+# ifndef ROOT_TVector3
+#  include <TVector3.h>
+# endif
+# ifndef ROOT_TString
+#  include <TString.h>
+# endif
+# ifndef ROOT_TRandom
+#  include <TRandom.h>
+# endif
+# ifndef ROOT_TBrowser
+#  include <TBrowser.h>
+# endif
+# ifndef ROOT_TCanvas
+#  include <TCanvas.h>
+# endif
+# ifndef __CSTDARG__
+#  include <cstdarg>
+# endif
+# ifndef __IOSTREAM__
+#  include <iostream>
+# endif
+#endif
+
+//____________________________________________________________________
+/** @brief Simple geometry 
+    @ingroup geo_geom
+    @code 
+    gSystem->Load("libPhysics.so");
+    gSystem->Load("libGeom.so");
+    gROOT->LoadMacro("GeoGeometry.C+");
+    Geometry g;
+    g.Exec();
+    @endcode 
+ */
+class Geometry 
+{
+public:
+  /** Constructor */
+  Geometry();
+  /** Destructor */
+  virtual ~Geometry() {}
+  /** Initialize  */
+  virtual void Initialize();
+  /** Register  */
+  virtual void Register();
+  /** @e Do-It member function  */
+  virtual void Exec();
+  /** Try to align */
+  virtual void Align();
+  /** Convert detector coordinates to spatial coordinates 
+      @param sector Sector number 
+      @param strip  Strip number
+      @param xyz    Spatial coordinates
+  */
+  virtual void Detector2XYZ(UInt_t sector, UInt_t strip, TVector3& xyz);
+  /** Debug level */
+  static Int_t fgDebug;
+  /** Print debugging messages 
+      @param lvl Acceptance level 
+      @param where Where it happened
+      @param format Message format. */
+  static void Debug(Int_t lvl, const char* where, const char* format, ...) 
+  {
+    if (lvl > fgDebug) return;
+    static char buf[2048];
+    va_list ap;
+    va_start(ap, format);
+    vsnprintf(buf, 2048, format, ap);
+    std::cout << "D: " << where << ": " << buf << std::endl;
+    va_end(ap);
+  }
+protected:
+  /** List of translations */
+  TObjArray*    fMatricies;
+  /** Shape parameter */
+  TVector2     fA;
+  /** Shape parameter */
+  TVector2     fB;
+  /** Shape parameter */
+  TVector2     fC;
+  /** Spacing between sensor and hybrid */
+  Double_t      fSpacer;
+  /** Thickness of aluminum plates in honey comb */
+  Double_t     fAlThickness;
+  /** Width of bonding pads */
+  Double_t      fBondingWidth;
+  /** chip layer thickenss */
+  Double_t      fChipThickness;
+  /** Copper layer thickness */
+  Double_t      fCopperThickness;
+  /** Upper radious */
+  Double_t     fHighR;
+  /** Thickness of honey comb */
+  Double_t     fHoneycombThickness;
+  /** Inner honey comb inner radius  */
+  Double_t     fInnerHoneyHighR;
+  /** Inner honey comb outer radius  */
+  Double_t     fInnerHoneyLowR;
+  /** Z coordinate of inner ring  */
+  Double_t     fInnerZ;
+  /** Length of support legs */
+  Double_t     fLegLength;
+  /** Offset from edge of legs */
+  Double_t     fLegOffset;
+  /** Radius of support legs */
+  Double_t     fLegRadius;
+  /** Inner radius */
+  Double_t     fLowR;
+  /** Spacing between modules */
+  Double_t     fModuleSpacing;
+  /** Number of strps */
+  Int_t                fNStrips;
+  /** Outer honey comb inner radius */
+  Double_t     fOuterHoneyHighR;
+  /** Outer honey comb outer radius */
+  Double_t     fOuterHoneyLowR;
+  /** Z coordinate of outer */
+  Double_t     fOuterZ;
+  /** Thickness of print board */
+  Double_t     fPrintboardThickness;
+  /** Cache of ring depth */
+  Double_t     fRingDepth;
+  /** Thickness of silicon sensor */
+  Double_t     fSiThickness;
+  /** ??? */
+  Double_t     fSpacerHeight;
+  /** Opening angle of sensor */
+  Double_t     fTheta;
+  /** Radius of wafer the sensor is cut out of */
+  Double_t     fWaferRadius;
+  /** Air tracking medium */
+  TGeoMedium*  fAir;
+  /** Aluminum tracking medium */
+  TGeoMedium*  fAl;
+  /** Carbon tracking medium */
+  TGeoMedium*   fCarbon;
+  /** Chip tracking medium */
+  TGeoMedium*   fChip;
+  /** Copper tracking medium */
+  TGeoMedium*   fCopper;
+  /** Kapton tracking medium */
+  TGeoMedium*   fKapton;
+  /** PCB tracking medium */
+  TGeoMedium*  fPCB;
+  /** Plastic tracking medium */
+  TGeoMedium*  fPlastic;
+  /** Active silicon tracking medium */
+  TGeoMedium*  fSi;
+  /** Vacuum  tracking medium */
+  TGeoMedium*  fVacuum;
+  /** Use assemblies */
+  Bool_t        fDoubleAssembly;
+  /** Make a detailed geometry */
+  void MakeDetailed(TGeoVolume* mother);
+};
+
+//____________________________________________________________________
+Int_t Geometry::fgDebug = 10;
+
+//____________________________________________________________________
+Geometry::Geometry()
+  : fMatricies(0),
+    fA( 4.3000,   1.3972),
+    fB(17.2000,   2.1452),
+    fC(15.3327,   4.9819)
+{
+  fBondingWidth           =   0.5000;
+  fWaferRadius            =   6.7000;
+  fSiThickness            =   0.0300;
+  fLowR                   =   4.3000;
+  fHighR                  =  17.2000;
+  fTheta                  =  18.0000;
+  fNStrips                =      512;
+  fRingDepth              =   2.1600;
+  fLegRadius              =   0.5000;
+  fLegLength              =   1.0000;
+  fLegOffset              =   2.0000;
+  fModuleSpacing          =   1.0000;
+  fPrintboardThickness    =   0.1000;
+  fSpacerHeight           =   0.0300;
+
+  fInnerZ                 =  83.4000;
+  fOuterZ                 =  75.2000;
+  fHoneycombThickness     =   1.0000;
+  fAlThickness            =   0.1000;
+  fInnerHoneyLowR         =   5.3000;
+  fInnerHoneyHighR        =  18.2000;
+  fOuterHoneyLowR         =  16.6000;
+  fOuterHoneyHighR        =  29.0000;
+
+  fCopperThickness        =   0.01;
+  fChipThickness          =   0.01;
+  fAlThickness            =   0.10;
+  fHoneycombThickness     =   1.00;
+  fSpacer                 =   0.10;
+  fLegLength              =   1.00;
+}
+
+//____________________________________________________________________
+void 
+Geometry::Exec()
+{
+  Initialize();
+  Register();
+  // gGeoManager->DefaultColors();
+  gGeoManager->ViewLeaves(kFALSE);
+  TCanvas* c = new TCanvas;
+  c->SetFillColor(0);
+  gGeoManager->GetTopVolume()->Draw();
+  new TBrowser("b", "Browser");
+}
+
+//____________________________________________________________________
+void 
+Geometry::Initialize()
+{
+  Double_t mMax  = 10;
+  Double_t mType = 2;
+  
+  // Air 
+  Double_t pAir[]     = { 0., mType, mMax, 1.,  .001, 1., .001, .001 };
+  TGeoMixture* fmdAir = new TGeoMixture("FMD air", 4, .00120479);
+  fmdAir->DefineElement(0, 12.0107,  6., 0.000124);
+  fmdAir->DefineElement(1, 14.0067,  7., 0.755267);
+  fmdAir->DefineElement(2, 15.9994,  8., 0.231781);
+  fmdAir->DefineElement(3, 39.948,  18., 0.012827);
+  fmdAir->SetTransparency('0');
+  fAir = new TGeoMedium("FMD Air", 1, fmdAir,  pAir);
+
+  // Silicon 
+  Double_t pSi[]      = { 1., mType, mMax, 1.,  .001, 1., .001, .001 };
+  TGeoMaterial* fmdSi = new TGeoMaterial("FMD Si", 28.0855, 14, 2.33);
+  fmdSi->SetFillColor(2);
+  fSi = new TGeoMedium("FMD Si", 1, fmdSi,  pSi);
+
+  // Vacumm 
+  Double_t pVacuum[]  = { 0., mType, mMax, 10,  .01, .1, .003, .003 };
+  TGeoMaterial* fmdVacuum = new TGeoMaterial("FMD Vacuum",1e-16,1e-16,1e-16);
+  fmdVacuum->SetTransparency('0');
+  fVacuum = new TGeoMedium("FMD Vacuum", 1, fmdVacuum,pVacuum);
+
+
+  // PCB 
+  Double_t pPCB[]     = { 0., mType, mMax, 1.,  .001, 1., .001, .001 };
+  TGeoMixture* fmdPCB = new TGeoMixture("FMD PCB", 14, 1.8);
+  fmdPCB->DefineElement(0,  28.0855,   14, 0.15144894);
+  fmdPCB->DefineElement(1,  40.078,    20, 0.08147477);
+  fmdPCB->DefineElement(2,  26.981538, 13, 0.04128158);
+  fmdPCB->DefineElement(3,  24.305,    12, 0.00904554);
+  fmdPCB->DefineElement(4,  10.811,     5, 0.01397570);
+  fmdPCB->DefineElement(5,  47.867,    22, 0.00287685);
+  fmdPCB->DefineElement(6,  22.98977,  11, 0.00445114);
+  fmdPCB->DefineElement(7,  39.0983,   19, 0.00498089);
+  fmdPCB->DefineElement(8,  55.845,    26, 0.00209828);
+  fmdPCB->DefineElement(9,  18.9984,    9, 0.00420000);
+  fmdPCB->DefineElement(10, 15.9994,    8, 0.36043788);
+  fmdPCB->DefineElement(11, 12.0107,    6, 0.27529425);
+  fmdPCB->DefineElement(12, 14.0067,    7, 0.01415852);
+  fmdPCB->DefineElement(13,  1.00794,   1, 0.03427566);
+  fmdPCB->SetFillColor(7);
+  fPCB = new TGeoMedium("FMD PCB", 1, fmdPCB,  pPCB);
+
+  // Chip 
+  Double_t pChip[]  = { 0., mType, mMax, 10., .01,  1., .003, .003 };
+  TGeoMixture* fmdChip = new TGeoMixture("FMD Chip", 6, 2.36436);
+  fmdChip->DefineElement(0,  12.0107,   6, 0.039730642);
+  fmdChip->DefineElement(1,  14.0067,   7, 0.001396798);
+  fmdChip->DefineElement(2,  15.9994,   8, 0.01169634);
+  fmdChip->DefineElement(3,   1.00794,  1, 0.004367771);
+  fmdChip->DefineElement(4,  28.0855,  14, 0.844665);
+  fmdChip->DefineElement(5, 107.8682,  47, 0.0981434490);
+  fmdChip->SetFillColor(4);
+  fChip = new TGeoMedium("FMD Chip", 1, fmdChip, pChip);
+
+  // Carbon 
+  Double_t pC[]       = { 0., mType, mMax, 10., .01,  1., .003, .003 };
+  TGeoMaterial* fmdC = new TGeoMaterial("FMD C", 12.011, 6, 2.265);
+  fmdC->SetFillColor(6);
+  fCarbon = new TGeoMedium("FMD C", 1, fmdC, pC);
+  
+  // Kapton (inside of Honeycombs)
+  Double_t pKapton[]  = { 0., mType, mMax, 1.,  .001, 1., .001, .001 };
+  TGeoMaterial* fmdKapton = new TGeoMaterial("FMD Kapton", 12.011,   6.,0.01);
+  fmdKapton->SetFillColor(3);
+  fKapton                 = new TGeoMedium("FMD Kapton",1,fmdKapton,pKapton);
+
+  // Plastic 
+  Double_t pPlastic[] = { 0., mType, mMax, 10., .01,  1., .003, .003 };
+  TGeoMixture* fmdPlastic = new TGeoMixture("FMD Plastic", 2, 1.03);
+  fmdPlastic->DefineElement(0,  1.01,   1, .5);
+  fmdPlastic->DefineElement(1,  12.01,  6, .5);
+  fmdPlastic->SetFillColor(4);
+  fPlastic = new TGeoMedium("FMD Plastic", 1, fmdPlastic,  pPlastic);
+
+  // Aluminium 
+  Double_t pAl[]  = { 0., mType, mMax, 10., .001,  -1., .003, .003 };
+  TGeoMaterial* fmdAl = new TGeoMaterial("FMD Al", 26.981539, 13, 2.7);
+  fmdAl->SetFillColor(3);
+  fAl = new TGeoMedium("FMD Al", 1, fmdAl, pAl);
+
+  // Copper 
+  Double_t pCopper[]  = { 0., mType, mMax, 10., .01,  1., .003, .003 };
+  TGeoMaterial* fmdCopper = new TGeoMaterial("FMD Copper", 63.546,  29.,8.96);
+  fmdCopper->SetFillColor(3);
+  fCopper = new TGeoMedium("FMD Copper",  1, fmdCopper, pCopper);
+}
+
+#define DEGRAD TMath::Pi()/180.
+//____________________________________________________________________
+void 
+Geometry::Detector2XYZ(UInt_t sector, UInt_t strip, TVector3& xyz)
+{
+  UInt_t   mod      = sector / 2;
+  if (!fMatricies) {
+    Warning("Detector2XYZ", "No matricies");
+    return;
+  }
+  TGeoMatrix* m = static_cast<TGeoMatrix*>(fMatricies->At(mod));
+  if (!m) {
+    Warning("Detector2XYZ", "No matrix found for module %d", mod);
+    return;
+  }
+  Debug(10, "Detector2XYZ", "Transforming (%2d,%3d)", sector, strip);
+  Double_t rmax     = fB.Mod();
+  Double_t stripoff = fA.Mod();
+  Double_t dstrip   = (rmax - stripoff) / fNStrips;
+  Double_t r        = (strip + .5) * dstrip + stripoff; // fLowR
+  Double_t theta    = ((sector % 2) - .5) * fTheta;
+  Double_t modThick = (fSiThickness 
+                      + fPrintboardThickness 
+                      + fCopperThickness
+                      + fChipThickness 
+                      + fSpacer);
+  Debug(10,"Detector2XYZ", "Radius %7.3f, angle %7.3f (%f, %f)", r, theta, 
+       fLowR, stripoff);
+  Double_t local[] = {
+    r * TMath::Cos(theta * DEGRAD), 
+    r * TMath::Sin(theta * DEGRAD), 
+    -modThick + fSiThickness / 2
+  };
+  Double_t master[3];
+  Debug(10, "Detector2XYZ", "Local (%7.3f,%7.3f,%7.3f)", 
+       local[0], local[1], local[2]);
+  m->LocalToMaster(local, master);
+  Debug(10, "Detector2XYZ", "Master (%7.3f,%7.3f,%7.3f)", 
+       master[0], master[1], master[2]);
+  xyz.SetXYZ(master[0], master[1], master[2]);
+}
+//____________________________________________________________________
+void
+Geometry::MakeDetailed(TGeoVolume* caveVolume)
+{
+  Info("MakeSimple", "Using a detailed geometry");
+  Double_t xv[6] = { fA.X(), fC.X(), fB.X(),  fB.X(),  fC.X(),  fA.X() };
+  Double_t yv[6] = { fA.Y(), fC.Y(), fB.Y(), -fB.Y(), -fC.Y(), -fA.Y() };
+  Double_t rmax     = fB.Mod();
+  Double_t stripoff = fA.Mod();
+  Double_t dstrip   = (rmax - stripoff) / fNStrips;
+
+  // Double_t hybridThick = (fPrintboardThickness + fCopperThickness 
+  //                        + fChipThickness);
+  // Double_t modThick    = fSiThickness + fSpacer + hybridThick;
+  // Double_t fmdWidth    = (modThick + fHoneycombThickness + fModuleSpacing
+  //                        + fLegLength + fSpacer);
+  
+  // Top
+  // TGeoTube*   fmdShape     = new TGeoTube(fLowR-.1,fHighR+.1,fmdWidth/2+.1);
+  // TGeoPcon*   fmdShape     = new TGeoPcon(0, 360, 2);
+  // fmdShape->DefineSection(0, -fmdWidth / 2 - .1, fLowR-.1, fHighR+.1);
+  // fmdShape->DefineSection(1,  fmdWidth / 2 + .1, fLowR-.1, fHighR+.1);
+  // TGeoVolume* fmdVolume    = new TGeoVolume("FMD1", fmdShape, fVacuum);
+  
+  
+  // Sensor 
+  TGeoXtru* sensorShape = new TGeoXtru(2);
+  sensorShape->DefinePolygon(6, xv, yv);
+  sensorShape->DefineSection(0, - fSiThickness/2);
+  sensorShape->DefineSection(1, fSiThickness/2);
+  TGeoVolume* sensorVolume= new TGeoVolume("FISE",sensorShape,fSi);
+  sensorVolume->SetLineColor(2);
+  // sensorVolume->VisibleDaughters(kFALSE);
+  // sensorVolume->SetVisibility(kTRUE);
+
+  // Virtual volume shape to divide 
+  TGeoTubeSeg* activeShape   = new TGeoTubeSeg(fLowR, rmax, fSiThickness/2, 
+                                              - fTheta, fTheta);
+  TGeoVolume*  activeVolume  = new TGeoVolume("FIAC", activeShape,fSi);
+  activeVolume->SetLineColor(3);
+  sensorVolume->AddNodeOverlap(activeVolume, 0);
+  // activeVolume->SetTransparency(0x3f);
+  TGeoVolume* sectorVolume   = activeVolume->Divide("FISC",2,2,-fTheta,
+                                                   0,0,"N");
+  TGeoVolume* stripVolume    = sectorVolume->Divide("FIST",1,fNStrips,
+                                                   stripoff,dstrip,0,"SX");
+  (void)stripVolume;
+  
+  // Position
+  Double_t x, y, z;
+  
+  // Make PCB volume 
+  for (Int_t i = 0; i < 3; i++) yv[i] -= fBondingWidth;
+  for (Int_t i = 3; i < 6; i++) yv[i] += fBondingWidth;
+  Double_t off = (TMath::Tan(TMath::Pi() * fTheta / 180) * fBondingWidth);
+
+  // PCB layer 
+  TGeoXtru* pcbShape      = new TGeoXtru(2);
+  pcbShape->DefinePolygon(6, xv, yv);
+  pcbShape->DefineSection(0, - fPrintboardThickness/2);
+  pcbShape->DefineSection(1, fPrintboardThickness/2);
+  TGeoVolume* pcbVolume   = new TGeoVolume("FPCB",pcbShape,fPCB);
+  pcbVolume->SetLineColor(4);
+
+  // Copper layer
+  TGeoXtru* cuShape       = new TGeoXtru(2);
+  cuShape->DefinePolygon(6, xv, yv);
+  cuShape->DefineSection(0, - fCopperThickness/2);
+  cuShape->DefineSection(1, fCopperThickness/2);
+  TGeoVolume* cuVolume    = new TGeoVolume("FCOP",cuShape,fCopper);
+  cuVolume->SetLineColor(4);
+
+  // Chip layer
+  TGeoXtru*   chipShape   = new TGeoXtru(2);
+  chipShape->DefinePolygon(6, xv, yv);
+  chipShape->DefineSection(0, - fChipThickness/2);
+  chipShape->DefineSection(1, fChipThickness/2);
+  TGeoVolume* chipVolume = new TGeoVolume("FCHI",chipShape,fChip);
+  chipVolume->SetLineColor(4);
+
+  // Short leg shape 
+  TGeoTube*   shortLegShape  = new TGeoTube(0, fLegRadius, fLegLength / 2);
+  TGeoVolume* shortLegVolume = new TGeoVolume("FISL", shortLegShape, fPlastic);
+  shortLegVolume->SetLineColor(2);
+  // Long leg shape
+  TGeoTube*   longLegShape   = new TGeoTube(0, fLegRadius, 
+                                           (fLegLength + fModuleSpacing) / 2);
+  TGeoVolume* longLegVolume  = new TGeoVolume("FILL", longLegShape, fPlastic);
+  longLegVolume->SetLineColor(2);
+  
+  // Make a front volume 
+  TGeoVolume* frontVolume    =  new TGeoVolumeAssembly("FIFV");
+  z                          =  fPrintboardThickness / 2;
+  frontVolume->AddNode(pcbVolume, 0, new TGeoTranslation(0, 0, z));
+  z                          += (fPrintboardThickness + fCopperThickness) / 2;
+  frontVolume->AddNode(cuVolume, 0, new TGeoTranslation(0, 0, z));
+  z                          += (fCopperThickness + fChipThickness) / 2;
+  frontVolume->AddNode(chipVolume, 0, new TGeoTranslation(0, 0, z));
+  z                          += (fChipThickness + fModuleSpacing 
+                                + fLegLength) / 2;
+  x                          =  fA.X() + fLegOffset + fLegRadius;
+  y                          =  0;
+  frontVolume->AddNode(longLegVolume, 0, new TGeoTranslation(x, y, z));
+  x                          =  fC.X();
+  y                          =  fC.Y() - fLegOffset - fLegRadius - off;
+  frontVolume->AddNode(longLegVolume, 1, new TGeoTranslation(x,y,z));
+  y                          =  -y;
+  frontVolume->AddNode(longLegVolume, 2, new TGeoTranslation(x,y,z));
+
+  // Make a back volume 
+  TGeoVolume* backVolume     =  new TGeoVolumeAssembly("FIBV");
+  z                          =  fPrintboardThickness / 2;
+  backVolume->AddNode(pcbVolume, 0, new TGeoTranslation(0, 0, z));
+  z                          += (fPrintboardThickness + fCopperThickness) / 2;
+  backVolume->AddNode(cuVolume, 0, new TGeoTranslation(0, 0, z));
+  z                          += (fCopperThickness + fChipThickness) / 2;
+  backVolume->AddNode(chipVolume, 0, new TGeoTranslation(0, 0, z));
+  z                          += (fChipThickness + fLegLength) / 2;
+  x                          =  fA.X() + fLegOffset + fLegRadius;
+  y                          =  0;
+  backVolume->AddNode(shortLegVolume, 0, new TGeoTranslation(x, y, z));
+  x                          =  fC.X();
+  y                          =  fC.Y() - fLegOffset - fLegRadius - off;
+  backVolume->AddNode(shortLegVolume, 1, new TGeoTranslation(x,y,z));
+  y                          =  -y;
+  backVolume->AddNode(shortLegVolume, 2, new TGeoTranslation(x,y,z));
+
+#if 1
+  TGeoVolume* ringTopVolume  = new TGeoVolumeAssembly("FITV");
+  TGeoVolume* ringBotVolume  = new TGeoVolumeAssembly("FIBV");
+  Double_t dz = 0;
+#else
+  Double_t dz = (fSiThickness + fSpacer + fModuleSpacing + 
+                fPrintboardThickness + fCopperThickness + 
+                fChipThickness + fLegLength);
+  TGeoTubeSeg* topRingShape = new TGeoTubeSeg(fA.X(), rmax, dz/2, 0, 180);
+  TGeoTubeSeg* botRingShape = new TGeoTubeSeg(fA.X(), rmax, dz/2, 180, 360);
+  TGeoVolume* ringTopVolume = new TGeoVolume("FITV", topRingShape, fAir);
+  TGeoVolume* ringBotVolume = new TGeoVolume("FIBV", botRingShape, fAir);
+#endif
+
+  TGeoVolume* half = ringTopVolume;
+  // Place in mother 
+  for (Int_t i = 0; i < 10; i++) {
+    if (i > 4)  half  = ringBotVolume;
+    Bool_t      front = ((i % 2) == 0);
+    z                 =  -dz/2 + fSiThickness / 2 + (i % 2) * fModuleSpacing;
+    TGeoMatrix* m1    = new TGeoCombiTrans(0,0,z,0);
+    m1->RotateZ((2 * i + 1) * fTheta);
+    half->AddNode(sensorVolume, i, m1);
+    TGeoVolume* vol   = (front ? frontVolume : backVolume);
+    z                 += fSpacer + fSiThickness / 2;
+    TGeoMatrix* m2    = new TGeoCombiTrans(0,0,z,0);
+    m2->RotateZ((2 * i + 1) * fTheta);
+    half->AddNode(vol, i, m2);    
+  }
+
+  TGeoVolume* fmdTopVolume  = new TGeoVolumeAssembly("FMT1");
+  TGeoVolume* fmdBotVolume  = new TGeoVolumeAssembly("FMB1");
+  fmdTopVolume->AddNode(ringTopVolume, 0, new TGeoTranslation(0,0,dz/2));
+  fmdBotVolume->AddNode(ringBotVolume, 0, new TGeoTranslation(0,0,dz/2));
+  
+  // Top of Honeycomb
+  TGeoTubeSeg* hcShape       = new TGeoTubeSeg(fLowR, fHighR,
+                                              fHoneycombThickness / 2,
+                                              0, 180);
+  TGeoVolume*  hcVolume      = new TGeoVolume("F1II", hcShape, fAl);
+  hcVolume->VisibleDaughters(kFALSE);
+  hcVolume->SetLineColor(5);
+  // Air in top of honeycomb
+  TGeoTubeSeg* ihcShape      = new TGeoTubeSeg(fLowR + fAlThickness, 
+                                              fHighR - fAlThickness, 
+                                              (fHoneycombThickness 
+                                               - fAlThickness) / 2,
+                                              0, 180);
+  TGeoVolume*  ihcVolume     = new TGeoVolume("F1IK", ihcShape, fKapton);
+  hcVolume->AddNode(ihcVolume, 0);
+
+  // Add honey comb to mothers. 
+  z = (fSiThickness + fSpacer + fPrintboardThickness + fCopperThickness
+       + fChipThickness + fModuleSpacing + fLegLength + fHoneycombThickness/2);
+  fmdTopVolume->AddNode(hcVolume, 0, new TGeoTranslation(0,0,z));
+  TGeoMatrix* r = new TGeoCombiTrans(0,0,z, 0); r->RotateZ(180);
+  fmdBotVolume->AddNode(hcVolume, 1, r);
+
+  z = 0;
+  caveVolume->AddNode(fmdTopVolume, 1, new TGeoTranslation(0,0,z));
+  caveVolume->AddNode(fmdBotVolume, 1, new TGeoTranslation(0,0,z));
+}
+
+//____________________________________________________________________
+void
+Geometry::Register()
+{
+  // Framework::Task::Register(option);
+  TGeoShape* caveShape = new TGeoBBox(fHighR+10, fHighR+10, fHighR+10);
+  TGeoVolume* caveVolume = new TGeoVolume("Cave", caveShape, fVacuum);
+  gGeoManager->SetTopVolume(caveVolume);
+
+  if (!gGeoManager->GetVolume("FM1T")) MakeDetailed(caveVolume);
+  gGeoManager->CloseGeometry();
+}
+
+//____________________________________________________________________
+void
+Geometry::Align()
+{
+  if (!gGeoManager)                 return;
+  if (!gGeoManager->GetTopVolume()) return;
+  if (!gGeoManager->IsClosed())     return;
+  
+  if (!fMatricies) fMatricies = new TObjArray;
+  TGeoIterator next(gGeoManager->GetTopVolume());
+  TGeoNode* node = 0;
+
+  while ((node = static_cast<TGeoNode*>(next()))) {
+    // node->Print();
+    if (node->GetName()[0] == 'F' && node->GetName()[2] == 'S' && 
+       node->GetName()[3] == 'E') {
+      // We go an FMD module 
+      TString path("/");
+      path.Append(gGeoManager->GetNode(0)->GetName());
+      Int_t nlevel = next.GetLevel();
+      Debug(10, "Exec", "Level is %d", nlevel);
+      for (int lvl = 0; lvl <= nlevel; lvl++) {
+       TGeoNode* p = next.GetNode(lvl);
+       if (!p) continue;
+       Debug(10, "Exec", "Adding '%s' to path '%s'",
+             p->GetName(), path.Data());
+       if (!path.IsNull()) path.Append("/");
+       path.Append(p->GetName());
+      }
+      TGeoPhysicalNode* pnode = gGeoManager->MakePhysicalNode(path.Data());
+      TGeoHMatrix* matrix = new TGeoHMatrix(*node->GetMatrix());
+      TGeoRotation pertub;
+      Double_t angles[] = { gRandom->Uniform(0, 3), 
+                           gRandom->Uniform(0, 3), 
+                           gRandom->Uniform(0, 3) };
+      pertub.RotateX(angles[0]);
+      pertub.RotateY(angles[1]);
+      pertub.RotateZ(angles[1]);
+      *matrix *= pertub;
+      Debug(5, "Exec", "Aliging %s (%f,%f,%f)",
+           pnode->GetName(), angles[0], angles[1], angles[2]);
+      pnode->Align(matrix);
+    }
+  }
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+  
index 21e7432..daedba8 100644 (file)
 #include <TArrayI.h>
 #include <AliRun.h>
 #include <TNtuple.h>
-
+#include <AliDetector.h>
+#include <TFile.h>
+#include <TGeoManager.h>
+#include <TGeoVolume.h>
+#include <TGeoMedium.h>
+#include <TParticle.h>
+#include <THStack.h>
 /** @class Media 
     @brief Media description 
     @ingroup FMD_script
  */
 struct Media : public TNamed
 {
-  TArrayI*        fMeds;
+  TArrayI         fMeds;
   TNtuple*        fCount;
+  TH1D*           fHist;
   Media(const char* name, const char* title) 
     : TNamed(name, title), fMeds(0) 
   {
-    fCount = new TNtuple(GetName(), "E/I:DET:RNG:SEC:STR:X/F:Y:Z:EDEP:PDG/I");
+    fHist  = new TH1D(Form("h%s",name), title, 120, -6, 6);
+    // fHist->SetFillStyle(3001);
+    fCount = new TNtuple(GetName(), GetTitle(), 
+                        "E:DET:RNG:SEC:STR:X:Y:Z:EDEP:PDG");
   }
-  Media(const char* name) 
-    : TNamed(name, "Media information"), fMeds(0), fCount(0)
+  Media(const char* name, TArrayI* a) 
+    : TNamed(name, "Media information"), fMeds(a->fN, a->fArray)
   {
-    AliDetector* det = gAlice->GetDetector(name);
-    if (!det) {
-      Warning("Media", "Detector %s not found in gAlice", name);
-      return;
+    fHist  = new TH1D(Form("h%s",name), GetTitle(), 120, -6, 6);
+    fHist->SetFillStyle(3001);
+    fHist->SetFillColor(fMeds.fN > 0 ? fMeds[0] : 0);
+    fCount = new TNtuple(GetName(), GetTitle(), 
+                        "E:DET:RNG:SEC:STR:X:Y:Z:EDEP:PDG");
+    std::cout << "Media list " << name << ":" << std::endl;
+    for (Int_t i = 0; i < fMeds.fN; i++) {
+      std::cout << " " << fMeds[i] << std::flush;
+      if (fMeds[i] == 0 && i > 0) break;
     }
-    fMeds = det->GetIdtmed();
-    if (!fMeds) {
-      Warning("Media", "No mediums for detector %s found", name);
-      return;
+    std::cout << std::endl;
+  }
+  Bool_t HasMedia(Int_t id) 
+  {
+    if (fMeds.fN <= 0) return kFALSE;
+    for (Int_t i = 0; i < fMeds.fN; i++) {
+      if (fMeds[i] == id) return kTRUE;
+      if (fMeds[i] == 0 && i > 0) break;
     }
-    fCount = new TNtuple(GetName(), "E/I:DET:RNG:SEC:STR:X/F:Y:Z:EDEP:PDG/I");
+    return kFALSE;
+  }
+  void Fill(Int_t ev, AliFMDHit* hit) 
+  {
+    Float_t x[10];
+    x[0] = ev;
+    x[1] = hit->Detector();
+    x[2] = Int_t(hit->Ring());
+    x[3] = hit->Sector();
+    x[4] = hit->Strip();
+    x[5] = hit->X();
+    x[6] = hit->Y();
+    x[7] = hit->Z();
+    x[8] = hit->Edep();
+    x[9] = hit->Pdg();
+    // return;
+    fCount->Fill(x);    
+    // r = TMath::Sqrt(hit->X() * hit->X() + hit->Y() * hit->Y());
+    // if(r!=0) Float_t wgt=1./(2. * 3.1415*r);
+    Double_t r = TMath::Sqrt(hit->X()*hit->X()+hit->Y()*hit->Y());
+    Double_t t = TMath::ATan2(r, hit->Z());
+    Double_t e = -TMath::Log(TMath::Tan(t/2));
+    fHist->Fill(e);
   }
 };
 
@@ -52,6 +93,12 @@ struct Media : public TNamed
 /** @class GetMedia
     @brief Get media where a particle is produced
     @code 
+    aliroot
+    Root> gAlice->Init("$ALICE_ROOT/FMD/Config.C");
+    Root> .x $ALICE_ROOT/FMD/scriptsWriteMedArrays.C
+    Root> gAlice->Run();
+    Root> .q
+    aliroot
     Root> .L Compile.C
     Root> Compile("GetMedia.C")
     Root> GetMedia c
@@ -66,53 +113,99 @@ private:
   TObjArray  fMedia;
   Media*     fOther;
   Media*     fAll;
+  Media*     fPrim;
   Int_t      fEv;
+  THStack*   fStack;
   TFile*     fOutput;
 public:
   //__________________________________________________________________
   GetMedia(const char* modlist="FMD:ITS:BODY:ABSO:START:PIPE", 
           const char* output="media.root") 
+    :  fModList(modlist)
   { 
+    AddLoad(kGeometry);
+    AddLoad(kKinematics);
+
     fOutput = TFile::Open(output, "RECREATE");
-    fOther  = new Media("other", "Unknown media"),
-    fAll    = new Media("All", "All media")
+    fStack  = new THStack("summed", "Summed hits");
+    fOther  = AddMedia("other", "Unknown media");
+    fAll    = AddMedia("All", "All media");
+    fPrim   = AddMedia("Primary", "Primary particles");
     fEv     = 0;
-    AddLoad(kKinematics);
+
+    
+    fAll->fHist->SetFillColor(0);
+    fAll->fHist->SetFillStyle(3004);
+    fPrim->fHist->SetFillColor(1);
+    fPrim->fHist->SetFillStyle(3001);
+    fStack->Add(fPrim->fHist);
   }
   //__________________________________________________________________
-  Media* FindMedia(Int_t med) 
+  Media* AddMedia(const char* name, const char* title, TArrayI* a=0)
   {
-    TIter next(&fMedia);
     Media* media = 0;
-    while ((media == static_cast<Media*>(next()))) {
-      if (!media->fMeds) continue;
-      TArrayI& ids = *(media->fMeds);
-      for (Int_t i = 0; i < ids.fN; i++) 
-       if (med == ids[i]) return media;
+    if (!a) media = new Media(name, title);
+    else    media = new Media(name, a);
+    if (a)  {
+      fMedia.Add(media);
+      fStack->Add(media->fHist);
+    }
+    return media;
+  }
+  //__________________________________________________________________
+  Media* FindMedia(Int_t med) 
+  {
+    for (Int_t i = 0; i < fMedia.GetEntries(); i++) {
+      Media* media = static_cast<Media*>(fMedia.At(i));
+      if (!media) continue;
+      if (media->HasMedia(med)) return media;
     }
     return 0;
   }
   //__________________________________________________________________
   Bool_t Init()
   {
+    Bool_t ret = AliFMDInputHits::Init();
     if (!gGeoManager) {
-      Error("GetMedia", "No geometry defined - make sure you have that");
+      Error("Init", "No geometry defined - make sure you have that");
       return kFALSE;
     }
-    if (!gAlice) {
-      Error("GetMedia", "gAlice not defined");
+    if (!fRun) {
+      Error("Init", "gAlice not defined");
       return kFALSE;
     }
     Int_t now = 0;
     Int_t colon;
+    TFile* file = TFile::Open("medid.root", "READ");
+    if (!file) {
+      Error("Init", "couldn't open medid.root");
+      return kFALSE;
+    }
+    fOutput->cd();
     while ((colon = fModList.Index(":", now)) >= 0) {
-      fMedia.Add(new Media(fModList(now, colon-now)));
+      TString d(fModList(now, colon-now));
       now = colon+1;
+      TArrayI* a = 0;
+      file->GetObject(d.Data(), a);
+      if (!a) {
+       Warning("Init", "No medium array for %s", d.Data());
+       continue;
+      }
+      AddMedia(d.Data(),0, a);
     }
-    if (now < fModList.Length()) 
-      fMedia.Add(new Media(now, fModList.Length()-now));
+    if (now < fModList.Length()) {
+      colon = fModList.Length();
+      TString d(fModList(now, colon-now));
+      TArrayI* a = 0;
+      file->GetObject(d.Data(), a);
+      if (!a) 
+       Warning("Init", "No medium array for %s", d.Data());
+      else 
+       AddMedia(d.Data(),0, a);
+    }
+    file->Close();
     if (fMedia.GetEntries() <= 0) return kFALSE;
-    return AliFMDInputHits::Init();
+    return ret;
   }  
   //__________________________________________________________________
   /** Begining of event
@@ -131,43 +224,62 @@ public:
                << track << ")" << std::endl;
       return kFALSE;
     }
+    fAll->Fill(fEv, hit);
+    if (track->IsPrimary()) { 
+      fPrim->Fill(fEv, hit);
+      return kTRUE;
+    }
+      
     // Get production vertex 
     Double_t vx = track->Vx();
     Double_t vy = track->Vy();
     Double_t vz = track->Vz();
-    
-    fAll->fCounts->Fill(fEv, hit->Detector(),Int_t(hit->Ring()),
-                       hit->Sector(), hit->Strip(), 
-                       hit->X(), hit->Y(), hit->Z(), hit->Edep(), 
-                       hit->Pdg());
     // Get node 
     TGeoNode* prodNode = gGeoManager->FindNode(vx,vy,vz);
-    if (!prodNode) return kTRUE;
-    
+    if (!prodNode) { 
+      Warning("ProcessHit", "didn't find a node for production vertex");
+      return kTRUE;
+    }
+    //Info("ProcessHit", "Production vertex in node %s", prodNode->GetName());
+
     // Get volume 
     TGeoVolume* prodVol = prodNode->GetVolume();
-    if (!prodVol) return kTRUE;
+    if (!prodVol) { 
+      Warning("ProcessHit", "didn't find a volume for production vertex");
+      return kTRUE;
+    }
+    //Info("ProcessHit", "Production vertex in volume %s",prodVol->GetName());
     
     // Med medium 
     TGeoMedium* prodMed = prodVol->GetMedium();
-    if (!prodMed) return kTRUE;
+    if (!prodMed) { 
+      Warning("ProcessHit", "didn't find a medium for production vertex");
+      return kTRUE;
+    }
+    // Info("ProcessHit", "Production vertex in medium %s %d", 
+    //      prodMed->GetName(), prodMed->GetId());
     
-    Media* media = FindMedia(prodMed->GetUniqueID());
-    if (media) media = fOther; 
     
-    // r = TMath::Sqrt(hit->X() * hit->X() + hit->Y() * hit->Y());
-    // if(r!=0) Float_t wgt=1./(2. * 3.1415*r);
-    media->fCounts->Fill(fEv, hit->Detector(),Int_t(hit->Ring()),
-                        hit->Sector(), hit->Strip(), 
-                        hit->X(), hit->Y(), hit->Z(), hit->Edep(), 
-                        hit->Pdg());
+    Media* media = FindMedia(prodMed->GetId());
+    if (!media) { 
+      Warning("ProcessHit", "Media not found %s (%d)", 
+             prodMed->GetName(), prodMed->GetId());
+      media = fOther; 
+    }
+    // Info("ProcessHit", "Adding to %s", media->GetName());
+    media->Fill(fEv, hit);
     return kTRUE;
   }
   //__________________________________________________________________
   Bool_t Finish()
   {
+    fOutput->cd();
+    fStack->Write();
+    fStack->Draw();
     fOutput->Write();
     fOutput->Close();
+    fOutput = 0;
+    fMedia.Delete();
     return kTRUE;
   }
 };
diff --git a/FMD/scripts/NodeGeometry.C b/FMD/scripts/NodeGeometry.C
new file mode 100644 (file)
index 0000000..19c1ca1
--- /dev/null
@@ -0,0 +1,363 @@
+//____________________________________________________________________
+//
+//
+// $Id$
+//
+// Script I used for rapid prototyping of the FMD3 geometry - in
+// particular the support cone 
+//
+/** @defgroup node_geom Simple geometry
+    @ingroup FMD_script
+*/
+#include <TGeometry.h>
+#include <TNode.h>
+#include <TXTRU.h>
+#include <TTUBE.h>
+#include <TTUBS.h>
+#include <TPCON.h>
+#include <TBRIK.h>
+#include <TCanvas.h>
+#include <vector>
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+
+//____________________________________________________________________
+/** @brief A 2D point
+    @ingroup node_geom
+ */
+struct point_t
+{
+  point_t(double x=0, double y=0) : first(x), second(y) {}
+  double first;
+  double second;
+};
+
+//____________________________________________________________________
+/** @brief Shape of a ring
+    @ingroup node_geom
+ */
+struct Ring 
+{
+  // typedef std::pair<double,double> point_t;
+  typedef std::vector<point_t>     points_t;
+  /** Constructor 
+      @param rL         Lower radius
+      @param rH         Higer radius
+      @param theta      Opening angle
+      @param waferR     Wafer radius
+      @param siThick    Silicon thickness 
+      @param staggering Staggering of modules */
+  Ring(double rL, double rH, double theta,  double waferR, 
+       double siThick, double staggering)
+    : fStaggering(staggering),
+      fInnerRadius(rL), 
+      fOuterRadius(rH), 
+      fAngle(theta), 
+      fRadius(waferR),
+      fThickness(siThick), 
+      fVerticies(6)
+  {
+    double tan_theta  = tan(fAngle * TMath::Pi() / 180.);
+    double tan_theta2 = pow(tan_theta,2);
+    double r2         = pow(fRadius,2);
+    double ir2        = pow(fInnerRadius,2);
+    double or2        = pow(fOuterRadius,2);
+    double y_A        = tan_theta * fInnerRadius;
+    double x_D        = fInnerRadius + sqrt(r2 - tan_theta2 * ir2);
+    double x_D2       = pow(x_D,2);
+    double y_B        = sqrt(r2 - or2 + 2 * fOuterRadius * x_D - x_D2);
+    double x_C        = ((x_D + sqrt(-tan_theta2 * x_D2 
+                                    + r2 * (1 + tan_theta2))) 
+                        / (1 + tan_theta2));
+    double y_C        = tan_theta * x_C;
+    
+    fVerticies[0] = point_t(fInnerRadius,  y_A);
+    fVerticies[1] = point_t(x_C,           y_C);
+    fVerticies[2] = point_t(fOuterRadius,  y_B);
+    fVerticies[3] = point_t(fOuterRadius, -y_B);
+    fVerticies[4] = point_t(x_C,          -y_C);
+    fVerticies[5] = point_t(fInnerRadius, -y_A);
+  }
+  /** Destructor */
+  virtual ~Ring() 
+  {
+    fVerticies.clear();
+  }
+  /** Create a shape 
+      @return pointer to new shape  */
+  TShape* CreateShape()
+  {
+    std::cout << "Creating Module shape" << std::flush;
+    TXTRU* moduleShape = new TXTRU("Module","Module", "", 6, 2);
+    for (Int_t i = 0; i  < 6; i++) {
+      std::cout << "." << std::flush;
+      point_t& p = fVerticies[i];
+      moduleShape->DefineVertex(i, p.first, p.second);
+    }
+    moduleShape->DefineSection(0, -fThickness/2, 1, 0, 0);
+    moduleShape->DefineSection(1,  fThickness/2, 1, 0, 0);
+    std::cout << std::endl;
+    return (TShape*)moduleShape;
+  }
+  /** Create a node that represents a ring. 
+      @return Node */
+  TNode* CreateRing(const char* name, double z) 
+  {
+    std::cout << "Creating Ring node for " << name << std::flush;
+    double bredth = fStaggering + fThickness;
+    TShape* ringShape   = new TTUBE(Form("%sShape", name), "Ring Shape", 
+                                   "", fInnerRadius, 
+                                   fOuterRadius,bredth/2);
+    TNode*  ringNode    = new TNode(Form("%sNode", name), "Ring Node", 
+                                   ringShape, 0, 0, z+bredth/2, 0);
+    TShape* moduleShape = CreateShape();
+    Int_t n = Int_t(360 / 2 / fAngle);
+    for (Int_t i = 0; i < n; i++) {
+      std::cout << "." << std::flush;
+      ringNode->cd();
+      Double_t theta  = 2  * fAngle * i;
+      Double_t z      = -(bredth+fThickness)/2+(i%2?0:fStaggering);
+      TRotMatrix* rot = new TRotMatrix(Form("%sRotation%02d", name, i), 
+                                      "Rotation", 90, theta, 90, 
+                                      fmod(90 + theta, 360), 0, 0);
+      TNode* moduleNode = new TNode(Form("%sModule%02d", name, i), 
+                                   "Module", moduleShape, 0, 0, z,
+                                   rot);
+      moduleNode->SetFillColor(2);
+      moduleNode->SetLineColor(2);
+      moduleNode->SetLineWidth(2);
+    }
+    std::cout << std::endl;
+    ringNode->SetVisibility(0);
+    return ringNode;
+  }
+  double fStaggering;
+  /** Inner radius */
+  double fInnerRadius;
+  /** Outer radius */
+  double fOuterRadius;
+  /** Opening angle (in degrees) */
+  double fAngle;
+  /** Radius (in centimeters) */
+  double fRadius;
+  /** Thickness */
+  double fThickness;
+  /** List of verticies */
+  points_t fVerticies;
+};
+
+//____________________________________________________________________
+/** @brief Shape of a detector
+    @ingroup node_geom
+ */
+struct Detector
+{
+  /** Constructor 
+      @param id 
+      @param inner 
+      @param outer */
+  Detector(Ring* inner, double iZ, Ring* outer=0, double oZ=0) 
+    : fInner(inner), fInnerZ(iZ), fOuter(outer), fOuterZ(oZ)
+  {}
+  /** Destructor */
+  virtual ~Detector() {}
+  /** Create rings */
+  virtual void CreateRings() 
+  {
+    if (fInner) fInner->CreateRing("inner", fInnerZ);
+    if (fOuter) fOuter->CreateRing("outer", fOuterZ);
+  }
+  /** Create a node that represents the support */
+  virtual void CreateSupport(double) { }
+  /** Pointer to inner ring */
+  Ring* fInner;
+  /** Position in z of inner ring */ 
+  double fInnerZ;
+  /** Pointer to outer ring */
+  Ring* fOuter;
+  /** Position in z of inner ring */ 
+  double fOuterZ;
+};
+
+//____________________________________________________________________
+/** @brief FMD3 simple node geometry 
+    @ingroup node_geom
+ */
+struct FMD3 : public Detector
+{
+  /** Constructor 
+      @param inner Inner ring representation  
+      @param outer Outer ring representation */
+  FMD3(Ring* inner, Ring* outer)
+    : Detector(inner, -62.8,outer, -75.2)
+  {
+    fNoseRl     = 5.5;
+    fNoseRh     = 6.7;
+    fNoseDz     = 2.8 / 2;
+    fNoseZ      = -46;
+    fConeL      = 30.9;
+    fBackRl     = 61 / 2;
+    fBackRh     = 66.8 /2;
+    fBackDz     = 1.4 / 2;
+    fBeamDz     = .5 / 2;
+    fBeamW      = 6;
+    fFlangeR    = 49.25;
+  }
+  virtual ~FMD3() {}
+  void CreateRings()
+  {
+    double zdist      = fConeL - 2 * fBackDz - 2 * fNoseDz;
+    double tdist      = fBackRh - fNoseRh;
+    double alpha      = tdist / zdist;
+    double x, rl, rh, z;
+    z  = fNoseZ - fConeL / 2;
+    TPCON* fmd3Shape = new TPCON("fmd3Shape", "FMD 3 Shape", "", 0, 360, 7);
+    x  = fNoseZ;
+    rl = fNoseRl;
+    rh = fNoseRh;
+    fmd3Shape->DefineSection(0, x - z, rl, rh);
+    x  = fNoseZ-2*fNoseDz;
+    fmd3Shape->DefineSection(1, x - z, rl, rh);
+    x  = fInnerZ - fInner->fStaggering - fInner->fThickness; 
+    rl = fInner->fInnerRadius;
+    rh = fNoseRh + alpha * TMath::Abs(x-fNoseZ + 2 * fNoseDz);
+    fmd3Shape->DefineSection(2, x - z, rl, rh);
+    x  = fOuterZ;
+    rl = fOuter->fInnerRadius;
+    rh = fBackRh;
+    fmd3Shape->DefineSection(3, x - z, rl, rh);
+    x  = fNoseZ - zdist - 2 * fNoseDz;
+    rl = fOuter->fInnerRadius;
+    rh = fBackRh;
+    fmd3Shape->DefineSection(4, x - z, rl, rh);
+    x  = fNoseZ - zdist - 2 * fNoseDz;
+    rl = fOuter->fInnerRadius;
+    rh = fFlangeR;
+    fmd3Shape->DefineSection(5, x - z, rl, rh);
+    x  = fNoseZ - fConeL;
+    rl = fOuter->fInnerRadius;
+    rh = fFlangeR;
+    fmd3Shape->DefineSection(6, x - z, rl, rh);
+
+    TNode* fmd3Node = new TNode("fmd3Node", "FMD3 Node", fmd3Shape,0,0,z,0);
+    fmd3Node->SetLineColor(11);
+    fmd3Node->SetFillColor(11);
+    fmd3Node->SetVisibility(1);
+    fmd3Node->cd();
+    if (fInner) fInner->CreateRing("inner", fInnerZ-z);
+    fmd3Node->cd();
+    if (fOuter) fOuter->CreateRing("outer", fOuterZ-z);
+    fmd3Node->cd();
+    CreateSupport(fNoseZ - z);
+  }
+  
+  /** Create support volumes  */
+  void CreateSupport(double noseZ) 
+  {
+    TShape* noseShape = new TTUBE("noseShape", "Nose Shape", "", 
+                                 fNoseRl, fNoseRh, fNoseDz);
+    TNode*  noseNode  = new TNode("noseNode", "noseNode", noseShape, 
+                                 0, 0, noseZ - fNoseDz, 0);
+    noseNode->SetLineColor(0);
+    double  zdist = fConeL - 2 * fBackDz - 2 * fNoseDz;
+    double  tdist = fBackRh - fNoseRh;
+    double  beamL = TMath::Sqrt(zdist * zdist + tdist * tdist);
+    double  theta = -TMath::ATan2(tdist, zdist);
+    TShape* backShape = new TTUBE("backShape", "Back Shape", "", 
+                                 fBackRl, fBackRh, fBackDz);
+    TNode*  backNode  = new TNode("backNode", "backNode", backShape, 
+                                 0, 0, noseZ-2*fNoseDz-zdist-fBackDz, 0);
+    backNode->SetLineColor(0);
+    TShape* beamShape = new TBRIK("beamShape", "beamShape", "", 
+                                 fBeamDz, fBeamW / 2 , beamL / 2);
+    Int_t    n = 8;
+    Double_t r = fNoseRl + tdist / 2;
+    for (Int_t i = 0; i < n; i++) {
+      Double_t phi   = 360. / n * i;
+      Double_t t     = 180. * theta / TMath::Pi();
+      TRotMatrix* beamRotation = new TRotMatrix(Form("beamRotation%d", i), 
+                                               Form("beamRotation%d", i),
+                                               180-t,phi,90,90+phi,t,phi);
+      TNode* beamNode = new TNode(Form("beamNode%d", i), 
+                                 Form("beamNode%d", i), beamShape, 
+                                 r * TMath::Cos(phi / 180 * TMath::Pi()), 
+                                 r * TMath::Sin(phi / 180 * TMath::Pi()), 
+                                 noseZ-2*fNoseDz-zdist/2, beamRotation);
+      beamNode->SetLineColor(0);
+    }
+    Double_t flangel    = (fFlangeR - fBackRh) / 2;
+    TShape* flangeShape = new TBRIK("flangeShape", "FlangeShape", "", 
+                                   flangel, fBeamW / 2, fBackDz);
+    n = 4;
+    r = fBackRh + flangel;
+    for (Int_t i = 0; i < n; i++) {
+      Double_t phi = 360. / n * i + 180. / n;
+      TRotMatrix* flangeRotation = new TRotMatrix(Form("flangeRotation%d", i),
+                                                 Form("Flange Rotation %d",i),
+                                                 90,phi,90,90+phi,0,0);
+      TNode* flangeNode = new TNode(Form("flangeNode%d", i), 
+                                   Form("flangeNode%d", i), 
+                                   flangeShape,
+                                   r * TMath::Cos(phi / 180 * TMath::Pi()), 
+                                   r * TMath::Sin(phi / 180 * TMath::Pi()),
+                                   noseZ-2*fNoseDz-zdist-fBackDz, 
+                                   flangeRotation);
+      flangeNode->SetLineColor(0);
+                                 
+    }
+  }
+  /** Nose inner radius */
+  double  fNoseRl;
+  /** Nose outer radius */
+  double  fNoseRh;
+  /** Nose depth */
+  double  fNoseDz;
+  /** Nose start position */
+  double  fNoseZ;
+  /** Length of whole support structure */
+  double  fConeL;
+  /** Inner radius of back ring */
+  double  fBackRl;
+  /** Outer radius of back ring */
+  double  fBackRh;
+  /** Thickness of back ring */
+  double  fBackDz;
+  /** Thickness of beams */
+  double  fBeamDz;
+  /** Width of beams */
+  double  fBeamW;
+  /** Ending radius of flanges */
+  double  fFlangeR;
+};
+
+//____________________________________________________________________
+/** @brief Create a node geometry 
+    @ingroup node_geom
+    @code 
+    .x NodeGeometry.C++
+    @endcode 
+ */
+void
+NodeGeometry() 
+{
+  TGeometry* geometry = new TGeometry("geometry","geometry");
+  TShape* topShape = new TBRIK("topShape", "topShape", "", 40, 40, 150);
+  TNode* topNode = new TNode("topNode", "topNode", topShape, 0, 0, 0, 0);
+  topNode->SetVisibility(0);
+  topNode->cd();
+  
+  Ring inner( 4.3, 17.2, 18, 13.4 / 2, .03, 1);
+  Ring outer(15.6, 28.0,  9, 13.4 / 2, .03, 1);
+  FMD3 fmd3(&inner, &outer);
+  fmd3.CreateRings();
+  
+  TCanvas* c = new TCanvas("c", "c", 800, 800);
+  c->SetFillColor(1);
+  geometry->Draw();
+  // c->x3d("ogl");
+}
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/scripts/SimpleGeometry.C b/FMD/scripts/SimpleGeometry.C
deleted file mode 100644 (file)
index 4df0ae7..0000000
+++ /dev/null
@@ -1,330 +0,0 @@
-//____________________________________________________________________
-//
-//
-// $Id$
-//
-// Script I used for rapid prototyping of the FMD3 geometry - in
-// particular the support cone 
-//
-/** @defgroup simple_geom Simple geometry
-    @ingroup FMD_script
-*/
-/** Calculate wafer parameters 
-    @param dl      Lower diameter
-    @param dh      Higer diameter
-    @param theta   Opening angle
-    @param r       Wafer radius
-    @param verbose Maybe be verbose
-    @return  Corners of shape
-    @ingroup simple_geom
-*/
-//____________________________________________________________________
-TObjArray*
-waferParameters(double dl, double dh, double theta, double r, 
-               Bool_t verbose=kFALSE)
-{      
-  double tan_theta  = tan(theta * TMath::Pi() / 180.);
-  double tan_theta2 = pow(tan_theta,2);
-  double r2         = pow(r,2);
-  double y_A        = tan_theta * dl;
-  double x_D        = dl + sqrt(r2 - tan_theta2 * pow(dl,2));
-  double x_D_2      = dl - sqrt(r2 - tan_theta2 * pow(dl,2));
-  
-  double y_B       = sqrt(r2 - pow(dh,2) + 2 * dh * x_D - pow(x_D,2));
-  double x_C       = (x_D + sqrt(-tan_theta2 * pow(x_D,2) + r2 
-                                + r2 * tan_theta2)) / (1 + tan_theta2);
-  double y_C       = tan_theta * x_C;
-
-  if (verbose) {
-    cout << "A: (" << dl << "," << y_A << ")" << endl;
-    cout << "B: (" << dh << "," << y_B << ")" << endl;
-    cout << "C: (" << x_C << "," << y_C << ")" << endl;
-    cout << "D: (" << x_D << ",0)" << endl;
-    
-    cout << "Recentred at D:"  << endl;
-    cout << "A: (" << dl - x_D  << "," << y_A << ")" << endl;
-    cout << "B: (" << dh - x_D  << "," << y_B << ")" << endl;
-    cout << "C: (" << x_C - x_D << "," << y_C << ")" << endl;
-  }
-  
-  TObjArray* verticies = new TObjArray(6);
-  verticies->AddAt(new TVector2(dl,   y_A), 0);
-  verticies->AddAt(new TVector2(x_C,  y_C), 1);
-  verticies->AddAt(new TVector2(dh,   y_B), 2);
-  verticies->AddAt(new TVector2(dh,  -y_B), 3);
-  verticies->AddAt(new TVector2(x_C, -y_C), 4);
-  verticies->AddAt(new TVector2(dl,  -y_A), 5);
-  
-  return verticies;
-}
-
-/** Create the sensor shape 
-    @param name 
-    @param rl 
-    @param rh 
-    @param th 
-    @param r 
-    @param dz 
-    @ingroup simple_geom
-    @return  */
-//____________________________________________________________________
-TShape* 
-createModuleShape(const Char_t* name, double rl, double rh, double th, 
-                 double r, double dz) 
-{
-  std::cout << "Creating Module shape for " << name << std::flush;
-  // TShape* virtualShape = new TTUBS(Form("%sVirtual", name), 
-  //                                  Form("Virtual %s", name), 
-  //                                  "", rl, rh, 1, -th, th);
-  
-  TObjArray* v = waferParameters(rl, rh, th, r);
-  TXTRU* moduleShape = new TXTRU(Form("%sModule", name),
-                                Form("Module %s", name), 
-                                "", 6, 2);
-  for (Int_t i = 0; i  < 6; i++) {
-    std::cout << "." << std::flush;
-    TVector2* vv = static_cast<TVector2*>(v->At(i));
-    moduleShape->DefineVertex(i, vv->X(), vv->Y());
-  }
-  moduleShape->DefineSection(0, -dz, 1, 0, 0);
-  moduleShape->DefineSection(1,  dz, 1, 0, 0);
-  std::cout << std::endl;
-
-  return (TShape*)moduleShape;
-}
-
-/** 
-    @param name 
-    @param rl 
-    @param rh 
-    @param th 
-    @param siThick 
-    @param waferR 
-    @param staggering 
-    @param z 
-    @ingroup simple_geom
-    @return  
-*/
-//____________________________________________________________________
-TNode* 
-createRing(const char* name, double rl, double rh, double th, 
-          double siThick, double waferR, double staggering, double z) 
-{
-  std::cout << "Creating Ring node for " << name << std::flush;
-  TShape* ringShape   = new TTUBE(Form("%sShape", name), "Ring Shape", 
-                                 "", rl, rh, staggering + siThick);
-  TNode*  ringNode    = new TNode(Form("%sNode", name), "Ring Node", 
-                                 ringShape, 0, 0, z, 0);
-  TShape* moduleShape = createModuleShape(name, rl, rh, th, waferR, siThick);
-  Int_t n = 360 / 2 / th;
-  for (Int_t i = 0; i < n; i++) {
-    std::cout << "." << std::flush;
-    ringNode->cd();
-    Double_t theta = 2  * th * i;
-    Double_t z     = (i % 2 ? 0 : staggering);
-    TRotMatrix* rot = new TRotMatrix(Form("%sRotation%02d", name, i), 
-                                    "Rotation", 90, theta, 90, 
-                                    fmod(90 + theta, 360), 0, 0);
-    TNode* moduleNode = new TNode(Form("%sModule%02d", name, i), 
-                                 "Module", moduleShape, 0, 0, z,
-                                 rot);
-    moduleNode->SetFillColor(5);
-    moduleNode->SetLineColor(5);
-    moduleNode->SetLineWidth(2);
-  }
-  std::cout << std::endl;
-  ringNode->SetVisibility(0);
-  return ringNode;
-}
-
-//____________________________________________________________________
-/** 
-    @param noseRl 
-    @param noseRh 
-    @param noseDz 
-    @param noseZ 
-    @param backRl 
-    @param backRh 
-    @param backDz 
-    @param coneL 
-    @param beamW 
-    @param beamDz 
-    @param flangeR 
-    @return  
-    @ingroup simple_geom
-*/
-TNode*
-createSupport(double noseRl, double noseRh, double noseDz, double noseZ, 
-             double backRl, double backRh, double backDz, double coneL, 
-             double beamW, double beamDz,  double flangeR) 
-{
-  TShape* noseShape = new TTUBE("noseShape", "Nose Shape", "", 
-                               noseRl, noseRh, noseDz);
-  TNode*  noseNode  = new TNode("noseNode", "noseNode", noseShape, 
-                               0, 0, noseZ - noseDz, 0);
-  noseNode->SetLineColor(0);
-  
-  Double_t zdist = coneL - 2 * backDz - 2 * noseDz;
-  Double_t tdist = backRh - noseRh;
-  Double_t beamL = TMath::Sqrt(zdist * zdist + tdist * tdist);
-  Double_t theta = -TMath::ATan2(tdist, zdist);
-  
-
-  TShape* backShape = new TTUBE("backShape", "Back Shape", "", 
-                                backRl, backRh, backDz);
-  TNode*  backNode  = new TNode("backNode", "backNode", backShape, 
-                               0, 0, noseZ - 2 * noseDz - zdist - backDz, 0);
-  backNode->SetLineColor(0);
-
-
-  TShape* beamShape = new TBRIK("beamShape", "beamShape", "", 
-                               beamDz, beamW / 2 , beamL / 2);
-  Int_t    n = 8;
-  Double_t r = noseRl + tdist / 2;
-  for (Int_t i = 0; i < n; i++) {
-    Double_t phi   = 360. / n * i;
-    Double_t t     = 180. * theta / TMath::Pi();
-    TRotMatrix* beamRotation = new TRotMatrix(Form("beamRotation%d", i), 
-                                               Form("beamRotation%d", i),
-                                               180-t, phi, 90, 90+phi, 
-                                               t, phi);
-    TNode* beamNode = new TNode(Form("beamNode%d", i), 
-                               Form("beamNode%d", i), beamShape, 
-                               r * TMath::Cos(phi / 180 * TMath::Pi()), 
-                               r * TMath::Sin(phi / 180 * TMath::Pi()), 
-                               noseZ - 2 * noseDz - zdist / 2,  beamRotation);
-    beamNode->SetLineColor(0);
-  }
-  
-  Double_t flangel = (flangeR - backRh) / 2;
-  TShape* flangeShape = new TBRIK("flangeShape", "FlangeShape", "", 
-                                 flangel, beamW / 2, backDz);
-  n = 4;
-  r = backRh + flangel;
-  for (Int_t i = 0; i < n; i++) {
-    Double_t phi = 360. / n * i + 180. / n;
-    TRotMatrix* flangeRotation = new TRotMatrix(Form("flangeRotation%d", i),
-                                               Form("Flange Rotation %d", i),
-                                               90, phi, 90, 90+phi, 0, 0);
-    TNode* flangeNode = new TNode(Form("flangeNode%d", i), 
-                                 Form("flangeNode%d", i), 
-                                 flangeShape,
-                                 r * TMath::Cos(phi / 180 * TMath::Pi()), 
-                                 r * TMath::Sin(phi / 180 * TMath::Pi()),
-                                 noseZ - 2 * noseDz - zdist - backDz, 
-                                 flangeRotation);
-    flangeNode->SetLineColor(0);
-                                 
-  }
-  return 0;
-}
-
-                                
-
-
-//____________________________________________________________________
-/** 
-    @ingroup simple_geom
- */
-void
-SimpleGeometry() 
-{
-  TGeometry* geometry = new TGeometry("geometry","geometry");
-  TShape* topShape = new TBRIK("topShape", "topShape", "", 40, 40, 150);
-  TNode* topNode = new TNode("topNode", "topNode", topShape, 0, 0, 0, 0);
-  topNode->SetVisibility(0);
-  topNode->cd();
-  
-  Double_t waferR     = 13.4 / 2;
-  Double_t siThick    = .03;
-  Double_t staggering = 1;
-  Double_t innerRh    = 17.2;
-  Double_t innerRl    = 4.3;
-  Double_t innerTh    = 18;
-  Double_t innerZ     = -62.8;
-  Double_t outerRh    = 28;
-  Double_t outerRl    = 15.6;
-  Double_t outerTh    = 9;
-  Double_t outerZ     = -75.2;
-  Double_t noseRl     = 5.5;
-  Double_t noseRh     = 6.7;
-  Double_t noseDz     = 2.8 / 2;
-  Double_t noseZ      = -46;
-  Double_t coneL      = 30.9;
-  Double_t backRl     = 61 / 2;
-  Double_t backRh     = 66.8 /2;
-  Double_t backDz     = 1.4 / 2;
-  Double_t beamDz     = .5 / 2;
-  Double_t beamW      = 6;
-  Double_t flangeR    = 49.25;
-
-  Double_t zdist      = coneL - 2 * backDz - 2 * noseDz;
-  Double_t tdist      = backRh - noseRh;
-  Double_t alpha      = tdist / zdist;
-
-  Double_t x, rl, rh, z;
-  z  = noseZ - coneL / 2;
-  TPCON* fmd3Shape = new TPCON("fmd3Shape", "FMD 3 Shape", "", 0, 360, 7);
-  
-  x  = noseZ;
-  rl = noseRl;
-  rh = noseRh;
-  fmd3Shape->DefineSection(0, x - z, rl, rh);
-
-  x  = noseZ-2*noseDz;
-  fmd3Shape->DefineSection(1, x - z, rl, rh);
-
-  x  = innerZ - staggering - siThick; 
-  rl = innerRl;
-  rh = noseRh + alpha * TMath::Abs(x-noseZ + 2 * noseDz);
-  fmd3Shape->DefineSection(2, x - z, rl, rh);
-
-  x  = outerZ;
-  rl = outerRl;
-  rh = backRh;
-  fmd3Shape->DefineSection(3, x - z, rl, rh);
-
-  x  = noseZ - zdist - 2 * noseDz;
-  rl = outerRl;
-  rh = backRh;
-  fmd3Shape->DefineSection(4, x - z, rl, rh);
-
-  x  = noseZ - zdist - 2 * noseDz;
-  rl = outerRl;
-  rh = flangeR;
-  fmd3Shape->DefineSection(5, x - z, rl, rh);
-
-  x  = noseZ - coneL;
-  rl = outerRl;
-  rh = flangeR;
-  fmd3Shape->DefineSection(6, x - z, rl, rh);
-
-  TNode* fmd3Node = new TNode("fmd3Node", "FMD3 Node", fmd3Shape, 
-                             0, 0, z, 0);
-  fmd3Node->SetLineColor(3);
-  fmd3Node->SetVisibility(1);
-
-  fmd3Node->cd();
-  TNode* innerNode = createRing("inner", innerRl, innerRh, innerTh, 
-                               siThick, waferR, staggering, innerZ-z);
-
-
-  fmd3Node->cd();
-  TNode* outerNode = createRing("outer", outerRl, outerRh, outerTh, 
-                               siThick, waferR, staggering, outerZ-z);
-  
-
-  fmd3Node->cd();
-  TNode* supportNode = createSupport(noseRl, noseRh, noseDz, noseZ-z, 
-                                    backRl, backRh, backDz, coneL,
-                                    beamW, beamDz, flangeR);
-  
-  TCanvas* c = new TCanvas("c", "c", 800, 800);
-  c->SetFillColor(1);
-  geometry->Draw();
-  // c->x3d("ogl");
-}
-//____________________________________________________________________
-//
-// EOF
-//
diff --git a/FMD/scripts/TestIndex.C b/FMD/scripts/TestIndex.C
new file mode 100644 (file)
index 0000000..20c2140
--- /dev/null
@@ -0,0 +1,246 @@
+//____________________________________________________________________
+//
+// $Id$
+//
+/** @file    TestIndex.C
+    @author  Christian Holm Christensen <cholm@nbi.dk>
+    @date    Thu Mar 30 01:20:02 2006
+    @brief   Test AliFMDIndex and AliFMDObjIndex
+*/
+#ifndef __CINT__
+#include <AliFMDIndex.h>
+#include <iostream>
+#include <TFile.h>
+#endif
+/** @defgroup index_test Test of AliFMDIndex and AliFMDObjIndex 
+    @ingroup FMD_script
+*/
+/** Write an AliFMDIndex object to output stream
+    @ingroup index_test
+    @param o Stream
+    @param i Index object 
+    @return @a o */
+std::ostream&
+operator << (std::ostream& o, const AliFMDIndex& i) 
+{
+  return o << i.Name();
+}
+
+/** Do the comparision, and print to standard out
+    @ingroup index_test
+    @param lhs Left hand side
+    @param rhs Right hand side
+    @return @f$ lhs < rhs @f$ */
+bool 
+cmp(const AliFMDIndex& lhs, const AliFMDIndex& rhs)
+{
+  bool ret = lhs < rhs;
+  std::cout << "    (" << lhs << " <  " << rhs << "): " << ret << std::endl;
+  return ret;
+}
+
+/** Test if two index objects are equivilant (are the same).
+    Equivilance is defined as 
+    @f[
+    lhs \equiv rhs: \not (lhs < rhs \vee rhs < lhs)
+    @f]
+    @ingroup index_test
+    @param lhs Left hand side 
+    @param rhs Right hand side 
+    @return @c true if @a lhs and @a rhs are equivilant */
+bool 
+equiv(const AliFMDIndex& lhs, const AliFMDIndex& rhs)
+{
+  bool ret = !(cmp(lhs,rhs) || cmp(rhs,lhs));
+  std::cout << "    (" << lhs << " <> " << rhs << "): " << ret << std::endl;
+  return ret;
+}
+
+/** Check that @f$ \not (x < x)@f$
+    @ingroup index_test
+    @param x Object to test
+    @return @c true if @a x is not less than itself */
+bool
+self(const AliFMDIndex& x)
+{
+  bool ret = !cmp(x,x);
+  std::cout << "  !(" << x << " < " << x << "): " << ret << std::endl;
+  return ret;
+}
+
+/** Check if @f$ a \wedge b \Rightarrow c@f$
+    @ingroup index_test
+    @param a First condition
+    @param b Second condition
+    @param c Implication
+    @return @c true if the implication is valid */
+bool
+imply(bool a, bool b, bool c) 
+{
+  bool ret = ((a && b) && c) || (!(a && b));
+  return ret;
+}
+
+/** Check if the comparison operator is transitive, that is 
+    @f[
+    (x < y \wedge y < z) \Rightarrow x < z
+    @f]
+    @ingroup index_test
+    @param x First object
+    @param y Second object
+    @param z Third object
+    @return @c true if the implication is met. */
+bool
+trans(const AliFMDIndex& x, const AliFMDIndex& y, const AliFMDIndex& z)
+{
+  bool ret  = imply(cmp(x,y), cmp(y,z), cmp(x,z));
+  std::cout << "  (" << x << " < " << y << " && " << y << " < " << z 
+           << ") => " << x << " < " << z << " " 
+           << (ret ? "holds" : "violated") << std::endl;
+  return ret;
+  
+}
+
+/** Check that the comparison operator preserves equivilance, that is 
+    @f[ 
+    (x \equiv y \wedge y \equiv z) \Rightarrow (x \equiv z)
+    @f]
+    @ingroup index_test
+    @param x First object
+    @param y Second argument
+    @param z Third object
+    @return @c true if the implication holds. */
+bool
+equiv(const AliFMDIndex& x, const AliFMDIndex& y, const AliFMDIndex& z)
+{
+  bool ret  = imply(equiv(x,y), equiv(y,z), equiv(x,z));
+  std::cout << "  (" << x << " <> " << y << " && " << y << " <> " << z 
+           << ") => " << x << " <> " << z << " " 
+           << (ret ? "holds" : "violated") << std::endl;
+  return ret;
+}
+
+
+/** Check if the comparison operator is a @e strictly @e weak @e ordering
+    @ingroup index_test
+ */
+void
+TestIndex() 
+{
+  AliFMDIndex i1(1, 'I', 5, 63);
+  AliFMDIndex i2(i1);  
+  AliFMDIndex i3(i1);  
+  AliFMDIndex i4(1, 'I', 5, 127);
+  AliFMDIndex i5(1, 'I', 15, 63);
+  AliFMDIndex i6(2, 'O', 15, 60);
+  std::cout << "Is !(x < x): " << std::endl;
+  self(i1);
+  std::cout << "Does (x < y && y < z) imply x < z: " << std::endl;
+  // true, true
+  trans(i1, i4, i5);
+  // true, false
+  trans(i1, i6, i5);
+  // false, true
+  trans(i4, i1, i5);
+  // false, false
+  trans(i6, i5, i1);
+  std::cout << "Does !(x < y || x > y) && !(y < z || z > y) imply "
+           << "!(x < z || x > z)" << std::endl;
+  // true, true
+  equiv(i1, i2, i3);
+  // true, false
+  equiv(i1, i2, i4);
+  // false, true
+  equiv(i4, i1, i2);
+  // false, false
+  equiv(i1, i4, i5);
+  
+
+
+  TFile* file = TFile::Open("index.root", "RECREATE");
+  file->WriteObject(&i1,"i1");
+  file->Write();
+  file->Close();
+  
+  AliFMDIndex* i7 = 0;
+  file = TFile::Open("index.root", "READ");
+  file->GetObject("i1", i7);
+  file->Close();
+  std::cout << *i7 << " == " << i1 << ": " << (*i7 == i1) << std::endl;
+  
+}
+
+/** Check if the comparison operator is a @e strictly @e weak @e ordering
+    @ingroup index_test
+ */
+void
+TestObjIndex() 
+{
+  AliFMDObjIndex i1(1, 'I', 5, 63);
+  AliFMDObjIndex i2(i1);  
+  AliFMDObjIndex i3(i1);  
+  AliFMDObjIndex i4(1, 'I', 5, 127);
+  AliFMDObjIndex i5(1, 'I', 15, 63);
+  AliFMDObjIndex i6(2, 'O', 15, 60);
+  std::cout << "Is !(x < x): " << std::endl;
+  self(i1);
+  std::cout << "Does (x < y && y < z) imply x < z: " << std::endl;
+  // true, true
+  trans(i1, i4, i5);
+  // true, false
+  trans(i1, i6, i5);
+  // false, true
+  trans(i4, i1, i5);
+  // false, false
+  trans(i6, i5, i1);
+  std::cout << "Does !(x < y || x > y) && !(y < z || z > y) imply "
+           << "!(x < z || x > z)" << std::endl;
+  // true, true
+  equiv(i1, i2, i3);
+  // true, false
+  equiv(i1, i2, i4);
+  // false, true
+  equiv(i4, i1, i2);
+  // false, false
+  equiv(i1, i4, i5);
+  
+
+
+  TFile* file = TFile::Open("index.root", "RECREATE");
+  i1.Write("i1");
+  file->Write();
+  file->Close();
+  
+  file = TFile::Open("index.root", "READ");
+  AliFMDObjIndex* i7 = (AliFMDObjIndex*)(file->Get("i1"));
+  file->Close();
+  std::cout << *i7 << " == " << i1 << ": " << (*i7 == i1) << std::endl;
+  
+}
+
+/** Check that we can sort an array of index objects
+    @ingroup index_test
+ */
+void
+SortIndex()
+{
+  TList l;
+  for (size_t i = 0; i < 30; i++) {
+    UShort_t det  = gRandom->Integer(3)+1;
+    Char_t   ring = (gRandom->Uniform() > .5 ? 'O' : 'I');
+    UShort_t sec  = gRandom->Integer(ring == 'I' ?  20 :  40);
+    UShort_t str  = gRandom->Integer(ring == 'I' ? 512 : 256);
+    l.AddAt(new AliFMDObjIndex(det, ring, sec, str), i);
+  }
+  std::cout << "Before sort" << std::endl;
+  l.ls();
+  l.Sort();
+  std::cout << "After sort" << std::endl;
+  l.ls();
+}
+
+    
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/scripts/WriteMedArrays.C b/FMD/scripts/WriteMedArrays.C
new file mode 100644 (file)
index 0000000..e639fe6
--- /dev/null
@@ -0,0 +1,30 @@
+void
+WriteMedArrays()
+{
+  TFile* file = TFile::Open("medid.root", "RECREATE");
+  if (!file) {
+    Warning("WriteMedArrays", "failed to open medid.root");
+    return;
+  }
+  TObjArray* modules = gAlice->Modules();
+  if (!modules) {
+    Warning("WriteMedArrays", "failed to get modules");
+    return;
+  }
+  TIter next(modules);
+  AliModule* module = 0;
+  while ((module = static_cast<AliModule*>(next()))) {
+    Info("WriteMedArrays", "Getting medium id's for %s", module->GetName());
+    TArrayI* mediumIds = module->GetIdtmed();
+    if (!mediumIds) {
+      Warning("WriteMedArrays", "No medium id's for %s", module->GetName());
+      continue;
+    }
+    file->WriteObject(mediumIds,module->GetName());
+  }
+  file->Write();
+  file->Close();
+}
+
+  
+