AliFMD&
AliFMD::operator=(const AliFMD& other)
{
+ // Assignment operator
AliDetector::operator=(other);
fSDigits = other.fSDigits;
fNsdigits = other.fNsdigits;
Float_t pz,
Float_t edep,
Int_t pdg,
- Float_t t)
+ Float_t t,
+ Float_t l,
+ Bool_t stop)
{
//
// Add a hit to the list
// edep Energy deposited by track
// pdg Track's particle Id #
// t Time when the track hit
+ // l Track length through the material.
+ // stop Whether track was stopped or disappeared
//
TClonesArray& a = *(HitsArray());
// Search through the list of already registered hits, and see if we
}
// If hit wasn't already registered, do so know.
hit = new (a[fNhits]) AliFMDHit(fIshunt, track, detector, ring, sector,
- strip, x, y, z, px, py, pz, edep, pdg, t);
+ strip, x, y, z, px, py, pz, edep, pdg, t,
+ l, stop);
fNhits++;
return hit;
}
//____________________________________________________________________
//
// Manager class for the FMD - Base class.
+// AliFMDv1, AliFMDv0, and AliFMDAlla
+// provides concrete implementations.
+// This class is sooooo crowded
//
#ifndef ALIDETECTOR_H
# include <AliDetector.h>
Float_t pz=0,
Float_t edep=0,
Int_t pdg=0,
- Float_t t=0);
+ Float_t t=0,
+ Float_t len=0,
+ Bool_t stopped=kFALSE);
virtual void AddDigit(Int_t *digits, Int_t* notused=0);
virtual void AddDigitByFields(UShort_t detector=0,
Char_t ring='\0',
// Concrete implementation of AliFMDDetector
//
// This implements the geometry for FMD1
+// FMD1 has only one ring, of type `inner'.
+// It is sitting at z=320.
+// It is the FMD ring with highest eta
//
#include "AliFMD1.h" // ALIFMD1_H
#include "AliFMDRing.h" // ALIFMDRING_H
class AliFMDRing;
//__________________________________________________________________
-/** Geometry description and parameters of the FMD1
- detector.
-
- The FMD1 only has one ring.
-*/
+// Geometry description and parameters of the FMD1
+// detector.
+//
+// The FMD1 only has one ring.
+//
class AliFMD1 : public AliFMDDetector
{
public:
// Concrete implementation of AliFMDDetector
//
// This implements the geometry for FMD2
+// The FMD2 has two ring, one of both types.
+// FMD2 is mounted on the space-frame via 4 flanges
+// Support is not fleshed ot yet.
//
#include "AliFMD2.h" // ALIFMD2_H
#include "AliFMDRing.h" // ALIFMDRING_H
# include "AliFMDDetector.h"
#endif
-/** Geometry description and parameters of the FMD2
- detector.
-*/
+// Geometry description and parameters of the FMD2
+// detector.
+// This has two rings.
class AliFMD2 : public AliFMDDetector
{
protected:
//
// Concrete implementation of AliFMDDetector
//
-// This implements the geometry for FMD3
+// This implements the geometry for FMD3.
+// This has 2 rings.
+// The support of the FMD3 is a carbon-fibre cone, attached to the ITS
+// support via flanges. The cone also supports the beam-pipe.
//
#include "AliFMD3.h" // ALIFMD3_H
#include "AliLog.h" // ALILOG_H
AliFMD3::AliFMD3(AliFMDRing* inner, AliFMDRing* outer)
: AliFMDDetector(3, inner, outer)
{
+ // Constructor.
SetInnerZ(-62.8);
SetOuterZ(-75.2);
SetNoseZ();
void
AliFMD3::Init()
{
+ // Initialize
AliFMDDetector::Init();
SetInnerHoneyHighR(GetOuterHoneyHighR());
Double_t zdist = fConeLength - fBackLength - fNoseLength;
{
// Calculate the cone radius at Z
if (fAlpha < 0) {
- Warning("ConeR", "alpha not set: %lf", fAlpha);
+ AliWarning(Form("alpha not set: %lf", fAlpha));
return -1;
}
if (z > fNoseZ) {
- Warning("ConeR", "z=%lf is before start of cone %lf", z, fNoseZ);
+ AliWarning(Form("z=%lf is before start of cone %lf", z, fNoseZ));
return -1;
}
if (z < fOuterZ - fOuter->GetRingDepth() - fHoneycombThickness) {
- Warning("ConeR", "z=%lf is after end of cone %lf", z,
- fOuterZ - fOuter->GetRingDepth() - fHoneycombThickness);
+ AliWarning(Form("z=%lf is after end of cone %lf", z,
+ fOuterZ - fOuter->GetRingDepth() - fHoneycombThickness));
return -1;
}
Double_t e = fBeamThickness / TMath::Cos(TMath::ATan(fAlpha));
# include "AliFMDDetector.h"
#endif
-/** Geometry description and parameters of the FMD3 detector.
-
- FMD3 has a fairly complicated support structure
-*/
+// Geometry description and parameters of the FMD3 detector.
+// FMD3 has a fairly complicated support structure.
+// The cone also supports the beam-pipe.
+//
class AliFMD3 : public AliFMDDetector
{
public:
/* $Id$ */
//__________________________________________________________
//
-// Map of per strip Bool_t information
-//
+// Map of Bool_t for each FMD strip
+// Used in calibration and the like classes.
+// Used amoung other things for dead-channel map
+// Can also be used for other stuff too
// Created Mon Nov 8 12:51:51 2004 by Christian Holm Christensen
//
#include "AliFMDBoolMap.h" //ALIFMDBOOLMAP_H
//__________________________________________________________
//
// Map of Bool_t for each FMD strip
+// Used in calibration and the like classes.
+// Used amoung other things for dead-channel map
//
#ifndef ALIFMDMAP_H
# include <AliFMDMap.h>
//
#include "AliFMDDetector.h" // ALIFMDSUBDETECTOR_H
#include "AliFMDRing.h" // ALIFMDRING_H
-#include <AliLog.h> // ALILOG_H
//====================================================================
ClassImp(AliFMDDetector)
fInner(inner),
fOuter(outer)
{
+ // Constructor
+ //
+ // ID Id of detector (1,2, or 3)
+ // INNER Inner ring geometry
+ // OUTER Outer ring geometry (if any)
+ //
SetHoneycombThickness();
SetAlThickness();
SetInnerHoneyLowR(0);
}
//____________________________________________________________________
+AliFMDDetector::AliFMDDetector(const AliFMDDetector& other)
+ : TNamed(other),
+ fId(other.fId),
+ fInner(other.fInner),
+ fOuter(other.fOuter)
+{
+ // Copy constructor
+ SetHoneycombThickness(other.GetHoneycombThickness());
+ SetAlThickness(other.GetAlThickness());
+ SetInnerHoneyLowR(other.GetInnerHoneyLowR());
+ SetInnerHoneyHighR(other.GetInnerHoneyHighR());
+ SetInnerZ(other.GetInnerZ());
+ SetOuterZ(other.GetOuterZ());
+ SetOuterHoneyLowR(other.GetOuterHoneyLowR());
+ SetOuterHoneyHighR(other.GetOuterHoneyHighR());
+}
+
+//____________________________________________________________________
+AliFMDDetector&
+AliFMDDetector::operator=(const AliFMDDetector& other)
+{
+ // Assignment operator
+ SetName(other.GetName());
+ SetTitle(other.GetTitle());
+ fId = other.fId;
+ fInner = other.fInner;
+ fOuter = other.fOuter;
+ SetHoneycombThickness(other.GetHoneycombThickness());
+ SetAlThickness(other.GetAlThickness());
+ SetInnerHoneyLowR(other.GetInnerHoneyLowR());
+ SetInnerHoneyHighR(other.GetInnerHoneyHighR());
+ SetInnerZ(other.GetInnerZ());
+ SetOuterZ(other.GetOuterZ());
+ SetOuterHoneyLowR(other.GetOuterHoneyLowR());
+ SetOuterHoneyHighR(other.GetOuterHoneyHighR());
+ return *this;
+}
+
+//____________________________________________________________________
void
AliFMDDetector::Init()
{
+ // Initialize.
if (fInner) {
SetInnerHoneyLowR(fInner->GetLowR() + 1.);
SetInnerHoneyHighR(fInner->GetHighR() + 1.);
AliFMDRing*
AliFMDDetector::GetRing(Char_t id) const
{
+ // Get the specified ring
+ //
+ // ID Id of ring ('I' or 'O')
+ //
switch (id) {
case 'i':
case 'I': return GetInner();
Double_t
AliFMDDetector::GetRingZ(Char_t id) const
{
+ // Get the z-coordinate specified ring
+ //
+ // ID Id of ring ('I' or 'O')
+ //
switch (id) {
case 'i':
case 'I': return GetInnerZ();
Double_t& y,
Double_t& z) const
{
+ // Translate detector coordinates (this,ring,sector,strip) into
+ // (x,y,z) coordinates (in global reference frame)
AliFMDRing* r = GetRing(ring);
if (!r) return;
z = GetRingZ(ring);
UShort_t& sector,
UShort_t& strip) const
{
+ // Translate (x,y,z) coordinates (in global reference frame) into
+ // detector coordinates (this,ring,sector,strip).
AliFMDRing* rng = 0;
ring = -1;
for (int j = 0; j < 2; j++) {
{
public:
AliFMDDetector(Int_t id, AliFMDRing* inner, AliFMDRing* outer);
+ AliFMDDetector(const AliFMDDetector& other);
+ AliFMDDetector& operator=(const AliFMDDetector& other);
virtual ~AliFMDDetector() {}
/** Initialize the geometry */
virtual void Init();
#define ALIFMDDIGIT_H
//___________________________________________________________________
//
-// Digits classes for the FMD
+// Digits classes for the FMD
+// AliFMDBaseDigit - base class
+// AliFMDDigit - Normal (smeared) digit
+// AliFMDSDigit - Summable (non-smeared) digit
//
#ifndef ROOT_TObject
# include <TObject.h>
//____________________________________________________________________
class AliFMDDigit : public AliFMDBaseDigit
{
-protected:
- UShort_t fCount1; // Digital signal
- Short_t fCount2; // Digital signal (-1 if not used)
- Short_t fCount3; // Digital signal (-1 if not used)
public:
AliFMDDigit();
AliFMDDigit(UShort_t detector,
Short_t Count3() const { return fCount3; }
UShort_t Counts() const;
void Print(Option_t* opt="") const;
+protected:
+ UShort_t fCount1; // Digital signal
+ Short_t fCount2; // Digital signal (-1 if not used)
+ Short_t fCount3; // Digital signal (-1 if not used)
ClassDef(AliFMDDigit,1) // Normal FMD digit
};
//____________________________________________________________________
class AliFMDSDigit : public AliFMDBaseDigit
{
-protected:
- Float_t fEdep; // Energy deposited
- UShort_t fCount1; // Digital signal
- Short_t fCount2; // Digital signal (-1 if not used)
- Short_t fCount3; // Digital signal (-1 if not used)
public:
AliFMDSDigit();
AliFMDSDigit(UShort_t detector,
Float_t Edep() const { return fEdep; }
UShort_t Counts() const;
void Print(Option_t* opt="") const;
+protected:
+ Float_t fEdep; // Energy deposited
+ UShort_t fCount1; // Digital signal
+ Short_t fCount2; // Digital signal (-1 if not used)
+ Short_t fCount3; // Digital signal (-1 if not used)
ClassDef(AliFMDSDigit,1) // Summable FMD digit
};
//
// = E + (l - E) * ext(-B * t)
//
- const Float_t mipEnergy = 1.664 * siThickness * siDensity;
- const Float_t convf = (1 / mipEnergy * Float_t(fAltroChannelSize)
- / fVA1MipRange);
- UShort_t ped = MakePedestal();
+ Float_t mipEnergy = 1.664 * siThickness * siDensity;
+ Float_t convF = (1/mipEnergy*Float_t(fAltroChannelSize)/fVA1MipRange);
+ UShort_t ped = MakePedestal();
// In case we don't oversample, just return the end value.
if (fSampleRate == 1) {
- counts[0] = UShort_t(TMath::Min(edep * convf + ped,
+ counts[0] = UShort_t(TMath::Min(edep * convF + ped,
Float_t(fAltroChannelSize)));
return;
}
-
+
// Create a pedestal
Int_t n = fSampleRate;
- Float_t B = fShapingTime;
+ Float_t b = fShapingTime;
for (Ssiz_t i = 0; i < n; i++) {
Float_t t = Float_t(i) / n;
- Float_t s = edep + (last - edep) * TMath::Exp(-B * t);
- counts[i] = UShort_t(TMath::Min(s * convf + ped,
+ Float_t s = edep + (last - edep) * TMath::Exp(-b * t);
+ counts[i] = UShort_t(TMath::Min(s * convF + ped,
Float_t(fAltroChannelSize)));
}
}
UShort_t
AliFMDDigitizer::MakePedestal() const
{
+ // Make a pedestal
return UShort_t(TMath::Max(gRandom->Gaus(fPedestal, fPedestalWidth), 0.));
}
Short_t count2,
Short_t count3) const
{
+ // Add a digit
fmd->AddDigitByFields(detector, ring, sector, strip, count1, count2, count3);
}
UShort_t nhits,
const TArrayI& counts)
{
+ // Check that digit is consistent
Int_t integral = counts[0];
if (counts[1] >= 0) integral += counts[1];
if (counts[2] >= 0) integral += counts[2];
integral -= Int_t(fPedestal + 2 * fPedestalWidth);
if (integral < 0) integral = 0;
- Float_t convf = Float_t(fVA1MipRange) / fAltroChannelSize;
- Float_t mips = integral * convf;
+ Float_t convF = Float_t(fVA1MipRange) / fAltroChannelSize;
+ Float_t mips = integral * convF;
if (mips > Float_t(nhits) + .5 || mips < Float_t(nhits) - .5)
Warning("CheckDigit", "Digit -> %4.2f MIPS != %d +/- .5 hits",
mips, nhits);
//____________________________________________________________________
AliFMDSDigitizer::~AliFMDSDigitizer()
{
+ // Destructor
AliLoader* loader = fRunLoader->GetLoader("FMDLoader");
loader->CleanSDigitizer();
}
Short_t count2,
Short_t count3) const
{
+ // Add a summable digit
fmd->AddSDigitByFields(detector, ring, sector, strip, edep,
count1, count2, count3);
}
UShort_t fSampleRate; // Times the ALTRO samples pre-amp.
Float_t fShapingTime; // Shaping profile parameter
+ AliFMDBaseDigitizer(const AliFMDBaseDigitizer& o)
+ : AliDigitizer(o) {}
+ AliFMDBaseDigitizer& operator=(const AliFMDBaseDigitizer&) { return *this; }
ClassDef(AliFMDBaseDigitizer,0) // Base class for FMD digitizers
};
other.fMaxStrips),
fData(0)
{
+ // Copy constructor
fData = new AliFMDEdepHitPair[fMaxDetectors * fMaxRings *
fMaxSectors * fMaxStrips];
for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
AliFMDEdepMap&
AliFMDEdepMap::operator=(const AliFMDEdepMap& other)
{
+ // Assignment operator
fMaxDetectors = other.fMaxDetectors;
fMaxRings = other.fMaxRings;
fMaxSectors = other.fMaxSectors;
//____________________________________________________________________
void
+AliFMDEdepMap::Reset()
+{
+ // Reset to zero
+ for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
+ i++) { fData[i].fEdep = 0; fData[i].fN = 0; };
+}
+
+//____________________________________________________________________
+void
AliFMDEdepMap::Reset(const AliFMDEdepHitPair& val)
{
+ // Reset to val
for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
- i++) fData[i] = val;
+ i++) { fData[i].fEdep = val.fEdep; fData[i].fN = val.fN; };
}
//____________________________________________________________________
#endif
//____________________________________________________________________
//
-// Cache of Energy deposited, hit information perr strip
+// Cache of Energy deposited, hit information per strip.
+// Contains a pair of fEdep and fN
+// fEdep is the summed energy deposition, and fN is the number of hits
//
//____________________________________________________________________
class AliFMDEdepHitPair
{
public:
- Float_t fEdep;
- UShort_t fN;
+ Float_t fEdep; // summed energy deposition
+ UShort_t fN; // Number of hits
AliFMDEdepHitPair() : fEdep(0), fN(0) {}
+ AliFMDEdepHitPair& operator=(const AliFMDEdepHitPair& o)
+ { fEdep = o.fEdep; fN = o.fN; return *this; }
+ AliFMDEdepHitPair(const AliFMDEdepHitPair& o) : fEdep(o.fEdep), fN(o.fN) {}
};
//____________________________________________________________________
size_t maxStr = kMaxStrips);
virtual ~AliFMDEdepMap() { delete [] fData; }
AliFMDEdepMap& operator=(const AliFMDEdepMap& other);
- virtual void Reset(const AliFMDEdepHitPair& val=AliFMDEdepHitPair());
+ virtual void Reset();
+ virtual void Reset(const AliFMDEdepHitPair& val);
virtual AliFMDEdepHitPair& operator()(UShort_t detector,
Char_t ring,
UShort_t sector,
AliFMDGeometry::SetActive(Int_t* active, Int_t n)
{
fActive.Set(n);
- for (Int_t i = 0; i < n; i++) fActive[i] = active[i];
+ for (Int_t i = 0; i < n; i++) {
+ AliDebug(1, Form("Active vol id # %d: %d", i, active[i]));
+ fActive[i] = active[i];
+ }
}
//____________________________________________________________________
#ifndef ALIFMDGEOMETRY_H
#define ALIFMDGEOMETRY_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights
+ * reserved.
+ *
+ * Latest changes by Christian Holm Christensen <cholm@nbi.dk>
+ *
+ * See cxx source for full Copyright notice
+ */
//____________________________________________________________________
//
-// $Id$
-//
// Forward Multiplicity Detector based on Silicon wafers.
//
// This class is a singleton that handles the geometry parameters of
fPlastic(0)
{
// Default constructor
- fActiveId.Set(4);
+ fActiveId.Set(2);
}
//____________________________________________________________________
// fmd Pointer to AliFMD object
// detailed Whether to make a detailed simulation or not
//
- fActiveId.Set(4);
+ fActiveId.Set(2);
}
case 'i':
case 'I': fActiveId[0] = sid; break;
case 'o':
- case 'O': fActiveId[2] = sid; break;
+ case 'O': fActiveId[1] = sid; break;
}
// Shape of Printed circuit Board
case 'i':
case 'I': fActiveId[0] = sid; break;
case 'o':
- case 'O': fActiveId[2] = sid; break;
+ case 'O': fActiveId[1] = sid; break;
}
// Shape of Printed circuit Board
//____________________________________________________________________
//
// Hits in the FMD
-//
-// Latest changes by Christian Holm Christensen
+// Contains information on:
+// Position of hit
+// Momentum of track
+// PID of track
+// Energy loss of track
+// Track #
+// Track path length
+// Track stopping status.
+// Latest changes by Christian Holm Christensen
//
#include "AliFMDHit.h" // ALIFMDHIT_H
#include "AliLog.h" // ALILOG_H
fPz(0),
fPdg(0),
fEdep(0),
- fTime(0)
+ fTime(0),
+ fLength(0),
+ fStop(0)
{
+ // Default CTOR
fX = fY = fZ = 0;
}
Float_t pz,
Float_t edep,
Int_t pdg,
- Float_t t)
+ Float_t t,
+ Float_t l,
+ Bool_t stop)
: AliHit(shunt, track),
fDetector(detector),
fRing(ring),
fPz(pz),
fPdg(pdg),
fEdep(edep),
- fTime(t)
+ fTime(t),
+ fLength(l),
+ fStop(stop)
{
// Normal FMD hit ctor
//
cout << "\tPDG:\t" << fPdg << " " << (pdg ? pdg->GetName() : "?") << "\n"
<< "\tP:\t(" << fPx << "," << fPy << "," << fPz << ") "<<P() << "\n"
<< "\tX:\t" << fX << "," << fY << "," << fZ << "\n"
- << "\tTrack #:\t" << fTrack << std::endl;
+ << "\tTrack #:\t" << fTrack << "\tLength:\t"
+ << fLength << "cm\t" << (IsStop() ? "stopped" : "") << std::endl;
}
}
Float_t pz=0,
Float_t edep=0,
Int_t pdg=0,
- Float_t t=0);
+ Float_t t=0,
+ Float_t l=0,
+ Bool_t stop=kFALSE);
virtual ~AliFMDHit() {}
UShort_t Detector() const { return fDetector; }
Float_t Q() const;
Int_t Pdg() const { return fPdg; }
Float_t Time() const { return fTime; }
+ Float_t Length() const { return fLength; }
+ Bool_t IsStop() const { return fStop; }
void Print(Option_t* opt="") const;
void SetEdep(Float_t edep) { fEdep = edep; }
Int_t fPdg; // Particles PDG code
Float_t fEdep; // Energy deposition
Float_t fTime; // Particle's time of flight
-
- ClassDef(AliFMDHit,1) //Hits for detector FMD
+ Float_t fLength; // Track length through material.
+ Bool_t fStop; // Whether track was stopped.
+
+ ClassDef(AliFMDHit,2) //Hits for detector FMD
};
#endif
//____________________________________________________________________
//____________________________________________________________________
//
-//
+// Base class for caches of per-strip information.
+// This is used to index a strip.
+// Data stored depends on derived class.
+// This class provides some common infra-structure.
+// Derived classes sould define Reset, and operator().
//
#include "AliFMDMap.h" // ALIFMDMAP_H
+#include "AliLog.h"
//____________________________________________________________________
ClassImp(AliFMDMap)
Int_t idx = CheckIndex(det, ring, sec, str);
if (idx < 0) {
size_t ringi = (ring == 'I' || ring == 'i' ? 0 : 1);
- Fatal("CalcIndex", "Index (%d,'%c',%d,%d) out of bounds, "
- "in particular the %s index ",
- det, ring, sec, str,
- (det >= fMaxDetectors ? "Detector" :
- (ringi >= fMaxRings ? "Ring" :
- (sec >= fMaxSectors ? "Sector" : "Strip"))));
+ AliFatal(Form("CalcIndex", "Index (%d,'%c',%d,%d) out of bounds, "
+ "in particular the %s index ",
+ det, ring, sec, str,
+ (det >= fMaxDetectors ? "Detector" :
+ (ringi >= fMaxRings ? "Ring" :
+ (sec >= fMaxSectors ? "Sector" : "Strip")))));
return 0;
}
return size_t(idx);
//____________________________________________________________________
//
// Base class for caches of per-strip information.
+// This is used to index a strip.
+// Data stored depends on derived class.
//
-
class AliFMDMap : public TObject
{
public:
//____________________________________________________________________
//
// Base class for reconstructed charged particle multiplicty in the
-// FMD. The class contains the field fMethod which is a flag set by
+// FMD.
+// The class contains the field fMethod which is a flag set by
// the AliFMDMultAlgorithm that created the object. The flag tells us
// which algorithm was used to create the data stored in the object.
//
#ifndef ALIFMDMULT_H
#define ALIFMDMULT_H
-/* Reconstracted Particles Class: has number of reconstructed
- * particles in sectors from NumOfMinSector to NumberOfMaxSector()
- * rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
- */
+// Reconstracted Particles Class: has number of reconstructed
+// particles in sectors from NumOfMinSector to NumberOfMaxSector()
+// rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
+//
#ifndef ROOT_TObject
# include <TObject.h>
#endif
fTreeR(0),
fMult(0),
fFMD(0)
-{}
+{
+ // Default CTOR
+}
//____________________________________________________________________
AliFMDMultAlgorithm::~AliFMDMultAlgorithm()
{
+ // DTOR
if (fMult) {
fMult->Delete();
delete fMult;
void
AliFMDMultAlgorithm::PreEvent(TTree* treeR, Float_t /* ipZ */)
{
+ // Executed before each event.
if (fMult) fMult->Clear();
fNMult = 0;
fTreeR = treeR;
* See cxx source for full Copyright notice
*/
/* $Id$ */
+
//____________________________________________________________________
//
-// Base class for FMD reconstruction algorithms.
-//
+// Base class for FMD reconstruction algorithms.
// Derived classes will implement various ways of reconstructing the
// charge particle multiplicity in the FMD.
//
// [See also the AliFMDReconstructor class]
//
// This class reconstructs the multiplicity based on the assumption
-// that all particles are minimum ionizing.
+// that all particles are minimum ionizing.
+// Hence, the name `naiive'
//
#include "AliFMD.h" // ALIFMD_H
#include "AliFMDMultNaiive.h" // ALIFMDMULTNAIIVE_H
//
// This class reconstructs the muliplicity in regions based on the
// ratio of empty to full strips.
+// Hence the name `poisson'
//
#include "AliFMD.h" // ALIFMD_H
#include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
fDeltaPhi(0),
fThreshold(0)
{
+ // CTOR
SetDeltaEta();
SetDeltaPhi();
SetThreshold();
#ifndef ALIFMDMULTREGION_H
#define ALIFMDMULTREGION_H
-/* Reconstracted Particles Class: has number of reconstructed
- * particles in sectors from NumOfMinSector to NumberOfMaxSector()
- * rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
- */
+// Reconstracted Particles Class: has number of reconstructed
+// particles in sectors from NumOfMinSector to NumberOfMaxSector()
+// rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
+//
#ifndef ALIFMDMULT_H
# include "AliFMDMult.h"
#endif
#ifndef ALIFMDMULTSTRIP_H
#define ALIFMDMULTSTRIP_H
-/* Reconstracted Particles Class: has number of reconstructed
- * particles in sectors from NumOfMinSector to NumberOfMaxSector()
- * rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
- */
+// Reconstracted Particles Class: has number of reconstructed
+// particles in sectors from NumOfMinSector to NumberOfMaxSector()
+// rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
+//
#ifndef ALIFMDMULT_H
# include "AliFMDMult.h"
#endif
//
// This class is a singleton that handles various parameters of
// the FMD detectors.
+// Eventually, this class will use the Conditions DB to get the
+// various parameters, which code can then request from here.
//
#include "AliFMDParameters.h" // ALIFMDPARAMETERS_H
#include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
//
// Singleton class to handle various parameters (not geometry) of the
// FMD
+// Should get ata fromm Conditions DB.
//
#ifndef ROOT_TNamed
# include <TNamed.h>
protected:
AliFMDParameters();
virtual ~AliFMDParameters() {}
- static AliFMDParameters* fgInstance;
+ static AliFMDParameters* fgInstance; // Static singleton instance
const Float_t fSiDeDxMip; // MIP dE/dx in Silicon
UShort_t fVA1MipRange; // # MIPs the pre-amp can do
//____________________________________________________________________
//
// Class to read ADC values from a AliRawReader object.
+// Note, that it uses an ALTRO reader, which is wrong.
+// Perhaps we need to implement it our selves
//
#ifndef ROOT_TTask
# include <TTask.h>
{
public:
AliFMDRawReader(AliRawReader* reader, TTree* array);
-
+ virtual ~AliFMDRawReader() {}
virtual void Exec(Option_t* option="");
protected:
TTree* fTree; //! Pointer to tree to read into
{
public:
AliFMDRawStream(AliRawReader* reader, UShort_t sampleRate=0);
+ virtual ~AliFMDRawStream() {}
Short_t Sector() const { return fRow; }
Char_t Ring() const { return (fSector == 0 ? 'I' : 'O'); }
//____________________________________________________________________
//
// Class to writer ADC values to a Raw File
+// Uses general ALTRO class - which is wrong
+// Should make it right!
//
#ifndef ROOT_TTask
# include <TTask.h>
{
public:
AliFMDRawWriter(AliFMD* fmd);
+ virtual ~AliFMDRawWriter() {}
virtual void Exec(Option_t* option="");
protected:
fId(id),
fVerticies(0)
{
+ // CTOR
SetBondingWidth();
SetWaferRadius();
SetSiThickness();
void
AliFMDRing::Init()
{
+ // Initialize
Double_t tanTheta = TMath::Tan(fTheta * TMath::Pi() / 180.);
Double_t tanTheta2 = TMath::Power(tanTheta,2);
Double_t r2 = TMath::Power(fWaferRadius,2);
TVector2*
AliFMDRing::GetVertex(Int_t i) const
{
+ // Get the i'th vertex of polygon shape
return static_cast<TVector2*>(fVerticies.At(i));
}
Double_t& y,
Double_t& z) const
{
+ // Translate detector coordinates (this,sector,strip) to global
+ // coordinates (x,y,z)
if (sector >= GetNSectors()) {
Error("Detector2XYZ", "Invalid sector number %d (>=%d) in ring %c",
sector, GetNSectors(), fId);
UShort_t& sector,
UShort_t& strip) const
{
+ // Translate global coordinates (x,y,z) to detector coordinates
+ // (this,sector,strip)
sector = strip = 0;
Double_t r = TMath::Sqrt(x * x + y * y);
Int_t str = Int_t((r - fMinR) / GetPitch());
//
// A map of per strip UShort_t information (for example ADC values,
// number of hits and so on).
+// Used for various calib information, and the like,
+// as well as in reconstruction.
+// Can be used elsewhere too.
//
#include "AliFMDUShortMap.h" // ALIFMDUSHORTMAP_H
other.fMaxStrips),
fData(0)
{
+ // CTOR
fData = new UShort_t[fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips];
for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
i++) fData[i] = other.fData[i];
AliFMDUShortMap&
AliFMDUShortMap::operator=(const AliFMDUShortMap& other)
{
+ // Assignment operator
fMaxDetectors = other.fMaxDetectors;
fMaxRings = other.fMaxRings;
fMaxSectors = other.fMaxSectors;
void
AliFMDUShortMap::Reset(const UShort_t& val)
{
+ // Reset to val
for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
i++) fData[i] = val;
}
detector = det;
//Double_t rz = fmd->GetDetector(detector)->GetRingZ(ring);
- Int_t n = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors();
+ AliFMDDetector* gdet = fmd->GetDetector(detector);
+ AliFMDRing* gring = gdet->GetRing(ring);
+ if (!gring)
+ AliFatal(Form("Ring %c not found (%s)", ring, mc->CurrentVolPath()));
+ Int_t n = gring->GetNSectors();
#if 0
if (rz < 0) {
Int_t s = ((n - sector + n / 2) % n) + 1;
}
#endif
if (sector < 1 || sector > n) {
- Warning("Step", "sector # %d out of range (0-%d)", sector-1, n-1);
+ AliWarning(Form("Step", "sector # %d out of range (0-%d)", sector-1, n-1));
return kFALSE;
}
sector--;
if (mc->IsTrackOut()) what.Append("out ");
Int_t mother = gAlice->GetMCApp()->GetPrimary(trackno);
- Warning("Step", "Track # %5d deposits a lot of energy\n"
- " Volume: %s\n"
- " Momentum: (%7.4f,%7.4f,%7.4f)\n"
- " PDG: %d (%s)\n"
- " Edep: %-14.7f keV (mother %d)\n"
- " p/m: %-7.4f/%-7.4f = %-14.7f\n"
- " Processes: %s\n"
- " What: %s\n",
- trackno, mc->CurrentVolPath(), p.X(), p.Y(), p.Z(),
- pdg, pname.Data(), edep, mother, p.P(), mass,
- poverm, processes.Data(), what.Data());
+ AliWarning(Form("Step", "Track # %5d deposits a lot of energy\n"
+ " Volume: %s\n"
+ " Momentum: (%7.4f,%7.4f,%7.4f)\n"
+ " PDG: %d (%s)\n"
+ " Edep: %-14.7f keV (mother %d)\n"
+ " p/m: %-7.4f/%-7.4f = %-14.7f\n"
+ " Processes: %s\n"
+ " What: %s\n",
+ trackno, mc->CurrentVolPath(), p.X(), p.Y(), p.Z(),
+ pdg, pname.Data(), edep, mother, p.P(), mass,
+ poverm, processes.Data(), what.Data()));
}
// Check that the track is actually within the active area
"Accumulated Edep=%g (%f,%f,%f)", trackno,
mc->CurrentVolPath(), edep, fCurrentDeltaE,
v.X(), v.Y(), v.Z()));
+ TVector3 cur(v.Vect());
+ cur -= fCurrentV.Vect();
+ Double_t len = cur.Mag();
AliFMDHit* h =
AddHitByFields(trackno, detector, ring, sector, strip,
fCurrentV.X(), fCurrentV.Y(), fCurrentV.Z(),
fCurrentP.X(), fCurrentP.Y(), fCurrentP.Z(),
- fCurrentDeltaE, fCurrentPdg, fCurrentV.T());
+ fCurrentDeltaE, fCurrentPdg, fCurrentV.T(),
+ len, mc->IsTrackDisappeared()||mc->IsTrackStop());
// Add a copy
if (isBad && fBad) {
new ((*fBad)[fBad->GetEntries()]) AliFMDHit(*h);
//____________________________________________________________________
//
// Manager class for the FMD - Detailed version.
-//
+// Implements the full geometry,
+// And does stepping
+//
#ifndef ALIFMD_H
# include "AliFMD.h"
#endif
kMuonCocktailCent1Single, //
kMuonCocktailPer1Single, //
kMuonCocktailPer4Single,
+ kFlatFMD1,
+ kFlatFMD2,
+ kFlatFMD3,
kEgMax
};
"kMuonCocktailPer4HighPt", //
"kMuonCocktailCent1Single", //
"kMuonCocktailPer1Single", //
- "kMuonCocktailPer4Single"
+ "kMuonCocktailPer4Single",
+ "kFMD1Flat",
+ "kFMD2Flat",
+ "kFMD3Flat"
};
//____________________________________________________________________
((TFluka*)gMC)->SetGeneratePemf(kFALSE);
else
((TFluka*)gMC)->SetGeneratePemf(kTRUE);
- TString neuxsc(gSystem->Which(".", "neuxsc.bin"));
- if (neusxc->IsNull()) {
- gSystem->Link("$(FLUPRO)/neuxsc.bin", "neuxsc.bin");
+ TString flupro(gSystem->Getenv("FLUPRO"));
+ if (flupro.IsNull())
+ Fatal("Config.C", "Environment variable FLUPRO not set");
+#if 0
+ char* files[] = { "brems_fin.bin", "cohff.bin", "elasct.bin",
+ "gxsect.bin", "nuclear.bin", "sigmapi.bin",
+ 0 };
+ char* file = files[0];
+ while (file) {
+ TString which(gSystem->Which(".", file));
+ if (which.IsNull()) {
+ if (gSystem->Symlink(Form("%s/%s", flupro.Data(), file), file)!=0)
+ Fatal("Config.C", "Couldn't link $(FLUPRO)/%s -> .", file);
+ }
+ file++;
}
+#endif
+ TString neuxsc(gSystem->Which(".", "neuxsc.bin"));
+ if (neuxsc.IsNull())
+ gSystem->Symlink(Form("%s/neuxsc_72.bin", flupro.Data()),
+ "neuxsc.bin");
gSystem->CopyFile("$(FLUPRO)/random.dat", "old.seed", kTRUE);
}
break;
gGener=gener;
}
break;
+ case kFlatFMD1:
+ {
+ comment = comment.Append(" Flat in FMD1 range");
+ AliGenBox* gener = AliGenBox(2000);
+ gener->SetPart(211);
+ gener->SetMomentumRange(3,4);
+ gener->SetPhiRange(0, 360);
+ gener->SetThetaRange(0.77, 3.08);
+ gGener = gener;
+ }
+ break;
+ case kFlatFMD2:
+ {
+ comment = comment.Append(" Flat in FMD2 range");
+ AliGenBox* gener = AliGenBox(2000);
+ gener->SetPart(211);
+ gener->SetMomentumRange(3,4);
+ gener->SetPhiRange(0, 360);
+ gener->SetThetaRange(2.95, 20.42);
+ gGener = gener;
+ }
+ break;
+ case kFlatFMD3:
+ {
+ comment = comment.Append(" Flat in FMD3 range");
+ AliGenBox* gener = AliGenBox(2000);
+ gener->SetPart(211);
+ gener->SetMomentumRange(3,4);
+ gener->SetPhiRange(0, 360);
+ gener->SetThetaRange(155.97, 176.73);
+ gGener = gener;
+ }
+ break;
+
default: break;
}
return gGener;
{
private:
TH2D* fElossVsPMQ; // Histogram
+ Int_t fPdg;
public:
- DrawHits()
- {
- Double_t emax = 1e4;
- Double_t emin = 1e-5;
- Int_t n = 901;
- TArrayF tkine(n);
- Float_t dp = 1/TMath::Log10(emax/emin)/10;
- Float_t pmin = TMath::Log10(emin);
- tkine[0] = emin;
- for (Int_t i=1; i < tkine.fN; i++) {
- Float_t el = pmin + i * dp;
- tkine[i] = TMath::Power(10, el);
- }
- Double_t dmin = .1;
- Double_t dmax = 10;
- Int_t nd = 1000;
- TArrayF eloss(nd+1);
- eloss[0] = dmin;
- for (Int_t i = 1; i < eloss.fN; i++){
- eloss[i] = dmin + i * (dmax-dmin)/(eloss.fN-1);
+ TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
+ {
+ TArrayF bins(n+1);
+ Float_t dp = n / TMath::Log10(max / min);
+ Float_t pmin = TMath::Log10(min);
+ bins[0] = min;
+ for (Int_t i = 1; i < n+1; i++) {
+ Float_t p = pmin + i / dp;
+ bins[i] = TMath::Power(10, p);
}
- fElossVsPMQ = new TH2D("bad", "#Delta E vs. p/(mq^{2})",
+ return bins;
+ }
+ DrawHits(Int_t m=1000, Double_t emin=1, Double_t emax=1000,
+ Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3,
+ Int_t pdg=211)
+ : fPdg(pdg)
+ {
+ TArrayF tkine(MakeLogScale(n, tmin, tmax));
+ TArrayF eloss(MakeLogScale(m, emin, emax));
+ fElossVsPMQ = new TH2D("bad", "#Delta E/#Delta x vs. p/(mq^{2})",
tkine.fN-1, tkine.fArray, eloss.fN-1, eloss.fArray);
- fElossVsPMQ->SetXTitle("p/(mq^{2})");
- fElossVsPMQ->SetYTitle("#Delta E [MeV]");
+ fElossVsPMQ->SetXTitle(Form("p%s", (fPdg == 0 ? "/(mq^{2})" : " [GeV]")));
+ fElossVsPMQ->SetYTitle("#Delta E/#Delta x [MeV/cm]");
}
Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
{
std::cout << "No hit" << std::endl;
return kFALSE;
}
- if (TMath::Abs(hit->Pdg()) != 211) return kTRUE;
- // Float_t pmq = 0;
- // Float_t q = hit->Q() / 3.;
- // if (hit->M() != 0 && hit->Q() != 0)
- // pmq = hit->P() / (hit->M()*q*q);
- Float_t pmq = hit->P();
- fElossVsPMQ->Fill(pmq, hit->Edep());
+
+ if (hit->IsStop()) return kTRUE;
+ Float_t x = hit->P();
+ Float_t y = hit->Edep()/hit->Length();
+ if (fPdg != 0) {
+ if (TMath::Abs(hit->Pdg()) != fPdg) return kTRUE;
+ }
+ else {
+ Float_t q = hit->Q() / 3.;
+ if (hit->M() == 0 || q == 0) return kTRUE;
+ x /= hit->M();
+ y /= q * q;
+ }
+ fElossVsPMQ->Fill(x, y);
return kTRUE;
}
Bool_t Finish()
{
gStyle->SetPalette(1);
gStyle->SetOptTitle(0);
+ gStyle->SetCanvasColor(0);
+ gStyle->SetCanvasBorderSize(0);
+ gStyle->SetPadColor(0);
+ gStyle->SetPadBorderSize(0);
fElossVsPMQ->SetStats(kFALSE);
fElossVsPMQ->Draw("COLZ");
return kTRUE;
--- /dev/null
+//
+// $Id$
+//
+// Script that contains a class to draw eloss from hits, versus ADC
+// counts from digits, using the AliFMDInputHits class in the util library.
+//
+// It draws the energy loss versus the p/(mq^2). It can be overlayed
+// with the Bethe-Bloc curve to show how the simulation behaves
+// relative to the expected.
+//
+// Use the script `Compile.C' to compile this class using ACLic.
+//
+#include <TH2D.h>
+#include <AliFMDHit.h>
+#include <AliFMDDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDEdepMap.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+
+class DrawHitsDigits : public AliFMDInputHits
+{
+private:
+ TH2D* fElossVsAdc; // Histogram
+ AliFMDEdepMap fMap;
+public:
+ TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
+ {
+ TArrayF bins(n+1);
+ Float_t dp = n / TMath::Log10(max / min);
+ Float_t pmin = TMath::Log10(min);
+ bins[0] = min;
+ for (Int_t i = 1; i < n+1; i++) {
+ Float_t p = pmin + i / dp;
+ bins[i] = TMath::Power(10, p);
+ }
+ return bins;
+ }
+ DrawHitsDigits(Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
+ Int_t m=1100, Double_t amin=-0.5, Double_t amax=1099.5)
+ {
+ AddLoad(kDigits);
+ TArrayF eloss(MakeLogScale(n, emin, emax));
+ TArrayF adcs(m+1);
+ adcs[0] = amin;
+ for (Int_t i = 1; i < m+1; i++) adcs[i] = adcs[i-1] + (amax-amin)/m;
+ fElossVsAdc = new TH2D("bad", "#Delta E vs. ADC",
+ eloss.fN-1, eloss.fArray, adcs.fN-1, adcs.fArray);
+ fElossVsAdc->SetXTitle("#Delta E/#Delta x [MeV/cm]");
+ fElossVsAdc->SetYTitle("ADC value");
+ }
+ Bool_t Begin(Int_t ev)
+ {
+ fMap.Reset();
+ return AliFMDInputHits::Begin(ev);
+ }
+ Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
+ {
+ // Cache the energy loss
+ if (!hit) return kFALSE;
+ UShort_t det = hit->Detector();
+ Char_t rng = hit->Ring();
+ UShort_t sec = hit->Sector();
+ UShort_t str = hit->Strip();
+ fMap(det, rng, sec, str).fEdep += hit->Edep();
+ fMap(det, rng, sec, str).fN++;
+ return kTRUE;
+ }
+ Bool_t Event()
+ {
+ if (!AliFMDInputHits::Event()) return kFALSE;
+ Int_t nEv = fTreeD->GetEntries();
+ for (Int_t i = 0; i < nEv; i++) {
+ Int_t digitRead = fTreeD->GetEntry(i);
+ if (digitRead <= 0) continue;
+ Int_t nDigit = fArrayD->GetEntries();
+ if (nDigit <= 0) continue;
+ for (Int_t j = 0; j < nDigit; j++) {
+ AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
+ if (!digit) continue;
+ UShort_t det = digit->Detector();
+ Char_t rng = digit->Ring();
+ UShort_t sec = digit->Sector();
+ UShort_t str = digit->Strip();
+ fElossVsAdc->Fill(fMap(det, rng, sec, str).fEdep, digit->Counts());
+ }
+ }
+ return kTRUE;
+ }
+ Bool_t Finish()
+ {
+ gStyle->SetPalette(1);
+ gStyle->SetOptTitle(0);
+ gStyle->SetCanvasColor(0);
+ gStyle->SetCanvasBorderSize(0);
+ gStyle->SetPadColor(0);
+ gStyle->SetPadBorderSize(0);
+ fElossVsAdc->SetStats(kFALSE);
+ fElossVsAdc->Draw("COLZ");
+ return kTRUE;
+ }
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
+//
+// Script to draw a X-section, LOSS, or range made with MakeXsection
+//
void
-DrawXsection(const char* filename="xsec.root",
+DrawXsection(Bool_t scale=kFALSE,
+ const char* filename="xsec.root",
const char* var="LOSS",
const char* medName="FMD_Si$",
Double_t thick=.03,
vb->SetAddress(&value);
Int_t n = tree->GetEntries();
- TDatabasePDG* pdgDb = TDatabasePDG::Instance();
- TParticlePDG* pdgP = pdgDb->GetParticle(pdgName);
- if (!pdgP) {
- std::cerr << "Couldn't find particle " << pdgName << std::endl;
- return;
+ Float_t xscale = 1;
+ Float_t yscale = 1;
+ if (scale) {
+ TDatabasePDG* pdgDb = TDatabasePDG::Instance();
+ TParticlePDG* pdgP = pdgDb->GetParticle(pdgName);
+ if (!pdgP) {
+ std::cerr << "Couldn't find particle " << pdgName << std::endl;
+ return;
+ }
+ Double_t m = pdgP->Mass();
+ Double_t q = pdgP->Charge() / 3;
+ if (m == 0 || q == 0) {
+ std::cerr << "Mass is 0" << std::endl;
+ return;
+ }
+ xscale = 1 / m;
+ yscale = 1 / (q * q);
}
- // Double_t m = pdgP->Mass();
- // Double_t q = pdgP->Charge() / 3;
- // std::cout << q << "\t" << m << std::endl;
- // if (m == 0) {
- /// std::cerr << "Mass is 0" << std::endl;
- // return;
- // }
- TGraph* graph = new TGraph(n);
+ TGraphErrors* graph = new TGraphErrors(n);
for (Int_t i = 0; i < n; i++) {
tree->GetEntry(i);
- graph->SetPoint(i, tkine, value*thick); // /(m*q*q)
+ Double_t x = tkine*xscale;
+ Double_t y = value*yscale;
+ graph->SetPoint(i, x, y);
+ // 5 sigma
+ graph->SetPointError(i, 0, 5 * .1 * y);
}
- graph->Draw("LP same");
+ graph->SetLineWidth(2);
+ graph->SetFillStyle(3001);
+ graph->SetFillColor(6);
+ graph->Draw("L3 same");
}
//____________________________________________________________________
--- /dev/null
+//
+// $Id$
+//
+// Script that contains a class to get the media where a certain track
+// was created. This is used for background studies.
+//
+// Use the script `Compile.C' to compile this class using ACLic.
+//
+#include <TH2D.h>
+#include <AliFMDHit.h>
+#include <AliFMDInput.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TArrayI.h>
+#include <AliRun.h>
+#include <TNtuple.h>
+
+struct Media : public TNamed
+{
+ TArrayI* fMeds;
+ TNtuple* fCount;
+ 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");
+ }
+ Media(const char* name)
+ : TNamed(name, "Media information"), fMeds(0), fCount(0)
+ {
+ AliDetector* det = gAlice->GetDetector(name);
+ if (!det) {
+ Warning("Media", "Detector %s not found in gAlice", name);
+ return;
+ }
+ fMeds = det->GetIdtmed();
+ if (!fMeds) {
+ Warning("Media", "No mediums for detector %s found", name);
+ return;
+ }
+ fCount = new TNtuple(GetName(), "E/I:DET:RNG:SEC:STR:X/F:Y:Z:EDEP:PDG/I");
+ }
+};
+
+
+class GetMedia : public AliFMDInputHits
+{
+private:
+ TString fModList;
+ TObjArray fMedia;
+ Media* fOther;
+ Media* fAll;
+ Int_t fEv;
+ TFile* fOutput;
+public:
+ GetMedia(const char* modlist="FMD:ITS:BODY:ABSO:START:PIPE",
+ const char* output="media.root")
+ {
+ fOutput = TFile::Open(output, "RECREATE");
+ fOther = new Media("other", "Unknown media"),
+ fAll = new Media("All", "All media")
+ fEv = 0;
+ AddLoad(kKinematics);
+ }
+ Media* FindMedia(Int_t med)
+ {
+ 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;
+ }
+ return 0;
+ }
+ Bool_t Init()
+ {
+ if (!gGeoManager) {
+ Error("GetMedia", "No geometry defined - make sure you have that");
+ return kFALSE;
+ }
+ if (!gAlice) {
+ Error("GetMedia", "gAlice not defined");
+ return kFALSE;
+ }
+ Int_t now = 0;
+ Int_t colon;
+ while ((colon = fModList.Index(":", now)) >= 0) {
+ fMedia.Add(new Media(fModList(now, colon-now)));
+ now = colon+1;
+ }
+ if (now < fModList.Length())
+ fMedia.Add(new Media(now, fModList.Length()-now));
+ if (fMedia.GetEntries() <= 0) return kFALSE;
+ return AliFMDInputHits::Init();
+ }
+ Bool_t Begin(Int_t ev)
+ {
+ fEv = ev;
+ return AliFMDInputHits::Begin(ev);
+ }
+ Bool_t ProcessHit(AliFMDHit* hit, TParticle* track)
+ {
+ if (!hit || !track) {
+ std::cout << "No hit or track (hit: " << hit << ", track: "
+ << track << ")" << std::endl;
+ return kFALSE;
+ }
+ // 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;
+
+ // Get volume
+ TGeoVolume* prodVol = prodNode->GetVolume();
+ if (!prodVol) return kTRUE;
+
+ // Med medium
+ TGeoMedium* prodMed = prodVol->GetMedium();
+ if (!prodMed) return kTRUE;
+
+ 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());
+ return kTRUE;
+ }
+ Bool_t Finish()
+ {
+ fOutput->Write();
+ fOutput->Close();
+ return kTRUE;
+ }
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff