&& hit->Sector() == sector
&& hit->Strip() == strip
&& hit->Track() == track) {
- Warning("AddHit", "already had a hit in FMD%d%c[%2d,%3d] for track # %d,"
- " adding energy (%f) to that hit (%f) -> %f",
- detector, ring, sector, strip, track, edep, hit->Edep(),
- hit->Edep() + edep);
+ AliDebug(1, Form("already had a hit in FMD%d%c[%2d,%3d] for track # %d,"
+ " adding energy (%f) to that hit (%f) -> %f",
+ detector, ring, sector, strip, track, edep, hit->Edep(),
+ hit->Edep() + edep));
hit->SetEdep(hit->Edep() + edep);
return hit;
}
#include <AliFMDHit.h>
#include "AliFMDv0.h"
#include "AliFMDAlla.h"
+#include "AliFMDDigitizerAlla.h"
#include "AliMagF.h"
#include "AliRun.h"
#include "AliMC.h"
if(id==fIdSens1||id==fIdSens2||id==fIdSens3||id==fIdSens4||id==fIdSens5) {
if(gMC->IsTrackEntering()) {
- vol[2]=copy;
+ vol[2]=copy-1;
gMC->CurrentVolOffID(1,copy1);
vol[1]=copy1;
gMC->CurrentVolOffID(2,copy2);
Char_t ring = (vol[0] % 2) == 0 ? 'I' : 'O';
UShort_t sector = vol[1];
UShort_t strip = vol[2];
+ gMC->TrackPosition(pos);
+ TVector3 cur(pos.Vect());
+ TVector3 old(hits[0], hits[1], hits[2]);
+ cur -= old;
+ Double_t len = cur.Mag();
AliFMDHit* h = new(lhits[fNhits++])
AliFMDHit(fIshunt,
gAlice->GetMCApp()->GetCurrentTrackNumber(),
detector, ring, sector, strip,
hits[0], hits[1], hits[2], hits[3], hits[4], hits[5],
- hits[6], iPart, hits[8]);
+ hits[6], iPart, hits[8], len,
+ gMC->IsTrackDisappeared()||gMC->IsTrackStop());
if (hits[6] > 1 && p/mass > 1) fBad->Add(h);
} // IsTrackExiting()
}
if (Edep>0)
charge=Int_t(gRandom->Gaus(chargeOnly,500));
}
+
+//____________________________________________________________________
+AliDigitizer*
+AliFMDAlla::CreateDigitizer(AliRunDigitizer* manager) const
+{
+ // Create a digitizer object
+ AliFMDDigitizerAlla* digitizer = new AliFMDDigitizerAlla(manager);
+ return digitizer;
+}
//--------------------------------------------------------------------------
//
// EOF
// virtual void Hit2Digits(Int_t bgrEvent, Option_t *opt1=" ",
// Option_t *opt2=" ",Text_t *name=" "); // hit to digit for v1 :test
virtual void Response( Float_t Edep);
-//private:
+ //private:
//Int_t fCharge;
-
-
+ AliDigitizer* CreateDigitizer(AliRunDigitizer* manager) const;
protected:
Int_t fIdSens1; // Sensetive volume in FMD
Int_t fIdSens2; // Sensetive volume in FMD
};
#endif
-
+//
+// Local Variables:
+// mode: C++
+// End:
+//
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, 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$ */
+
+//____________________________________________________________________
+//
+//
+//
+#include "AliFMDCalibSampleRate.h" // ALIFMDCALIBGAIN_H
+#include "AliFMDParameters.h" // ALIFMDPARAMETERS_H
+
+//____________________________________________________________________
+ClassImp(AliFMDCalibSampleRate)
+#if 0
+ ; // This is here to keep Emacs for indenting the next line
+#endif
+
+//____________________________________________________________________
+AliFMDCalibSampleRate::AliFMDCalibSampleRate()
+ : fRates(3)
+{
+ fRates.Reset(0);
+}
+
+//____________________________________________________________________
+AliFMDCalibSampleRate::AliFMDCalibSampleRate(const AliFMDCalibSampleRate& o)
+ : TObject(o), fRates(o.fRates)
+{}
+
+//____________________________________________________________________
+AliFMDCalibSampleRate&
+AliFMDCalibSampleRate::operator=(const AliFMDCalibSampleRate& o)
+{
+ fRates = o.fRates;
+ return (*this);
+}
+
+//____________________________________________________________________
+void
+AliFMDCalibSampleRate::Set(UShort_t ddl, UShort_t rate)
+{
+ if (ddl - AliFMDParameters::kBaseDDL < 0) return;
+ fRates[ddl - AliFMDParameters::kBaseDDL] = rate;
+}
+
+//____________________________________________________________________
+UShort_t
+AliFMDCalibSampleRate::Rate(UShort_t ddl) const
+{
+ if (ddl - AliFMDParameters::kBaseDDL < 0) return 0;
+ return fRates[ddl - AliFMDParameters::kBaseDDL];
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+#ifndef ALIFMDCALIBSAMPLERATE_H
+#define ALIFMDCALIBSAMPLERATE_H
+/* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights
+ * reserved.
+ *
+ * See cxx source for full Copyright notice
+ */
+#ifndef ROOT_TObject
+# include <TObject.h>
+#endif
+#ifndef ROOT_TArrayI
+# include <TArrayI.h>
+#endif
+//____________________________________________________________________
+//
+// Gain value and width for each strip in the FMD
+//
+class AliFMDCalibSampleRate : public TObject
+{
+public:
+ AliFMDCalibSampleRate();
+ AliFMDCalibSampleRate(const AliFMDCalibSampleRate& o);
+ AliFMDCalibSampleRate& operator=(const AliFMDCalibSampleRate& o);
+ void Set(UShort_t ddl, UShort_t rate);
+ UShort_t Rate(UShort_t ddl) const;
+protected:
+ TArrayI fRates; // Sample rates
+ ClassDef(AliFMDCalibSampleRate,1); // Sample rates
+};
+
+#endif
+//____________________________________________________________________
+//
+// Local Variables:
+// mode: C++
+// End:
+//
+
+
#include "AliFMDDigit.h" // ALIFMDDIGIT_H
#include "Riostream.h" // ROOT_Riostream
+#include <TString.h>
//====================================================================
ClassImp(AliFMDBaseDigit)
#if 0
- ; // This is here to keep Emacs for indenting the next line
+ ; // This is here to keep Emacs from indenting the next line
#endif
//____________________________________________________________________
<< flush;
}
+//____________________________________________________________________
+const char*
+AliFMDBaseDigit::GetName() const
+{
+ return Form("FMD%d%c[%2d,%3d]", fDetector, fRing, fSector, fStrip);
+}
+
//====================================================================
ClassImp(AliFMDDigit)
UShort_t Sector() const { return fSector; }
UShort_t Strip() const { return fStrip; }
virtual void Print(Option_t* opt="") const;
-
+ const char* GetName() const;
protected:
UShort_t fDetector; // (Sub) Detector # (1,2, or 3)
Char_t fRing; // Ring ID ('I' or 'O')
#include "AliFMDRing.h" // ALIFMDRING_H
#include "AliFMDHit.h" // ALIFMDHIT_H
#include "AliFMDDigit.h" // ALIFMDDIGIT_H
+#include "AliFMDParameters.h" // ALIFMDPARAMETERS_H
#include <AliRunDigitizer.h> // ALIRUNDIGITIZER_H
#include <AliRun.h> // ALIRUN_H
#include <AliLoader.h> // ALILOADER_H
{
// Normal CTOR
AliDebug(1," processed");
- SetVA1MipRange();
- SetAltroChannelSize();
- SetSampleRate();
SetShapingTime();
}
{
// Normal CTOR
AliDebug(1," processed");
- SetVA1MipRange();
- SetAltroChannelSize();
- SetSampleRate();
SetShapingTime();
}
//____________________________________________________________________
+UShort_t
+AliFMDBaseDigitizer::MakePedestal(UShort_t,
+ Char_t,
+ UShort_t,
+ UShort_t) const
+{
+ return 0;
+}
+
+//____________________________________________________________________
void
AliFMDBaseDigitizer::SumContributions(AliFMD* fmd)
{
// Get number of entries in the tree
Int_t ntracks = Int_t(hitsTree->GetEntries());
+ AliFMDParameters* param = AliFMDParameters::Instance();
Int_t read = 0;
// Loop over the tracks in the
for (Int_t track = 0; track < ntracks; track++) {
UShort_t sector = fmdHit->Sector();
UShort_t strip = fmdHit->Strip();
Float_t edep = fmdHit->Edep();
+
+ // Check if strip is `dead'
+ if (param->IsDead(detector, ring, sector, strip)) continue;
+
+ // Give warning in case of double hit
if (fEdep(detector, ring, sector, strip).fEdep != 0)
AliDebug(5, Form("Double hit in %d%c(%d,%d)",
detector, ring, sector, strip));
+ // Sum energy deposition
fEdep(detector, ring, sector, strip).fEdep += edep;
fEdep(detector, ring, sector, strip).fN += 1;
// Add this to the energy deposited for this strip
AliFMDDetector* det = geometry->GetDetector(detector);
if (!det) continue;
for (UShort_t ringi = 0; ringi <= 1; ringi++) {
+ Char_t ring = ringi == 0 ? 'I' : 'O';
// Get pointer to Ring
- AliFMDRing* r = det->GetRing((ringi == 0 ? 'I' : 'O'));
+ AliFMDRing* r = det->GetRing(ring);
if (!r) continue;
// Get number of sectors
// VA1_ALICE channel.
if (strip % 128 == 0) last = 0;
- Float_t edep = fEdep(detector, r->GetId(), sector, strip).fEdep;
- ConvertToCount(edep, last, r->GetSiThickness(),
- geometry->GetSiDensity(), counts);
+ Float_t edep = fEdep(detector, ring, sector, strip).fEdep;
+ ConvertToCount(edep, last, detector, ring, sector, strip, counts);
last = edep;
- AddDigit(fmd, detector, r->GetId(), sector, strip,
- edep, UShort_t(counts[0]),
- Short_t(counts[1]), Short_t(counts[2]));
+ AddDigit(fmd, detector, ring, sector, strip, edep,
+ UShort_t(counts[0]), Short_t(counts[1]),
+ Short_t(counts[2]));
#if 0
// This checks if the digit created will give the `right'
// number of particles when reconstructed, using a naiive
// approach. It's here only as a quality check - nothing
// else.
- CheckDigit(fEdep(detector, r->GetId(), sector, strip).fEdep,
- fEdep(detector, r->GetId(), sector, strip).fN,
+ CheckDigit(digit, fEdep(detector, ring, sector, strip).fN,
counts);
#endif
} // Strip
void
AliFMDBaseDigitizer::ConvertToCount(Float_t edep,
Float_t last,
- Float_t siThickness,
- Float_t siDensity,
+ UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip,
TArrayI& counts) const
{
// Convert the total energy deposited to a (set of) ADC count(s).
//
// = E + (l - E) * ext(-B * t)
//
- Float_t mipEnergy = 1.664 * siThickness * siDensity;
- Float_t convF = (1/mipEnergy*Float_t(fAltroChannelSize)/fVA1MipRange);
- UShort_t ped = MakePedestal();
-
+ AliFMDParameters* param = AliFMDParameters::Instance();
+ Float_t convF = 1/param->GetPulseGain(detector,ring,sector,strip);
+ UShort_t ped = MakePedestal(detector,ring,sector,strip);
+ UInt_t maxAdc = param->GetAltroChannelSize();
+ UShort_t rate = param->GetSampleRate(AliFMDParameters::kBaseDDL);
+ UShort_t size = param->GetAltroChannelSize();
+
// In case we don't oversample, just return the end value.
- if (fSampleRate == 1) {
- counts[0] = UShort_t(TMath::Min(edep * convF + ped,
- Float_t(fAltroChannelSize)));
+ if (rate == 1) {
+ counts[0] = UShort_t(TMath::Min(edep * convF + ped, Float_t(size)));
return;
}
// Create a pedestal
- Int_t n = fSampleRate;
Float_t b = fShapingTime;
- for (Ssiz_t i = 0; i < n; i++) {
- Float_t t = Float_t(i) / n;
+ for (Ssiz_t i = 0; i < rate; i++) {
+ Float_t t = Float_t(i) / rate;
Float_t s = edep + (last - edep) * TMath::Exp(-b * t);
- counts[i] = UShort_t(TMath::Min(s * convF + ped,
- Float_t(fAltroChannelSize)));
+ counts[i] = UShort_t(TMath::Min(s * convF + ped, Float_t(maxAdc)));
}
}
{
// Normal CTOR
AliDebug(1," processed");
- SetPedestal();
}
//____________________________________________________________________
if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
// Get the AliFMD object
- AliFMD* fmd = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
+ AliFMD* fmd =
+ static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
if (!fmd) {
AliError("Can not get FMD from gAlice");
return;
//____________________________________________________________________
UShort_t
-AliFMDDigitizer::MakePedestal() const
+AliFMDDigitizer::MakePedestal(UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip) const
{
// Make a pedestal
- return UShort_t(TMath::Max(gRandom->Gaus(fPedestal, fPedestalWidth), 0.));
+ AliFMDParameters* param =AliFMDParameters::Instance();
+ Float_t mean =param->GetPedestal(detector,ring,sector,strip);
+ Float_t width =param->GetPedestalWidth(detector,ring,sector,strip);
+ return UShort_t(TMath::Max(gRandom->Gaus(mean, width), 0.));
}
//____________________________________________________________________
//____________________________________________________________________
void
-AliFMDDigitizer::CheckDigit(Float_t /* edep */,
+AliFMDDigitizer::CheckDigit(AliFMDDigit* digit,
UShort_t nhits,
const TArrayI& counts)
{
// Check that digit is consistent
- Int_t integral = counts[0];
+ AliFMDParameters* param = AliFMDParameters::Instance();
+ UShort_t det = digit->Detector();
+ Char_t ring = digit->Ring();
+ UShort_t sec = digit->Sector();
+ UShort_t str = digit->Strip();
+ Float_t mean = param->GetPedestal(det,ring,sec,str);
+ Float_t width = param->GetPedestalWidth(det,ring,sec,str);
+ UShort_t range = param->GetVA1MipRange();
+ UShort_t size = param->GetAltroChannelSize();
+ 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);
+ integral -= Int_t(mean + 2 * width);
if (integral < 0) integral = 0;
- Float_t convF = Float_t(fVA1MipRange) / fAltroChannelSize;
+ Float_t convF = Float_t(range) / size;
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",
class AliFMD;
class AliLoader;
class AliRunLoader;
+class AliFMDDigit;
//====================================================================
virtual Bool_t Init();
// Extra member functions
- void SetVA1MipRange(UShort_t r=20) { fVA1MipRange = r; }
- void SetAltroChannelSize(UShort_t s=1024) { fAltroChannelSize = s;}
- void SetSampleRate(UShort_t r=1) { fSampleRate = (r>2 ? 2 : r); }
- void SetShapingTime(Float_t t=10) { fShapingTime = t; }
-
- UShort_t GetVA1MipRange() const { return fVA1MipRange; }
- UShort_t GetAltroChannelSize() const { return fAltroChannelSize; }
- UShort_t GetSampleRate() const { return fSampleRate; }
+ void SetShapingTime(Float_t t=10) { fShapingTime = t; }
Float_t GetShapingTime() const { return fShapingTime; }
protected:
virtual void SumContributions(AliFMD* fmd);
virtual void DigitizeHits(AliFMD* fmd) const;
virtual void ConvertToCount(Float_t edep,
Float_t last,
- Float_t siThickness,
- Float_t siDensity,
+ UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip,
TArrayI& counts) const;
- virtual UShort_t MakePedestal() const { return 0; }
+ virtual UShort_t MakePedestal(UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip) const;
virtual void AddNoise(TArrayI&) const {}
virtual void AddDigit(AliFMD* /* fmd */,
UShort_t /* detector */,
AliRunLoader* fRunLoader; // Run loader
AliFMDEdepMap fEdep; // Cache of Energy from hits
- UShort_t fVA1MipRange; // How many MIPs the pre-amp can do
- UShort_t fAltroChannelSize; // Largest # to store in 1 ADC chan.
- 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
+ ClassDef(AliFMDBaseDigitizer,1) // Base class for FMD digitizers
};
//====================================================================
AliFMDDigitizer(AliRunDigitizer * manager);
virtual ~AliFMDDigitizer() {}
virtual void Exec(Option_t* option=0);
-
-
- // Extra member functions
- void SetPedestal(Float_t mean=10, Float_t width=.5);
- void GetPedestal(Float_t& mean, Float_t& width) const;
protected:
virtual void AddDigit(AliFMD* fmd,
UShort_t detector,
UShort_t count1,
Short_t count2,
Short_t count3) const;
- virtual UShort_t MakePedestal() const;
- virtual void CheckDigit(Float_t edep,
+ virtual UShort_t MakePedestal(UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip) const;
+ virtual void CheckDigit(AliFMDDigit* digit,
UShort_t nhits,
const TArrayI& counts);
- Float_t fPedestal; // Mean of pedestal
- Float_t fPedestalWidth; // Width of pedestal
- ClassDef(AliFMDDigitizer,0) // Make Digits from Hits
+ ClassDef(AliFMDDigitizer,1) // Make Digits from Hits
};
-//____________________________________________________________________
-inline void
-AliFMDDigitizer::SetPedestal(Float_t mean, Float_t width)
-{
- fPedestal = mean;
- fPedestalWidth = width;
-}
-
-//____________________________________________________________________
-inline void
-AliFMDDigitizer::GetPedestal(Float_t& mean, Float_t& width) const
-{
- mean = fPedestal;
- width = fPedestalWidth;
-}
-
//====================================================================
class AliFMDSDigitizer : public AliFMDBaseDigitizer
--- /dev/null
+ /**************************************************************************
+ * Copyright(c) 1998-2000, 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. *
+ **************************************************************************/
+
+ //////////////////////////////////////////////////////////////////////////////
+// //
+// Forward Multiplicity Detector based on Silicon plates //
+// This class contains the procedures simulation ADC signal for //
+// the Forward Multiplicity detector : hits -> digits //
+// ADC signal consists //
+// - number of detector; //
+// - number of ring; //
+// - number of sector; //
+// - ADC signal in this channel //
+// //
+ //////////////////////////////////////////////////////////////////////////////
+
+#include <TTree.h>
+#include <TVector.h>
+#include <TObjArray.h>
+#include <TFile.h>
+#include <TDirectory.h>
+#include <TRandom.h>
+#include "AliLog.h"
+#include "AliFMDDigitizerAlla.h"
+#include "AliFMD.h"
+#include "AliFMDHit.h"
+#include "AliFMDDigit.h"
+#include "AliRunDigitizer.h"
+#include "AliRun.h"
+#include "AliLoader.h"
+#include "AliRunLoader.h"
+#include <stdlib.h>
+#include <Riostream.h>
+
+ClassImp(AliFMDDigitizerAlla)
+#if 0
+ ;
+#endif
+
+//___________________________________________
+AliFMDDigitizerAlla::AliFMDDigitizerAlla()
+ : AliDigitizer()
+{
+ // Default ctor - don't use it
+}
+
+//___________________________________________
+AliFMDDigitizerAlla::AliFMDDigitizerAlla(AliRunDigitizer* manager)
+ : AliDigitizer(manager)
+{
+ // ctor which should be used
+ // fDebug =0;
+ AliDebug(1," processed");
+}
+
+//------------------------------------------------------------------------
+AliFMDDigitizerAlla::~AliFMDDigitizerAlla()
+{
+ // Destructor
+}
+
+ //------------------------------------------------------------------------
+Bool_t AliFMDDigitizerAlla::Init()
+{
+ // Initialization
+ // cout<<"AliFMDDigitizerAlla::Init"<<endl;
+ return kTRUE;
+}
+
+
+//---------------------------------------------------------------------
+
+void AliFMDDigitizerAlla::Exec(Option_t * /*option*/)
+{
+
+ // Conver hits to digits:
+ // - number of detector;
+ // - number of ring;
+ // - number of sector;
+ // - ADC signal in this channel
+ AliDebug(1," start...");
+ Float_t de[3][2][50][520];
+ for (Int_t i=0; i<3; i++)
+ for(Int_t j=0; j<2; j++)
+ for(Int_t k=0; k < 50; k++)
+ for (Int_t l=0; l < 520; l++)
+ de[i][j][k][l]=0;
+ Int_t nstrips[] = { 512, 256 };
+ Int_t nsectors[] = { 20, 40 };
+
+ AliRunLoader* outRL =
+ AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
+ AliLoader* pOutFMD = outRL->GetLoader("FMDLoader");
+
+ AliRunLoader* inRL =
+ AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
+
+ if (!inRL) {
+ AliError("Can not find Run Loader for input stream 0");
+ return;
+ }
+ if (!inRL->GetAliRun()) inRL->LoadgAlice();
+
+ AliFMD * fFMD = static_cast<AliFMD*>(inRL->GetAliRun()->GetDetector("FMD"));
+ if (!fFMD) {
+ AliError("Can not get FMD from gAlice");
+ return;
+ }
+
+ // Loop over files to digitize
+ Int_t nFiles = GetManager()->GetNinputs();
+ for (Int_t inputFile=0; inputFile < nFiles; inputFile++) {
+ AliDebug(1,Form("Digitizing event number %d",
+ fManager->GetOutputEventNr()));
+ if (fFMD) {
+ AliRunLoader* inRL =
+ AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
+ AliLoader* pInFMD = inRL->GetLoader("FMDLoader");
+ pInFMD->LoadHits("READ");
+
+
+ TTree* tH = pInFMD->TreeH();
+ if (!tH) {
+ pInFMD->LoadHits("read");
+ tH = pInFMD->TreeH();
+ }
+ TBranch* brHits = tH->GetBranch("FMD");
+ if (brHits) fFMD->SetHitsAddressBranch(brHits);
+ else AliFatal("Branch FMD hit not found");
+ TClonesArray *fFMDhits = fFMD->Hits ();
+
+ Int_t ntracks = tH->GetEntries();
+ for (Int_t track = 0; track < ntracks; track++) {
+ brHits->GetEntry(track);
+ Int_t nhits = fFMDhits->GetEntries ();
+ for (Int_t hit = 0; hit < nhits; hit++) {
+ AliFMDHit* fmdHit =
+ static_cast<AliFMDHit*>(fFMDhits->UncheckedAt(hit));
+ Int_t detector = fmdHit->Detector();
+ Int_t iring = fmdHit->Ring() == 'I' ? 0 : 1;
+ Int_t sector = fmdHit->Sector();
+ Int_t strip = fmdHit->Strip();
+ de[detector-1][iring][sector][strip] += fmdHit->Edep();
+ } //hit loop
+ } //track loop
+ } //if FMD
+ }
+
+
+ // Put noise and make ADC signal
+ Float_t mipI = 1.664 * 0.04 * 2.33 / 22400; // = 6.923e-6;
+ for (Int_t detector = 1; detector <= 3; detector++){
+ for (Int_t iring = 0; iring < 2; iring++) {
+ if (detector == 1 && iring == 1) continue;
+ char ring = (iring == 0 ? 'I' : 'O');
+ for (Int_t sector = 0; sector < nsectors[iring]; sector++) {
+ for (Int_t strip = 0; strip < nstrips[iring]; strip++) {
+ Int_t signal = Int_t(de[detector-1][iring][sector][strip] / mipI);
+ Int_t pedestal = Int_t(gRandom->Gaus(500,250));
+ Int_t charge = signal + pedestal;
+ if(charge <= 500) charge = 500;
+ //dynamic range from MIP(0.155MeV) to 30MIP(4.65MeV)
+ //1024 ADC channels
+ Float_t channelWidth = (22400 * 50) / 1024;
+ Int_t adc = Int_t(charge / channelWidth);
+ if (adc > 1023) adc = 1023;
+ fFMD->AddDigitByFields(detector, ring, sector, strip, adc);
+ } //strip
+ } //sector
+ } //iring
+ } // detector
+
+ pOutFMD->LoadDigits("update");
+ TTree* treeD = pOutFMD->TreeD();
+ if (!treeD) {
+ pOutFMD->MakeTree("D");
+ treeD = pOutFMD->TreeD();
+ }
+
+ treeD->Reset();
+ TClonesArray* digits = fFMD->Digits();
+ fFMD->MakeBranchInTree(treeD, fFMD->GetName(), &(digits), 4000, 0);
+ if (!treeD->GetBranch("FMD")) AliFatal("No branch for FMD digits");
+ treeD->Fill();
+ pOutFMD->WriteDigits("OVERWRITE");
+ pOutFMD->UnloadHits();
+ pOutFMD->UnloadDigits();
+ fFMD->ResetDigits();
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+
+
+
+
--- /dev/null
+#ifndef ALIFMDDIGITIZERALLA_H
+#define ALIFMDDIGITIZERALLA_H
+/* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+//
+// Resurection of Alla's old digitizer
+// This is to investigate the changes we've seen in the reconstructed
+// particle multiplicity
+//
+#include <AliDigitizer.h>
+#include <AliRunDigitizer.h>
+class TClonesArray;
+
+class AliFMDDigitizerAlla : public AliDigitizer
+{
+public:
+ AliFMDDigitizerAlla();
+ AliFMDDigitizerAlla(AliRunDigitizer * manager);
+ virtual ~AliFMDDigitizerAlla();
+ virtual Bool_t Init();
+ // Do the main work
+ void Exec(Option_t* option=0) ;
+ Int_t PutNoise(Int_t charge) {return (Int_t)(gRandom->Gaus(charge,500));}
+ TClonesArray *Digits() const {return fDigits;}
+ TClonesArray *Hits() const {return fHits;}
+ enum {kBgTag = -1};
+private:
+ TClonesArray *fDigits; // ! array with digits
+ TClonesArray *fHits; // List of hits
+ AliRunDigitizer* GetManager(){return fManager;}
+ ClassDef(AliFMDDigitizerAlla,0)
+};
+#endif
+//
+// Local Variables:
+// mode: C++
+// End:
+//
+
+
+
+
+++ /dev/null
-/**************************************************************
- * Copyright(c) 1998-1999, 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$ */
-//__________________________________________________________
-//
-// Map of per strip Float_t information
-//
-// Created Mon Nov 8 12:51:51 2004 by Christian Holm Christensen
-//
-#include "AliFMDFloatMap.h" //ALIFMDFLOATMAP_H
-//__________________________________________________________
-ClassImp(AliFMDFloatMap)
-#if 0
- ; // This is here to keep Emacs for indenting the next line
-#endif
-//__________________________________________________________
-AliFMDFloatMap::AliFMDFloatMap(const AliFMDFloatMap& other)
- : AliFMDMap(other.fMaxDetectors,
- other.fMaxRings,
- other.fMaxSectors,
- other.fMaxStrips),
- fData(0)
-{
- // Copy constructor
- fData = new Float_t[fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips];
- for (size_t i = 0; i < fMaxDetectors * fMaxRings
- * fMaxSectors * fMaxStrips; i++)
- fData[i] = other.fData[i];
-}
-
-//__________________________________________________________
-AliFMDFloatMap::AliFMDFloatMap(size_t maxDet,
- size_t maxRing,
- size_t maxSec,
- size_t maxStr)
- : AliFMDMap(maxDet, maxRing, maxSec, maxStr),
- fData(0)
-{
- // Constructor.
- // Parameters:
- // maxDet Maximum number of detectors
- // maxRing Maximum number of rings per detector
- // maxSec Maximum number of sectors per ring
- // maxStr Maximum number of strips per sector
- fData = new Float_t[fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips];
- Reset();
-}
-
-//__________________________________________________________
-AliFMDFloatMap&
-AliFMDFloatMap::operator=(const AliFMDFloatMap& other)
-{
- // Assignment operator
- fMaxDetectors = other.fMaxDetectors;
- fMaxRings = other.fMaxRings;
- fMaxSectors = other.fMaxSectors;
- fMaxStrips = other.fMaxStrips;
- if (fData) delete [] fData;
- fData = new Float_t[fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips];
- for (size_t i = 0; i < fMaxDetectors * fMaxRings
- * fMaxSectors * fMaxStrips; i++)
- fData[i] = other.fData[i];
- return *this;
-}
-
-//__________________________________________________________
-void
-AliFMDFloatMap::Reset(const Float_t& val)
-{
- // Reset map to val
- for (size_t i = 0; i < fMaxDetectors * fMaxRings
- * fMaxSectors * fMaxStrips; i++)
- fData[i] = val;
-}
-
-//__________________________________________________________
-Float_t&
-AliFMDFloatMap::operator()(UShort_t det,
- Char_t ring,
- UShort_t sec,
- UShort_t str)
-{
- // Get data
- // Parameters:
- // det Detector #
- // ring Ring ID
- // sec Sector #
- // str Strip #
- // Returns appropriate data
- return fData[CalcIndex(det, ring, sec, str)];
-}
-
-//__________________________________________________________
-const Float_t&
-AliFMDFloatMap::operator()(UShort_t det,
- Char_t ring,
- UShort_t sec,
- UShort_t str) const
-{
- // Get data
- // Parameters:
- // det Detector #
- // ring Ring ID
- // sec Sector #
- // str Strip #
- // Returns appropriate data
- return fData[CalcIndex(det, ring, sec, str)];
-}
-
-//__________________________________________________________
-//
-// EOF
-//
-
+++ /dev/null
-#ifndef ALIFMDFLOATMAP_H
-#define ALIFMDFLOATMAP_H
-/* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights
- * reserved.
- *
- * See cxx source for full Copyright notice
- */
-#ifndef ALIFMDMAP_H
-# include <AliFMDMap.h>
-#endif
-//____________________________________________________________________
-//
-// Array of floats indexed by strip identifier.
-//
-class AliFMDFloatMap : public AliFMDMap
-{
-public:
- AliFMDFloatMap(size_t maxDet = kMaxDetectors,
- size_t maxRing= kMaxRings,
- size_t maxSec = kMaxSectors,
- size_t maxStr = kMaxStrips);
- AliFMDFloatMap(const AliFMDFloatMap& o);
- virtual ~AliFMDFloatMap() { delete [] fData; }
- AliFMDFloatMap& operator=(const AliFMDFloatMap& o);
- virtual void Reset(const Float_t& v=Float_t());
- virtual Float_t& operator()(UShort_t det,
- Char_t ring,
- UShort_t sec,
- UShort_t str);
- virtual const Float_t& operator()(UShort_t det,
- Char_t ring,
- UShort_t sec,
- UShort_t str) const;
-protected:
- Float_t* fData;
- ClassDef(AliFMDFloatMap,1) // Map of floats
-};
-
-#endif
-//____________________________________________________________________
-//
-// Local Variables:
-// mode: C++
-// End:
-//
-
//
#include "AliFMDInput.h" // ALIFMDHIT_H
#include "AliLog.h" // ALILOG_H
-#include "Riostream.h" // ROOT_Riostream
-#include <TDatabasePDG.h>
-#include <TMath.h>
-#include <TString.h>
-#include <TParticle.h>
+#include "AliLoader.h" // ALILOADER_H
+#include "AliRunLoader.h" // ALIRUNLOADER_H
+#include "AliRun.h" // ALIRUN_H
+#include "AliStack.h" // ALISTACK_H
+#include "AliFMD.h" // ALIFMD_H
#include "AliFMDHit.h" // ALIFMDHIT_H
-#include "AliFMDDigit.h" // ALIFMDDigit_H
-#include "AliFMDMultStrip.h" // ALIFMDMultStrip_H
-#include "AliFMDMultRegion.h" // ALIFMDMultRegion_H
+#include "AliFMDDigit.h" // ALIFMDDigit_H
+#include "AliFMDMultStrip.h" // ALIFMDMultStrip_H
+#include "AliFMDMultRegion.h" // ALIFMDMultRegion_H
+#include <TTree.h> // ROOT_TTree
+#include <TParticle.h> // ROOT_TParticle
+#include <TString.h> // ROOT_TString
+#include <TDatabasePDG.h> // ROOT_TDatabasePDG
+#include <TMath.h> // ROOT_TMath
+#include <TGeoManager.h> // ROOT_TGeoManager
+#include <Riostream.h> // ROOT_Riostream
//____________________________________________________________________
ClassImp(AliFMDInput)
return kFALSE;
}
fTreeE = fLoader->TreeE();
+
+ // Optionally, get the geometry
+ if (TESTBIT(fTreeMask, kGeometry)) {
+ TString fname(fRun->GetGeometryFileName());
+ if (fname.IsNull()) {
+ Warning("Init", "No file name for the geometry from AliRun");
+ fname = gSystem->DirName(fGAliceFile);
+ fname.Append("/geometry.root");
+ }
+ fGeoManager = TGeoManager::Import(fname.Data());
+ if (!fGeoManager) {
+ Fatal("Init", "No geometry manager found");
+ return kFALSE;
+ }
+ }
fIsInit = kTRUE;
return fIsInit;
if (fFMDLoader->LoadRecPoints()) return kFALSE;
fTreeR = fFMDLoader->TreeR();
if (!fArrayN) fArrayN = new TClonesArray("AliFMDMultStrip");
- if (!fArrayP) fArrayP = new TClonesArray("AliFMDMultRegion");
- fTreeR->SetBranchAddress("FMDNaiive", &fArrayN);
- fTreeR->SetBranchAddress("FMDPoisson", &fArrayP);
+ // if (!fArrayP) fArrayP = new TClonesArray("AliFMDMultRegion");
+ fTreeR->SetBranchAddress("FMD", &fArrayN);
+ // fTreeR->SetBranchAddress("FMDPoisson", &fArrayP);
}
return kTRUE;
}
// classes for customized class that do some sort of analysis on the
// various types of data produced by the FMD.
//
-#ifndef ALILOADER_H
-# include <AliLoader.h>
-#endif
-#ifndef ALIRUNLOADER_H
-# include <AliRunLoader.h>
-#endif
-#ifndef ALIRUN_H
-# include <AliRun.h>
-#endif
-#ifndef ALISTACK_H
-# include <AliStack.h>
-#endif
-#ifndef ALIFMD_H
-# include <AliFMD.h>
-#endif
-#ifndef ROOT_TTree
-# include <TTree.h>
-#endif
-#ifndef ROOT_TParticle
-# include <TParticle.h>
-#endif
+#include <TObject.h>
#ifndef ROOT_TString
# include <TString.h>
#endif
+class AliRunLoader;
+class AliLoader;
+class AliStack;
+class AliRun;
+class AliFMD;
+class AliFMDHit;
+class TString;
+class TClonesArray;
+class TTree;
+class TGeoManager;
+class TParticle;
+
//___________________________________________________________________
class AliFMDInput : public TObject
{
public:
enum ETrees {
- kHits = 1,
- kKinematics,
- kDigits,
- kSDigits,
- kHeader,
- kRecPoints
+ kHits = 1, // Hits
+ kKinematics, // Kinematics (from sim)
+ kDigits, // Digits
+ kSDigits, // Summable digits
+ kHeader, // Header information
+ kRecPoints, // Reconstructed points
+ kGeometry // Not really a tree
};
AliFMDInput();
AliFMDInput(const char* gAliceFile);
TClonesArray* fArrayS; // SDigit info array
TClonesArray* fArrayN; // Mult (single) info array
TClonesArray* fArrayP; // Mult (region) info array
+ TGeoManager* fGeoManager; // Geometry manager
Int_t fTreeMask; // Which tree's to load
Bool_t fIsInit;
ClassDef(AliFMDInput,0) //Hits for detector FMD
+++ /dev/null
-/**************************************************************************
- * Copyright(c) 1998-1999, 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$ */
-
-//____________________________________________________________________
-//
-// 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)
-#if 0
- ; // This is here to keep Emacs for indenting the next line
-#endif
-
-//____________________________________________________________________
-AliFMDMap::AliFMDMap(size_t maxDet,
- size_t maxRing,
- size_t maxSec,
- size_t maxStr)
- : fMaxDetectors(maxDet),
- fMaxRings(maxRing),
- fMaxSectors(maxSec),
- fMaxStrips(maxStr)
-{
- // Construct a map
- //
- // Parameters:
- // maxDet Maximum # of detectors
- // maxRinf Maximum # of rings
- // maxSec Maximum # of sectors
- // maxStr Maximum # of strips
-}
-
-//____________________________________________________________________
-Int_t
-AliFMDMap::CheckIndex(size_t det, Char_t ring, size_t sec, size_t str) const
-{
- // Check that the index supplied is OK. Returns true index, or -1
- // on error.
- size_t ringi = (ring == 'I' || ring == 'i' ? 0 : 1);
- size_t idx =
- (det + fMaxDetectors * (ringi + fMaxRings * (sec + fMaxSectors * str)));
- if (idx >= fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips)
- return -1;
- return idx;
-}
-
-
-//____________________________________________________________________
-size_t
-AliFMDMap::CalcIndex(size_t det, Char_t ring, size_t sec, size_t str) const
-{
- // Calculate index into storage from arguments.
- //
- // Parameters:
- // det Detector #
- // ring Ring ID
- // sec Sector #
- // str Strip #
- //
- // Returns appropriate index into storage
- //
- Int_t idx = CheckIndex(det, ring, sec, str);
- if (idx < 0) {
- size_t ringi = (ring == 'I' || ring == 'i' ? 0 : 1);
- 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);
-}
-
-
-//___________________________________________________________________
-//
-// EOF
-//
+++ /dev/null
-#ifndef ALIFMDMAP_H
-#define ALIFMDMAP_H
-/* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights
- * reserved.
- *
- * See cxx source for full Copyright notice
- */
-#ifndef ROOT_TObject
-# include <TObject.h>
-#endif
-//____________________________________________________________________
-//
-// 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:
- enum {
- kMaxDetectors = 3,
- kMaxRings = 2,
- kMaxSectors = 40,
- kMaxStrips = 512
- };
- AliFMDMap(size_t maxDet = kMaxDetectors,
- size_t maxRing= kMaxRings,
- size_t maxSec = kMaxSectors,
- size_t maxStr = kMaxStrips);
- virtual ~AliFMDMap() {}
- Int_t CheckIndex(size_t det, Char_t ring, size_t sec, size_t str) const;
-protected:
- size_t CalcIndex(size_t det, Char_t ring, size_t sec, size_t str) const;
- size_t fMaxDetectors; // Maximum # of detectors
- size_t fMaxRings; // Maximum # of rings
- size_t fMaxSectors; // Maximum # of sectors
- size_t fMaxStrips; // Maximum # of strips
- ClassDef(AliFMDMap, 1) // Cache of per strip information
-};
-
-#ifdef MAY_USE_TEMPLATES
-#ifndef __VECTOR__
-# include <vector>
-#endif
-//____________________________________________________________________
-//
-// Class template for classes that cache per strip information.
-// Access to the data goes via
-//
-// Type& AliFMDMap<Type>::operator()(size_t detector,
-// Char_t ring,
-// size_t sector,
-// size_t strip);
-//
-// (as well as a const version of this member function).
-// The elements can be reset to the default value by calling
-// AliFMDMap<Type>::Clear(). This resets the values to `Type()'.
-//
-template <typename Type>
-class AliFMDMap : public TObject
-{
-public:
- AliFMDMap(size_t maxDet=3, size_t maxRing=2, size_t maxSec=40,
- size_t maxStr=512);
- virtual ~AliFMDMap() {}
- void Clear(const Type& val=Type());
- Type& operator()(size_t det, Char_t ring, size_t sec, size_t str);
- const Type& operator()(size_t det, Char_t ring, size_t sec, size_t str)const;
-private:
- typedef std::vector<Type> ValueVector; // Type of container
- ValueVector fValues; // Contained values
- size_t fMaxDetectors; // Maximum # of detectors
- size_t fMaxRings; // Maximum # of rings
- size_t fMaxSectors; // Maximum # of sectors
- size_t fMaxStrips; // Maximum # of strips
-
- size_t CalcIndex(size_t det, Char_t ring, size_t sec, size_t str) const;
- ClassDef(AliFMDMap, 0); // Map of FMD index's to values
-};
-
-
-//____________________________________________________________________
-template <typename Type>
-inline
-AliFMDMap<Type>::AliFMDMap(size_t maxDet,
- size_t maxRing,
- size_t maxSec,
- size_t maxStr)
- : fValues(maxDet * maxRing * maxSec * maxStr),
- fMaxDetectors(maxDet),
- fMaxRings(maxRing),
- fMaxSectors(maxSec),
- fMaxStrips(maxStr)
-{
- // Construct a map
- //
- // Parameters:
- // maxDet Maximum # of detectors
- // maxRinf Maximum # of rings
- // maxSec Maximum # of sectors
- // maxStr Maximum # of strips
-}
-
-
-//____________________________________________________________________
-template <typename Type>
-inline size_t
-AliFMDMap<Type>::CalcIndex(size_t det, Char_t ring, size_t sec, size_t str) const
-{
- // Calculate index into storage from arguments.
- //
- // Parameters:
- // det Detector #
- // ring Ring ID
- // sec Sector #
- // str Strip #
- //
- // Returns appropriate index into storage
- //
- size_t ringi = (ring == 'I' || ring == 'i' ? 0 : 1);
- size_t idx =
- (det + fMaxDetectors * (ringi + fMaxRings * (sec + fMaxSectors * str)));
- if (idx >= fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips) {
- 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"))));
- return 0;
- }
- return idx;
-}
-
-//____________________________________________________________________
-template <typename Type>
-inline void
-AliFMDMap<Type>::Clear(const Type& val)
-{
- // Resets stored values to the default value for that type
- for (size_t i = 0; i < fValues.size(); ++i) fValues[i] = val;
-}
-
-//____________________________________________________________________
-template <typename Type>
-inline Type&
-AliFMDMap<Type>::operator()(size_t det, Char_t ring, size_t sec, size_t str)
-{
- // Parameters:
- // det Detector #
- // ring Ring ID
- // sec Sector #
- // str Strip #
- //
- // Returns data[det][ring][sec][str]
- return fValues[CalcIndex(det, ring, sec, str)];
-}
-
-//____________________________________________________________________
-template <typename Type>
-inline const Type&
-AliFMDMap<Type>::operator()(size_t det,
- Char_t ring,
- size_t sec,
- size_t str) const
-{
- // Parameters:
- // det Detector #
- // ring Ring ID
- // sec Sector #
- // str Strip #
- //
- // Returns data[det][ring][sec][str]
- return fValues[CalcIndex(det, ring, sec, str)];
-}
-
-
-//____________________________________________________________________
-//
-// Some specialisations
-//
-typedef AliFMDMap<UShort_t> AliFMDAdcMap;
-typedef AliFMDMap<std::pair<Float_t, UShort_t> > AliFMDEdepMap;
-
-#endif
-#endif
-//____________________________________________________________________
-//
-// Local Variables:
-// mode: C++
-// End:
-//
-// EOF
-//
-
-
AliFMDMultPoisson::PreEvent(TTree* tree, Float_t ipZ)
{
// Reset internal data
+ AliDebug(30, Form("before event with vertex %f", ipZ));
AliFMDMultAlgorithm::PreEvent(tree, ipZ);
fCurrentVertexZ = ipZ;
fEmpty.Reset(kFALSE);
// vertex of this event
//
if (!digit) return;
- if (count < fThreshold) fEmpty(digit->Detector() - 1,
+ AliDebug(30, Form("Processing digit %s (%s)", digit->GetName(),
+ count < fThreshold ? "empty" : "hit"));
+ if (count < fThreshold) fEmpty(digit->Detector(),
digit->Ring(),
digit->Sector(),
digit->Strip()) = kTRUE;
Int_t emptyStrips = 0;
for (Int_t sector = minSector; sector < maxSector; sector++)
for (Int_t strip = minStrip; strip < maxStrip; strip++)
- if (fEmpty(sub->GetId() - 1, r->GetId(), sector, strip))
+ if (fEmpty(sub->GetId(), r->GetId(), sector, strip))
emptyStrips++;
// The total number of strips
// 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
-#include "AliFMDRing.h" // ALIFMDRING_H
-#include "AliLog.h" // ALILOG_H
+#include "AliLog.h" // ALILOG_H
+#include "AliFMDParameters.h" // ALIFMDPARAMETERS_H
+#include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
+#include "AliFMDRing.h" // ALIFMDRING_H
+#include "AliFMDCalibGain.h" // ALIFMDCALIBGAIN_H
+#include "AliFMDCalibPedestal.h" // ALIFMDCALIBPEDESTAL_H
+#include "AliFMDCalibSampleRate.h" // ALIFMDCALIBPEDESTAL_H
#include <Riostream.h>
//====================================================================
//____________________________________________________________________
AliFMDParameters::AliFMDParameters()
- : fSiDeDxMip(1.664)
+ : fSiDeDxMip(1.664),
+ fFixedPulseGain(0),
+ fEdepMip(0),
+ fZeroSuppression(0),
+ fSampleRate(0),
+ fPedestal(0),
+ fPulseGain(0),
+ fDeadMap(0)
{
// Default constructor
SetVA1MipRange();
SetPedestal();
SetPedestalWidth();
SetPedestalFactor();
+ SetThreshold();
}
+//__________________________________________________________________
+Float_t
+AliFMDParameters::GetThreshold() const
+{
+ if (!fPulseGain) return fFixedThreshold;
+ return fPulseGain->Threshold();
+}
+
+//__________________________________________________________________
+Float_t
+AliFMDParameters::GetPulseGain(UShort_t detector, Char_t ring,
+ UShort_t sector, UShort_t strip) const
+{
+ // Returns the pulser calibrated gain for strip # strip in sector #
+ // sector or ring id ring of detector # detector.
+ //
+ // For simulation, this is normally set to
+ //
+ // VA1_MIP_Range
+ // ------------------ * MIP_Energy_Loss
+ // ALTRO_channel_size
+ //
+ if (!fPulseGain) {
+ if (fFixedPulseGain <= 0)
+ fFixedPulseGain = fVA1MipRange * GetEdepMip() / fAltroChannelSize;
+ return fFixedPulseGain;
+ }
+ return fPulseGain->Value(detector, ring, sector, strip);
+}
+
+//__________________________________________________________________
+Bool_t
+AliFMDParameters::IsDead(UShort_t detector, Char_t ring,
+ UShort_t sector, UShort_t strip) const
+{
+ if (!fDeadMap) return kFALSE;
+ return fDeadMap->operator()(detector, ring, sector, strip);
+}
+
+//__________________________________________________________________
+UShort_t
+AliFMDParameters::GetZeroSuppression(UShort_t detector, Char_t ring,
+ UShort_t sector, UShort_t strip) const
+{
+ if (!fZeroSuppression) return fFixedZeroSuppression;
+ // Need to map strip to ALTRO chip.
+ return fZeroSuppression->operator()(detector, ring, sector, strip/128);
+}
+
+//__________________________________________________________________
+UShort_t
+AliFMDParameters::GetSampleRate(UShort_t ddl) const
+{
+ if (!fSampleRate) return fFixedSampleRate;
+ // Need to map sector to digitizier card.
+ return fSampleRate->Rate(ddl);
+}
+//__________________________________________________________________
+Float_t
+AliFMDParameters::GetPedestal(UShort_t detector, Char_t ring,
+ UShort_t sector, UShort_t strip) const
+{
+ if (!fPedestal) return fFixedPedestal;
+ return fPedestal->Value(detector, ring, sector, strip);
+}
+
+//__________________________________________________________________
+Float_t
+AliFMDParameters::GetPedestalWidth(UShort_t detector, Char_t ring,
+ UShort_t sector, UShort_t strip) const
+{
+ if (!fPedestal) return fFixedPedestalWidth;
+ return fPedestal->Width(detector, ring, sector, strip);
+}
+
+
//__________________________________________________________________
Float_t
AliFMDParameters::GetEdepMip() const
{
// Get energy deposited by a MIP in the silicon sensors
- AliFMDGeometry* fmd = AliFMDGeometry::Instance();
- return (fSiDeDxMip
- * fmd->GetRing('I')->GetSiThickness()
- * fmd->GetSiDensity());
+ if (fEdepMip <= 0){
+ AliFMDGeometry* fmd = AliFMDGeometry::Instance();
+ fEdepMip = (fSiDeDxMip
+ * fmd->GetRing('I')->GetSiThickness()
+ * fmd->GetSiDensity());
+ }
+ return fEdepMip;
}
//____________________________________________________________________
#ifndef ROOT_TNamed
# include <TNamed.h>
#endif
+#ifndef ROOT_TArrayI
+# include <TArrayI.h>
+#endif
+#ifndef ALIFMDUSHORTMAP_H
+# include <AliFMDUShortMap.h>
+#endif
+#ifndef ALIFMDBOOLMAP_H
+# include <AliFMDBoolMap.h>
+#endif
+typedef AliFMDUShortMap AliFMDCalibZeroSuppression;
+typedef AliFMDBoolMap AliFMDCalibDeadMap;
+class AliFMDCalibPedestal;
+class AliFMDCalibGain;
+class AliFMDCalibSampleRate;
class AliFMDParameters : public TNamed
{
public:
static AliFMDParameters* Instance();
- // Set various parameters
- void SetVA1MipRange(UShort_t r=20) { fVA1MipRange = r; }
- void SetAltroChannelSize(UShort_t s=1024) { fAltroChannelSize = s;}
- void SetChannelsPerAltro(UShort_t size=128) { fChannelsPerAltro = size; }
- void SetZeroSuppression(UShort_t s=0) { fZeroSuppression = s; }
- void SetSampleRate(UShort_t r=1) { fSampleRate = (r>2?2:r);}
- void SetPedestal(Float_t p=10) { fPedestal = p; }
- void SetPedestalWidth(Float_t w=1) { fPedestalWidth = w; }
- void SetPedestalFactor(Float_t f=3) { fPedestalFactor = f; }
+ // Set various `Fixed' parameters
+ void SetVA1MipRange(UShort_t r=20) { fVA1MipRange = r; }
+ void SetAltroChannelSize(UShort_t s=1024) { fAltroChannelSize = s;}
+ void SetChannelsPerAltro(UShort_t size=128) { fChannelsPerAltro = size; }
+ void SetPedestalFactor(Float_t f=3) { fPedestalFactor = f; }
+
+ // Set various variable parameter defaults
+ void SetZeroSuppression(UShort_t s=0) { fFixedZeroSuppression = s; }
+ void SetSampleRate(UShort_t r=1) { fFixedSampleRate = (r>2?2:r);}
+ void SetPedestal(Float_t p=10) { fFixedPedestal = p; }
+ void SetPedestalWidth(Float_t w=1) { fFixedPedestalWidth = w; }
+ void SetThreshold(Float_t t=0) { fFixedThreshold = t; }
- // Get various parameters
+ // Get `Fixed' various parameters
UShort_t GetVA1MipRange() const { return fVA1MipRange; }
UShort_t GetAltroChannelSize() const { return fAltroChannelSize; }
UShort_t GetChannelsPerAltro() const { return fChannelsPerAltro; }
- UShort_t GetZeroSuppression() const { return fZeroSuppression; }
- UShort_t GetSampleRate() const { return fSampleRate; }
Float_t GetEdepMip() const;
- Float_t GetPedestal() const { return fPedestal; }
- Float_t GetPedestalWidth() const { return fPedestalWidth; }
Float_t GetPedestalFactor() const { return fPedestalFactor; }
+ // Get variable parameters
+ Bool_t IsDead(UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip) const;
+ Float_t GetThreshold() const;
+ Float_t GetPulseGain(UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip) const;
+ Float_t GetPedestal(UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip) const;
+ Float_t GetPedestalWidth(UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip) const;
+ UShort_t GetZeroSuppression(UShort_t detector,
+ Char_t ring,
+ UShort_t sector,
+ UShort_t strip) const;
+ UShort_t GetSampleRate(UShort_t ddl) const;
+
enum {
kBaseDDL = 0x1000 // DDL offset for the FMD
};
virtual ~AliFMDParameters() {}
static AliFMDParameters* fgInstance; // Static singleton instance
- const Float_t fSiDeDxMip; // MIP dE/dx in Silicon
- UShort_t fVA1MipRange; // # MIPs the pre-amp can do
- UShort_t fAltroChannelSize; // Largest # to store in 1 ADC ch.
- UShort_t fChannelsPerAltro; // Number of pre-amp. channels/adc chan.
- UShort_t fZeroSuppression; // Threshold for zero-suppression
- UShort_t fSampleRate; // Times the ALTRO samples pre-amp.
- Float_t fPedestal; // Pedestal to subtract
- Float_t fPedestalWidth; // Width of pedestal
- Float_t fPedestalFactor; // Number of pedestal widths
+ const Float_t fSiDeDxMip; // MIP dE/dx in Silicon
+ UShort_t fVA1MipRange; // # MIPs the pre-amp can do
+ UShort_t fAltroChannelSize; // Largest # to store in 1 ADC ch.
+ UShort_t fChannelsPerAltro; // Number of pre-amp. chan/adc chan.
+ Float_t fPedestalFactor; // Number of pedestal widths
+ Float_t fFixedPedestal; // Pedestal to subtract
+ Float_t fFixedPedestalWidth; // Width of pedestal
+ UShort_t fFixedZeroSuppression; // Threshold for zero-suppression
+ UShort_t fFixedSampleRate; // Times the ALTRO samples pre-amp.
+ Float_t fFixedThreshold; //
+ mutable Float_t fFixedPulseGain; //! Gain (cached)
+ mutable Float_t fEdepMip; //! Cache of energy loss for a MIP
+
+ AliFMDCalibZeroSuppression* fZeroSuppression; // Zero suppression from CDB
+ AliFMDCalibSampleRate* fSampleRate; // Sample rate from CDB
+ AliFMDCalibPedestal* fPedestal; // Pedestals
+ AliFMDCalibGain* fPulseGain; // Pulser gain
+ AliFMDCalibDeadMap* fDeadMap; // Pulser gain
+
- ClassDef(AliFMDParameters,1)
+ ClassDef(AliFMDParameters,2)
};
#endif
fSampleRate(1)
{
// Default CTOR
- AliFMDParameters* pars = AliFMDParameters::Instance();
- fSampleRate = pars->GetSampleRate();
}
TClonesArray* array = new TClonesArray("AliFMDDigit");
fTree->Branch("FMD", &array);
+ // Get sample rate
+ AliFMDParameters* pars = AliFMDParameters::Instance();
+ fSampleRate = pars->GetSampleRate(AliFMDParameters::kBaseDDL);
+
// Use AliAltroRawStream to read the ALTRO format. No need to
// reinvent the wheel :-)
AliFMDRawStream input(fReader, fSampleRate);
// Select FMD DDL's
- fReader->Select(AliFMDParameters::kBaseDDL >> 8);
+ fReader->Select(AliFMDParameters::kBaseDDL);
Int_t oldDDL = -1;
Int_t count = 0;
fFMD(fmd)
{
AliFMDParameters* pars = AliFMDParameters::Instance();
- fSampleRate = pars->GetSampleRate();
- fThreshold = pars->GetZeroSuppression();
+ fSampleRate = pars->GetSampleRate(AliFMDParameters::kBaseDDL);
fChannelsPerAltro = pars->GetChannelsPerAltro();
}
}
digitBranch->SetAddress(&digits);
+ AliFMDParameters* pars = AliFMDParameters::Instance();
Int_t nEvents = Int_t(digitTree->GetEntries());
for (Int_t event = 0; event < nEvents; event++) {
fFMD->ResetDigits();
Char_t ring = digit->Ring();
UShort_t sector = digit->Sector();
UShort_t strip = digit->Strip();
+ fThreshold = pars->GetZeroSuppression(det, ring, sector, strip);
if (det != prevDetector) {
AliDebug(15, Form("FMD: New DDL, was %d, now %d",
AliFMDParameters::kBaseDDL + prevDetector - 1,
#include "AliFMDReconstructor.h" // ALIFMDRECONSTRUCTOR_H
#include "AliFMDRawStream.h" // ALIFMDRAWSTREAM_H
#include "AliFMDRawReader.h" // ALIFMDRAWREADER_H
-#include "AliFMDMultAlgorithm.h" // ALIFMDMULTALGORITHM_H
-#include "AliFMDMultPoisson.h" // ALIFMDMULTPOISSON_H
-#include "AliFMDMultNaiive.h" // ALIFMDMULTNAIIVE_H
+#include "AliFMDMultStrip.h" // ALIFMDMULTNAIIVE_H
#include "AliESD.h" // ALIESD_H
+#include <AliESDFMD.h> // ALIESDFMD_H
+#include <TFile.h>
//____________________________________________________________________
ClassImp(AliFMDReconstructor)
//____________________________________________________________________
AliFMDReconstructor::AliFMDReconstructor()
: AliReconstructor(),
- fPedestal(0),
- fPedestalWidth(0),
- fPedestalFactor(0)
+ fMult(0),
+ fESDObj(0)
{
- // Make a new FMD reconstructor object - default CTOR.
- AliFMDParameters* pars = AliFMDParameters::Instance();
- fPedestal = pars->GetPedestal();
- fPedestalWidth = pars->GetPedestalWidth();
- fPedestalFactor = pars->GetPedestalFactor();
-
- fAlgorithms.Add(new AliFMDMultNaiive);
- fAlgorithms.Add(new AliFMDMultPoisson);
-
- TIter next(&fAlgorithms);
- AliFMDMultAlgorithm* algorithm = 0;
- while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
- algorithm->PreRun(0);
+ // Make a new FMD reconstructor object - default CTOR.
}
//____________________________________________________________________
AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
- : AliReconstructor(),
- fPedestal(0),
- fPedestalWidth(0),
- fPedestalFactor(0)
+ : AliReconstructor(),
+ fMult(other.fMult),
+ fESDObj(other.fESDObj)
{
// Copy constructor
- fPedestal = other.fPedestal;
- fPedestalWidth = other.fPedestalWidth;
- fPedestalFactor = other.fPedestalFactor;
-
-
- fAlgorithms.Delete();
- TIter next(&(other.fAlgorithms));
- AliFMDMultAlgorithm* algorithm = 0;
- while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
- fAlgorithms.Add(algorithm);
- fAlgorithms.SetOwner(kFALSE);
}
AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
{
// Assignment operator
- fPedestal = other.fPedestal;
- fPedestalWidth = other.fPedestalWidth;
- fPedestalFactor = other.fPedestalFactor;
-
- fAlgorithms.Delete();
- TIter next(&(other.fAlgorithms));
- AliFMDMultAlgorithm* algorithm = 0;
- while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
- fAlgorithms.Add(algorithm);
- fAlgorithms.SetOwner(kFALSE);
-
+ fMult = other.fMult;
+ fESDObj = other.fESDObj;
return *this;
}
AliFMDReconstructor::~AliFMDReconstructor()
{
// Destructor
- fAlgorithms.Delete();
+ if (fMult) fMult->Delete();
+ if (fMult) delete fMult;
+ if (fESDObj) delete fESDObj;
}
//____________________________________________________________________
{
// Initialize the reconstructor
AliDebug(1, Form("Init called with runloader 0x%x", runLoader));
+ // Current vertex position
fCurrentVertex = 0;
+ // Create array of reconstructed strip multiplicities
+ fMult = new TClonesArray("AliFMDMultStrip", 51200);
+ // Create ESD output object
+ fESDObj = new AliESDFMD;
+
+ // Check that we have a run loader
if (!runLoader) {
Warning("Init", "No run loader");
return;
}
+
+ // Check if we can get the header tree
AliHeader* header = runLoader->GetHeader();
if (!header) {
Warning("Init", "No header");
return;
}
+
+ // Check if we can get a simulation header
AliGenEventHeader* eventHeader = header->GenEventHeader();
if (eventHeader) {
TArrayF vtx;
TTree* digitsTree) const
{
// Convert Raw digits to AliFMDDigit's in a tree
+ AliDebug(1, "Reading raw data into digits tree");
AliFMDRawReader rawRead(reader, digitsTree);
// rawRead.SetSampleRate(fFMD->GetSampleRate());
rawRead.Exec();
AliDebug(1, "Reconstructing from digits in a tree");
#if 0
if (fESD) {
- // const AliESDVertex* vertex = fESD->GetVertex();
- // if (vertex) {
- // AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
- // fCurrentVertex = vertex->GetZv();
- // }
+ const AliESDVertex* vertex = fESD->GetVertex();
+ if (vertex) {
+ AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
+ fCurrentVertex = vertex->GetZv();
+ }
}
#endif
TBranch *digitBranch = digitsTree->GetBranch("FMD");
TClonesArray* digits = new TClonesArray("AliFMDDigit");
digitBranch->SetAddress(&digits);
- TIter next(&fAlgorithms);
- AliFMDMultAlgorithm* algorithm = 0;
- while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
- algorithm->PreEvent(clusterTree, fCurrentVertex);
+ // TIter next(&fAlgorithms);
+ // AliFMDMultAlgorithm* algorithm = 0;
+ // while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) {
+ // AliDebug(10, Form("PreEvent called for algorithm %s",
+ // algorithm->GetName()));
+ // algorithm->PreEvent(clusterTree, fCurrentVertex);
+ // }
+ if (fMult) fMult->Clear();
+ if (fESDObj) fESDObj->Clear();
+
+ fNMult = 0;
+ fTreeR = clusterTree;
+ fTreeR->Branch("FMD", &fMult);
+
+ AliDebug(5, "Getting entry 0 from digit branch");
digitBranch->GetEntry(0);
+ AliDebug(5, "Processing digits");
ProcessDigits(digits);
- next.Reset();
- algorithm = 0;
- while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
- algorithm->PostEvent();
- clusterTree->Fill();
+ // next.Reset();
+ // algorithm = 0;
+ // while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) {
+ // AliDebug(10, Form("PostEvent called for algorithm %s",
+ // algorithm->GetName()));
+ // algorithm->PostEvent();
+ // }
+ Int_t written = clusterTree->Fill();
+ AliDebug(10, Form("Filled %d bytes into cluster tree", written));
digits->Delete();
delete digits;
}
+
//____________________________________________________________________
void
AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
AliDebug(1, Form("Got %d digits", nDigits));
for (Int_t i = 0; i < nDigits; i++) {
AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
- AliFMDGeometry* fmd = AliFMDGeometry::Instance();
- AliFMDDetector* subDetector = fmd->GetDetector(digit->Detector());
- if (!subDetector) {
- Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
- continue;
- }
-
- AliFMDRing* ring = subDetector->GetRing(digit->Ring());
- Float_t ringZ = subDetector->GetRingZ(digit->Ring());
- if (!ring) {
- Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(),
- digit->Ring());
- break;
- }
+ AliFMDParameters* param = AliFMDParameters::Instance();
+ // Check that the strip is not marked as dead
+ if (param->IsDead(digit->Detector(), digit->Ring(),
+ digit->Sector(), digit->Strip())) continue;
+
+ // digit->Print();
+ // Get eta and phi
+ Float_t eta, phi;
+ PhysicalCoordinates(digit, eta, phi);
- Float_t realZ = fCurrentVertex + ringZ;
- Float_t stripR = ((ring->GetHighR() - ring->GetLowR())
- / ring->GetNStrips() * (digit->Strip() + .5)
- + ring->GetLowR());
- Float_t theta = TMath::ATan2(stripR, realZ);
- Float_t phi = (2 * TMath::Pi() / ring->GetNSectors()
- * (digit->Sector() + .5));
- Float_t eta = -TMath::Log(TMath::Tan(theta / 2));
+ // Substract pedestal.
UShort_t counts = SubtractPedestal(digit);
- TIter next(&fAlgorithms);
- AliFMDMultAlgorithm* algorithm = 0;
- while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
- algorithm->ProcessDigit(digit, eta, phi, counts);
+ // Gain match digits.
+ Double_t edep = Adc2Energy(digit, eta, counts);
+
+ // Make rough multiplicity
+ Double_t mult = Energy2Multiplicity(digit, edep);
+
+ AliDebug(10, Form("FMD%d%c[%2d,%3d]: "
+ "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
+ digit->Detector(), digit->Ring(), digit->Sector(),
+ digit->Strip(), digit->Counts(), counts, edep, mult));
+
+ // Create a `RecPoint' on the output branch.
+ AliFMDMultStrip* m =
+ new ((*fMult)[fNMult]) AliFMDMultStrip(digit->Detector(),
+ digit->Ring(),
+ digit->Sector(),
+ digit->Strip(),
+ eta, phi,
+ edep, mult,
+ AliFMDMult::kNaiive);
+ (void)m; // Suppress warnings about unused variables.
+ fNMult++;
+
+ fESDObj->SetMultiplicity(digit->Detector(), digit->Ring(),
+ digit->Sector(), digit->Strip(), mult);
+ fESDObj->SetEta(digit->Detector(), digit->Ring(),
+ digit->Sector(), digit->Strip(), eta);
}
}
-
+
//____________________________________________________________________
UShort_t
AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
// load this to subtract a pedestal that was given in a database or
// something like that.
- Int_t counts = 0;
- Float_t ped = fPedestal + fPedestalFactor * fPedestalWidth;
+ Int_t counts = 0;
+ AliFMDParameters* param = AliFMDParameters::Instance();
+ Float_t pedM = param->GetPedestal(digit->Detector(),
+ digit->Ring(),
+ digit->Sector(),
+ digit->Strip());
+ AliDebug(10, Form("Subtracting pedestal %f from signal %d",
+ pedM, digit->Counts()));
if (digit->Count3() > 0) counts = digit->Count3();
else if (digit->Count2() > 0) counts = digit->Count2();
else counts = digit->Count1();
- counts = TMath::Max(Int_t(counts - ped), 0);
+ counts = TMath::Max(Int_t(counts - pedM), 0);
+ if (counts > 0) AliDebug(10, "Got a hit strip");
+
return UShort_t(counts);
}
//____________________________________________________________________
+Float_t
+AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit,
+ Float_t /* eta */,
+ UShort_t count) const
+{
+ // Converts number of ADC counts to energy deposited.
+ // Note, that this member function can be overloaded by derived
+ // classes to do strip-specific look-ups in databases or the like,
+ // to find the proper gain for a strip.
+ //
+ // In this simple version, we calculate the energy deposited as
+ //
+ // EnergyDeposited = cos(theta) * gain * count
+ //
+ // where
+ //
+ // Pre_amp_MIP_Range
+ // gain = ----------------- * Energy_deposited_per_MIP
+ // ADC_channel_size
+ //
+ // is constant and the same for all strips.
+
+ // Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
+ // Double_t edep = TMath::Abs(TMath::Cos(theta)) * fGain * count;
+ AliFMDParameters* param = AliFMDParameters::Instance();
+ Float_t gain = param->GetPulseGain(digit->Detector(),
+ digit->Ring(),
+ digit->Sector(),
+ digit->Strip());
+ Double_t edep = count * gain;
+ AliDebug(10, Form("Converting counts %d to energy via factor %f",
+ count, gain));
+ return edep;
+}
+
+//____________________________________________________________________
+Float_t
+AliFMDReconstructor::Energy2Multiplicity(AliFMDDigit* /* digit */,
+ Float_t edep) const
+{
+ // Converts an energy signal to number of particles.
+ // Note, that this member function can be overloaded by derived
+ // classes to do strip-specific look-ups in databases or the like,
+ // to find the proper gain for a strip.
+ //
+ // In this simple version, we calculate the multiplicity as
+ //
+ // multiplicity = Energy_deposited / Energy_deposited_per_MIP
+ //
+ // where
+ //
+ // Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
+ //
+ // is constant and the same for all strips
+ AliFMDParameters* param = AliFMDParameters::Instance();
+ Double_t edepMIP = param->GetEdepMip();
+ Float_t mult = edep / edepMIP;
+ if (edep > 0)
+ AliDebug(10, Form("Translating energy %f to multiplicity via "
+ "divider %f->%f", edep, edepMIP, mult));
+ return mult;
+}
+
+//____________________________________________________________________
+void
+AliFMDReconstructor::PhysicalCoordinates(AliFMDDigit* digit,
+ Float_t& eta,
+ Float_t& phi) const
+{
+ // Get the eta and phi of a digit
+ //
+ // Get geometry.
+ AliFMDGeometry* fmd = AliFMDGeometry::Instance();
+ AliFMDDetector* subDetector = fmd->GetDetector(digit->Detector());
+ if (!subDetector) {
+ Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
+ return;
+ }
+
+ // Get the ring - we should really use geometry->Detector2XYZ(...)
+ // here.
+ AliFMDRing* ring = subDetector->GetRing(digit->Ring());
+ Float_t ringZ = subDetector->GetRingZ(digit->Ring());
+ if (!ring) {
+ Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(),
+ digit->Ring());
+ return;
+ }
+
+ // Correct for vertex offset.
+ Float_t realZ = fCurrentVertex + ringZ;
+ Float_t stripR = ((ring->GetHighR() - ring->GetLowR())
+ / ring->GetNStrips() * (digit->Strip() + .5)
+ + ring->GetLowR());
+ Float_t theta = TMath::ATan2(stripR, realZ);
+ phi = (2 * TMath::Pi() / ring->GetNSectors()
+ * (digit->Sector() + .5));
+ eta = -TMath::Log(TMath::Tan(theta / 2));
+}
+
+
+
+//____________________________________________________________________
void
AliFMDReconstructor::FillESD(TTree* /* digitsTree */,
TTree* /* clusterTree */,
- AliESD* /* esd*/) const
+ AliESD* esd) const
{
// nothing to be done
// FIXME: The vertex may not be known when Reconstruct is executed,
// so we may have to move some of that member function here.
+ AliDebug(1, Form("Calling FillESD with two trees and one ESD"));
+ // fESDObj->Print();
+
+ if (esd) {
+ AliDebug(1, Form("Writing FMD data to ESD tree"));
+ esd->SetFMDData(fESDObj);
+ // Let's check the data in the ESD
+#if 0
+ AliESDFMD* fromEsd = esd->GetFMDData();
+ if (!fromEsd) {
+ AliWarning("No FMD object in ESD!");
+ return;
+ }
+ for (UShort_t det = 1; det <= fESDObj->MaxDetectors(); det++) {
+ for (UShort_t ir = 0; ir < fESDObj->MaxRings(); ir++) {
+ Char_t ring = (ir == 0 ? 'I' : 'O');
+ for (UShort_t sec = 0; sec < fESDObj->MaxSectors(); sec++) {
+ for (UShort_t str = 0; str < fESDObj->MaxStrips(); str++) {
+ if (fESDObj->Multiplicity(det, ring, sec, str) !=
+ fromEsd->Multiplicity(det, ring, sec, str))
+ AliWarning(Form("Mult for FMD%d%c[%2d,%3d]",det,ring,sec,str));
+ if (fESDObj->Eta(det, ring, sec, str) !=
+ fromEsd->Eta(det, ring, sec, str))
+ AliWarning(Form("Eta for FMD%d%c[%2d,%3d]", det,ring,sec,str));
+ if (fESDObj->Multiplicity(det, ring, sec, str) > 0 &&
+ fESDObj->Multiplicity(det, ring, sec, str)
+ != AliESDFMD::kInvalidMult)
+ AliInfo(Form("Mult in FMD%d%c[%2d,%3d] is %f", det,ring,sec,str,
+ fESDObj->Multiplicity(det, ring, sec, str)));
+ }
+ }
+ }
+ }
+#endif
+ }
+
+#if 0
+ static Int_t evNo = -1;
+ evNo++;
+ if (esd) evNo = esd->GetEventNumber();
+ TString fname(Form("FMD.ESD.%03d.root", evNo));
+ TFile* file = TFile::Open(fname.Data(), "RECREATE");
+ if (!file) {
+ AliError(Form("Failed to open file %s", fname.Data()));
+ return;
+ }
+ fESDObj->Write();
+ file->Close();
+#endif
+
#if 0
TClonesArray* multStrips = 0;
TClonesArray* multRegions = 0;
{
// Cannot be used. See member function with same name but with 2
// TTree arguments. Make sure you do local reconstrucion
+ AliDebug(1, Form("Calling FillESD with loader and tree"));
AliError("MayNotUse");
}
//____________________________________________________________________
{
// Cannot be used. See member function with same name but with 2
// TTree arguments. Make sure you do local reconstrucion
+ AliDebug(1, Form("Calling FillESD with loader"));
AliError("MayNotUse");
}
//____________________________________________________________________
{
// Cannot be used. See member function with same name but with 2
// TTree arguments. Make sure you do local reconstrucion
+ AliDebug(1, Form("Calling FillESD with loader and raw reader"));
AliError("MayNotUse");
}
//____________________________________________________________________
{
// Cannot be used. See member function with same name but with 2
// TTree arguments. Make sure you do local reconstrucion
+ AliDebug(1, Form("Calling FillESD with raw reader, tree, and ESD"));
AliError("MayNotUse");
}
//____________________________________________________________________
{
// Cannot be used. See member function with same name but with 2
// TTree arguments. Make sure you do local reconstrucion
+ AliDebug(1, Form("Calling FillESD with loader and ESD"));
AliError("MayNotUse");
}
//____________________________________________________________________
{
// Cannot be used. See member function with same name but with 2
// TTree arguments. Make sure you do local reconstrucion
+ AliDebug(1, Form("Calling FillESD with loader, raw reader, and ESD"));
AliError("MayNotUse");
}
#ifndef ALIRECONSTRUCTOR_H
# include <AliReconstructor.h>
#endif
-#ifndef ROOT_TObjArray
-# include <TObjArray.h>
-#endif
//____________________________________________________________________
class TTree;
class AliRawReader;
class AliRunLoader;
class AliESD;
+class AliESDFMD;
//____________________________________________________________________
class AliFMDReconstructor: public AliReconstructor
protected:
virtual void ProcessDigits(TClonesArray* digits) const;
virtual UShort_t SubtractPedestal(AliFMDDigit* digit) const;
+ virtual Float_t Adc2Energy(AliFMDDigit* digit, Float_t eta,
+ UShort_t count) const;
+ virtual Float_t Energy2Multiplicity(AliFMDDigit* digit, Float_t edep) const;
+ virtual void PhysicalCoordinates(AliFMDDigit* digit, Float_t& eta,
+ Float_t& phi) const;
- TObjArray fAlgorithms; // Array of algorithms
- Float_t fPedestal; // Pedestal to subtract
- Float_t fPedestalWidth; // Width of pedestal
- Float_t fPedestalFactor;// Number of pedestal widths
+ mutable TClonesArray* fMult; // Cache of RecPoints
+ mutable Int_t fNMult; // Number of entries in fMult
+ mutable TTree* fTreeR; // Output tree
mutable Float_t fCurrentVertex; // Z-coordinate of primary vertex
+ mutable AliESDFMD* fESDObj; // ESD output object
AliESD* fESD;
ClassDef(AliFMDReconstructor, 0) // class for the FMD reconstruction
-#ifndef ALIFMDEDEPMAP_H
-#define ALIFMDEDEPMAP_H
+#ifndef ALIFMDUSHORTMAP_H
+#define ALIFMDUSHORTMAP_H
/* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights
* reserved.
*
}
#endif
if (sector < 1 || sector > n) {
- AliWarning(Form("Step", "sector # %d out of range (0-%d)", sector-1, n-1));
+ AliWarning(Form("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);
- AliWarning(Form("Step", "Track # %5d deposits a lot of energy\n"
+ AliWarning(Form("Track # %5d deposits a lot of energy\n"
" Volume: %s\n"
" Momentum: (%7.4f,%7.4f,%7.4f)\n"
" PDG: %d (%s)\n"
kMuonCocktailCent1Single, //
kMuonCocktailPer1Single, //
kMuonCocktailPer4Single,
- kFlatFMD1,
- kFlatFMD2,
- kFlatFMD3,
+ kFMD1Flat,
+ kFMD2Flat,
+ kFMD3Flat,
+ kFMDFlat,
kEgMax
};
"kMuonCocktailPer4Single",
"kFMD1Flat",
"kFMD2Flat",
- "kFMD3Flat"
+ "kFMD3Flat",
+ "kFMDFlat"
};
//____________________________________________________________________
gGener=gener;
}
break;
- case kFlatFMD1:
+ case kFMD1Flat:
{
comment = comment.Append(" Flat in FMD1 range");
- AliGenBox* gener = AliGenBox(2000);
+ AliGenBox* gener = new AliGenBox(2000);
gener->SetPart(211);
gener->SetMomentumRange(3,4);
gener->SetPhiRange(0, 360);
gGener = gener;
}
break;
- case kFlatFMD2:
+ case kFMD2Flat:
{
comment = comment.Append(" Flat in FMD2 range");
- AliGenBox* gener = AliGenBox(2000);
+ AliGenBox* gener = new AliGenBox(2000);
gener->SetPart(211);
gener->SetMomentumRange(3,4);
gener->SetPhiRange(0, 360);
gGener = gener;
}
break;
- case kFlatFMD3:
+ case kFMD3Flat:
{
comment = comment.Append(" Flat in FMD3 range");
- AliGenBox* gener = AliGenBox(2000);
+ AliGenBox* gener = new AliGenBox(2000);
gener->SetPart(211);
gener->SetMomentumRange(3,4);
gener->SetPhiRange(0, 360);
gGener = gener;
}
break;
+ case kFMDFlat:
+ {
+ comment = comment.Append(" Flat in FMD range");
+ AliGenCocktail* gener = AliGenCocktail("FMD cocktail");
+ gener->SetPart(211);
+ gener->SetMomentumRange(3,4);
+ gener->SetPhiRange(0, 360);
+ AliGenBox* gener3 = new AliGenBox(2000);
+ gener3->SetThetaRange(155.97, 176.73);
+ gener->AddGenerator(gener3, "FMD3", .33);
+ AliGenBox* gener2 = new AliGenBox(2000);
+ gener2->SetThetaRange(2.95, 20.42);
+ gener->AddGenerator(gener2, "FMD2", .33);
+ AliGenBox* gener1 = new AliGenBox(2000);
+ gener1->SetThetaRange(0.77, 3.08);
+ gener->AddGenerator(gener1, "FMD1", .34);
+ gGener = gener;
+ }
+ break;
default: break;
}
#pragma link C++ class AliFMDBaseDigit+;
#pragma link C++ class AliFMDDigit+;
#pragma link C++ class AliFMDSDigit+;
-#pragma link C++ class AliFMDMap+;
-#pragma link C++ class AliFMDFloatMap+;
#pragma link C++ class AliFMDBoolMap+;
#pragma link C++ class AliFMDUShortMap+;
#pragma link C++ class AliFMD1+;
#pragma link C++ class AliFMDParameters+;
#pragma link C++ class AliFMDCalibPedestal+;
#pragma link C++ class AliFMDCalibGain+;
+#pragma link C++ class AliFMDCalibSampleRate+;
#else
# error Not for compilation
// #pragma link C++ class AliFMDMap<UShort_t>;
// #pragma link C++ typedef AliFMDAdcMap;
#pragma link C++ class AliFMDReconstructor+;
-#pragma link C++ class AliFMDMultAlgorithm+;
-#pragma link C++ class AliFMDMultNaiive+;
-#pragma link C++ class AliFMDMultPoisson+;
+// #pragma link C++ class AliFMDMultAlgorithm+;
+// #pragma link C++ class AliFMDMultNaiive+;
+// #pragma link C++ class AliFMDMultPoisson+;
#pragma link C++ class AliFMDMult+;
#pragma link C++ class AliFMDMultRegion+;
#pragma link C++ class AliFMDMultStrip+;
#pragma link C++ class AliFMDBaseDigitizer+;
#pragma link C++ class AliFMDDigitizer+;
#pragma link C++ class AliFMDSDigitizer+;
+#pragma link C++ class AliFMDDigitizerAlla+;
#pragma link C++ class AliFMDRawWriter+;
#pragma link C++ class AliFMDAlla+;
void
Reconstruct()
{
+ AliLog::SetModuleDebugLevel("FMD", 2);
AliReconstruction rec;
rec.SetRunLocalReconstruction("FMD");
rec.SetRunVertexFinder(kFALSE);
- // rec.SetRunTracking(kFALSE);
+ rec.SetRunTracking("");
rec.SetFillESD("FMD");
- rec.SetInput("./");
+ // rec.SetInput("./");
rec.Run();
}
sim.SetConfigFile("$(ALICE_ROOT)/FMD/Config.C");
// sim.SetMakeSDigits("FMD");
sim.SetMakeDigits("FMD");
- // sim.SetWriteRawData("FMD");
+ sim.SetWriteRawData("FMD");
// sim.SetMakeDigitsFromHits("FMD");
TStopwatch w;
w.Start();
- sim.Run(10);
+ sim.Run(1);
w.Stop();
w.Print();
}
--- /dev/null
+ Things to do in the FMD offline code
+ ------------------------------------
+
+* Split rings/cone into half rings/cones for alignment.
+* In AliFMDGeometry, get the global transformation matrices for each
+ sensitive volume (use TGeoIterator).
+* Try implement alignment objects.
+* Get calibration parameters (in AliFMDParameters) from CDB.
+* Implement Detector2XYZ and XYZ2Detector via TGeoMatrix.
# $Id$
SRCS = AliFMDDigit.cxx \
- AliFMDMap.cxx \
- AliFMDFloatMap.cxx \
AliFMDBoolMap.cxx \
AliFMDUShortMap.cxx \
AliFMDCalibPedestal.cxx \
AliFMDCalibGain.cxx \
+ AliFMDCalibSampleRate.cxx \
AliFMDParameters.cxx \
AliFMDGeometry.cxx \
AliFMDRing.cxx \
AliFMDDetector.cxx \
AliFMD1.cxx \
AliFMD2.cxx \
- AliFMD3.cxx
+ AliFMD3.cxx
+
HDRS = $(SRCS:.cxx=.h)
DHDR := FMDbaseLinkDef.h
# Un comment to use old code with seperate simulator task
SRCS = AliFMDReconstructor.cxx \
AliFMDRawStream.cxx \
AliFMDRawReader.cxx \
- AliFMDMultAlgorithm.cxx \
- AliFMDMultNaiive.cxx \
- AliFMDMultPoisson.cxx \
AliFMDMultRegion.cxx \
AliFMDMult.cxx \
AliFMDMultStrip.cxx
# Un comment to use old code with seperate simulator task
# PACKCXXFLAGS := $(CXXFLAGS) -DUSE_PRE_MOVE
+# Not used any more
+# AliFMDMultAlgorithm.cxx \
+# AliFMDMultNaiive.cxx \
+# AliFMDMultPoisson.cxx \
+
#
# EOF
#
AliFMDG3OldSimulator.cxx \
AliFMDHit.cxx \
AliFMDDigitizer.cxx \
+ AliFMDDigitizerAlla.cxx \
AliFMDEdepMap.cxx \
AliFMDRawWriter.cxx \
AliFMDAlla.cxx
--- /dev/null
+#ifndef __CINT__
+#include <AliESD.h>
+#include <TObjArray.h>
+#include <TFile.h>
+#include <TSystemDirectory.h>
+#include <TString.h>
+#include <TChain.h>
+#include <iostream>
+#endif
+
+void
+CompareESD()
+{
+ TChain* chain = new TChain("esdTree");
+ TSystemDirectory* dir = new TSystemDirectory(".", ".");
+ TList* files = dir->GetListOfFiles();
+ files->Sort();
+ TIter next(files);
+ TSystemFile* file = 0;
+ TObjArray* other = new TObjArray;
+ while ((file = static_cast<TSystemFile*>(next()))) {
+ TString fname(file->GetName());
+ if (fname.Contains("AliESDs"))
+ chain->AddFile(fname.Data());
+ if (fname.Contains("FMD.ESD.")) {
+ TFile* o = TFile::Open(fname.Data());
+ if (!o) {
+ Warning("CompareESD", "Failed to open file %s", fname.Data());
+ other->Add(0);
+ continue;
+ }
+ AliESDFMD* fmdEsd = static_cast<AliESDFMD*>(o->Get("AliESDFMD"));
+ if (!fmdEsd) {
+ Warning("CompareESD", "Failed to get FMD object from reference %s",
+ fname.Data());
+ other->Add(0);
+ continue;
+ }
+ other->Add(fmdEsd);
+ }
+ }
+ delete dir;
+
+ AliESD* esd = 0;
+ chain->SetBranchAddress("ESD", &esd);
+ Int_t n = chain->GetEntries();
+ for (Int_t i = 0; i < n; i++) {
+ Int_t read = chain->GetEntry(i);
+ if (read <= 0) break;
+ std::cout << "Entry # " << i << std::endl;
+ if (!esd) break;
+ esd->Print();
+ AliESDFMD* esdObj = esd->GetFMDData();
+ AliESDFMD* refObj = static_cast<AliESDFMD*>(other->At(i));
+ if (!esdObj) {
+ Warning("CompareESD", "no FMD object in ESD");
+ continue;
+ }
+ if (!refObj) {
+ Warning("CompareESD", "no reference FMD object");
+ continue;
+ }
+ std::cout << " Comparing ... " << std::endl;
+ for (UShort_t det = 1; det <= esdObj->MaxDetectors(); det++) {
+ std::cout << "FMD" << det << std::endl;
+ for (UShort_t ir = 0; ir < esdObj->MaxRings(); ir++) {
+ Char_t ring = (ir == 0 ? 'I' : 'O');
+ std::cout << " Ring " << ring << std::endl;
+ for (UShort_t sec = 0; sec < esdObj->MaxSectors(); sec++) {
+ std::cout << " Sector # " << sec << std::endl;
+ for (UShort_t str = 0; str < esdObj->MaxStrips(); str++) {
+ if (esdObj->Multiplicity(det, ring, sec, str) !=
+ refObj->Multiplicity(det, ring, sec, str))
+ Warning("CompareESD",
+ "Mult for FMD%d%c[%2d,%3d]",det,ring,sec,str);
+ if (esdObj->Eta(det, ring, sec, str) !=
+ refObj->Eta(det, ring, sec, str))
+ Warning("CompareESD",
+ "Eta for FMD%d%c[%2d,%3d]", det,ring,sec,str);
+ if (esdObj->Multiplicity(det, ring, sec, str) > 0)
+ Info("CompareESD", "Mult in FMD%c%d[%2d,%3d] is %f",
+ det,ring,sec,str,
+ esdObj->Multiplicity(det, ring, sec, str));
+ }
+ }
+ }
+ }
+ }
+}
+
+
+
// include path to contain the relevant directories.
//
void
-Compile(const char* script)
+Compile(const char* script, Option_t* option="g")
{
gSystem->Load("libFMDutil.so");
gSystem->SetIncludePath("-I`root-config --incdir` "
"-I${ALICE_ROOT}/include "
"-I${ALICE_ROOT}/FMD "
"-I${ALICE_ROOT}/geant3/TGeant3");
- gROOT->ProcessLine(Form(".L %s++g", script));
+ gROOT->ProcessLine(Form(".L %s+%s", script, option));
}
//____________________________________________________________________
--- /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 <AliFMDUShortMap.h>
+#include <AliFMDFloatMap.h>
+#include <AliFMDMultStrip.h>
+#include <AliFMDMultRegion.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TCanvas.h>
+
+class DrawDigitsRecs : public AliFMDInputDigits
+{
+private:
+ TH2D* fAdcVsSingle; // Histogram
+ TH2D* fAdcVsRegion; // Histogram
+ TH2D* fSingleVsRegion; // Histogram
+ AliFMDUShortMap fMap;
+ AliFMDFloatMap fEta;
+ AliFMDFloatMap fPhi;
+ AliFMDFloatMap fMult;
+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;
+ }
+ DrawDigitsRecs(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1099.5,
+ Int_t n=105, Double_t mmin=-0.5, Double_t mmax=20.5)
+ {
+ AddLoad(kRecPoints);
+ fAdcVsSingle = new TH2D("adcVsSingle", "ADC vs. Multiplicity (strip)",
+ m, amin, amax, n, mmin, mmax);
+ fAdcVsSingle->SetXTitle("ADC value");
+ fAdcVsSingle->SetYTitle("Strip Multiplicity");
+
+ mmin *= 20;
+ mmax *= 20;
+ amin *= 50;
+ amax *= 50;
+ fAdcVsRegion = new TH2D("adcVsRegion", "ADC vs Muliplcity (region)",
+ m, amin, amax, n, mmin, mmax);
+ fAdcVsRegion->SetXTitle("ADC value");
+ fAdcVsRegion->SetYTitle("Region Multiplicity");
+
+ fSingleVsRegion = new TH2D("singleVsRegion", "Single vs Region",
+ n, mmin, mmax, n, mmin, mmax);
+ fSingleVsRegion->SetXTitle("Strip Multiplicity");
+ fSingleVsRegion->SetYTitle("Region Multiplicity");
+ }
+ Bool_t Begin(Int_t ev)
+ {
+ fMap.Reset();
+ fEta.Reset();
+ fPhi.Reset();
+ return AliFMDInputDigits::Begin(ev);
+ }
+ Bool_t ProcessDigit(AliFMDDigit* digit)
+ {
+ // Cache the energy loss
+ if (!digit) return kFALSE;
+ UShort_t det = digit->Detector();
+ Char_t rng = digit->Ring();
+ UShort_t sec = digit->Sector();
+ UShort_t str = digit->Strip();
+ if (str > 511) {
+ AliWarning(Form("Bad strip number %d in digit", str));
+ return kTRUE;
+ }
+ fMap(det, rng, sec, str) = digit->Counts();
+ return kTRUE;
+ }
+ Bool_t Event()
+ {
+ if (!AliFMDInputDigits::Event()) return kFALSE;
+ Int_t nEv = fTreeR->GetEntries();
+ for (Int_t i = 0; i < nEv; i++) {
+ Int_t recoRead = fTreeR->GetEntry(i);
+ if (recoRead <= 0) continue;
+ Int_t nSingle = fArrayN->GetEntries();
+ if (nSingle <= 0) continue;
+ for (Int_t j = 0; j < nSingle; j++) {
+ AliFMDMultStrip* single =
+ static_cast<AliFMDMultStrip*>(fArrayN->At(j));
+ if (!single) continue;
+ UShort_t det = single->Detector();
+ Char_t rng = single->Ring();
+ UShort_t sec = single->Sector();
+ UShort_t str = single->Strip();
+ if (str > 511) {
+ AliWarning(Form("Bad strip number %d in single", str));
+ continue;
+ }
+ fAdcVsSingle->Fill(fMap(det, rng, sec, str), single->Particles());
+ fEta(det, rng, sec, str) = single->Eta();
+ fPhi(det, rng, sec, str) = single->Phi() * 180 / TMath::Pi();
+ fMult(det, rng, sec, str) = single->Particles();
+ }
+
+ Int_t nRegion = fArrayP->GetEntries();
+ if (nRegion <= 0) continue;
+ for (Int_t j = 0; j < nRegion; j++) {
+ AliFMDMultRegion* region =
+ static_cast<AliFMDMultRegion*>(fArrayP->At(j));
+ if (!region) continue;
+ UShort_t det = region->Detector();
+ Char_t rng = region->Ring();
+ UInt_t suma = 0;
+ Float_t summ = 0;
+
+ for (size_t sec = 0; sec < AliFMDMap::kMaxSectors; sec++) {
+ Float_t phi = fPhi(det, rng, sec, 0);
+ Bool_t ok = (phi <= region->MaxPhi() && phi >= region->MinPhi());
+ if (!ok) continue;
+ for (size_t str = 0; str < AliFMDMap::kMaxStrips; str++) {
+ Float_t eta = fEta(det, rng, sec, str);
+ Float_t sign = eta < 0 ? -1. : 1.;
+ ok = (sign * eta <= sign * region->MaxEta()
+ && sign * eta >= sign * region->MinEta());
+ if (!ok) continue;
+ suma += fMap(det, rng, sec, str);
+ summ += fMult(det, rng, sec, str);
+ }
+ }
+ fAdcVsRegion->Fill(suma, region->Particles());
+ fSingleVsRegion->Fill(summ, region->Particles());
+ }
+ }
+ return kTRUE;
+ }
+ Bool_t Finish()
+ {
+ gStyle->SetPalette(1);
+ gStyle->SetOptTitle(0);
+ gStyle->SetCanvasColor(0);
+ gStyle->SetCanvasBorderSize(0);
+ gStyle->SetPadColor(0);
+ gStyle->SetPadBorderSize(0);
+
+ new TCanvas("c1", "ADC vs. Strip Multiplicity");
+ fAdcVsSingle->SetStats(kFALSE);
+ fAdcVsSingle->Draw("COLZ");
+
+ new TCanvas("c2", "ADC vs. Region Multiplicity");
+ fAdcVsRegion->SetStats(kFALSE);
+ fAdcVsRegion->Draw("COLZ");
+
+ new TCanvas("c3", "Single vs. Region Multiplicity");
+ fSingleVsRegion->SetStats(kFALSE);
+ fSingleVsRegion->Draw("COLZ");
+
+ return kTRUE;
+ }
+
+ ClassDef(DrawDigitsRecs,0);
+
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+//
+// $Id$
+//
+// Script that contains a class to draw hits, 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 <AliFMDInput.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TMarker3DBox.h>
+#include <TGeoManager.h>
+#include <TMath.h>
+#include <TApplication.h>
+#include <TButton.h>
+#include <TParticle.h>
+#include <TCanvas.h>
+#include <TView.h>
+#include <TVirtualX.h>
+
+class DrawHitPatterns : public AliFMDInputHits
+{
+private:
+ Bool_t fWait;
+ Bool_t fPrimary;
+ Int_t fPdg;
+ TObjArray* fMarkers;
+ TObjArray* fHits;
+ TCanvas* fCanvas;
+ TPad* fPad;
+ TButton* fButton;
+ TButton* fZoom;
+ TButton* fPick;
+ Bool_t fZoomMode;
+public:
+ static DrawHitPatterns* fgInstance;
+ DrawHitPatterns(Bool_t primary=kTRUE, Int_t pdg=211)
+ : fPrimary(primary), fPdg(pdg),
+ fCanvas(0), fPad(0), fButton(0)
+ {
+ AddLoad(kKinematics);
+ AddLoad(kGeometry);
+ fMarkers = new TObjArray;
+ fHits = new TObjArray;
+ fMarkers->SetOwner(kTRUE);
+ fHits->SetOwner(kFALSE);
+ fgInstance = this;
+ }
+ void Continue() { fWait = kFALSE; }
+ void Zoom() { fZoomMode = kTRUE; }
+ void Pick() { fZoomMode = kFALSE; }
+ void ExecuteEvent(Int_t event, Int_t px, Int_t py)
+ {
+ Info("ExecuteEvent", "Event %d, at (%d,%d)", px, py);
+ static Float_t x0, y0, x1, y1;
+ static Int_t px0, py0, pxOld,pyOld;
+ static Bool_t lineDraw;
+ if (px == 0 && py == 0) return;
+ if (!fZoomMode && fPad->GetView()) {
+ fPad->GetView()->ExecuteRotateView(event, px, py);
+ return;
+ }
+ gPad->SetCursor(kCross);
+ switch (event) {
+ case kButton1Down:
+ gPad->TAttLine::Modify();
+ x0 = fPad->AbsPixeltoX(px);
+ y0 = fPad->AbsPixeltoY(py);
+ px0 = pxOld = px;
+ py0 = pyOld = py;
+ lineDraw = kFALSE;
+ return;
+ case kButton1Motion:
+ if (lineDraw)
+ gVirtualX->DrawBox(px0, py0, pxOld, pyOld, TVirtualX::kHollow);
+ pxOld = px;
+ pyOld = py;
+ lineDraw = kTRUE;
+ gVirtualX->DrawBox(px0, py0, pxOld, pyOld, TVirtualX::kHollow);
+ return;
+ case kButton1Up:
+ fPad->GetCanvas()->FeedbackMode(kFALSE);
+ if (px == px0 || py == py0) return;
+ x1 = gPad->AbsPixeltoX(px);
+ y1 = gPad->AbsPixeltoY(py);
+ if (x1 < x0) std::swap(x0, x1);
+ if (y1 < y0) std::swap(y0, y1);
+ gPad->Range(x0, y0, x1, y1);
+ gPad->Modified();
+ return;
+ }
+ }
+ Int_t DistanceToPrimitive(Int_t px, Int_t py)
+ {
+ Info("DistanceToPrimitive", "@ (%d,%d)", px, py);
+ gPad->SetCursor(kCross);
+ Float_t xmin = gPad->GetX1();
+ Float_t xmax = gPad->GetX2();
+ Float_t dx = .02 * (xmax - xmin);
+ Float_t x = gPad->AbsPixeltoX(px);
+ if (x < xmin + dx || x > xmax - dx) return 9999;
+ return (fZoomMode ? 0 : 7);
+ }
+ Bool_t Begin(Int_t event)
+ {
+ if (!fCanvas) {
+ fCanvas = new TCanvas("display", "Display", 700, 700);
+ fCanvas->SetFillColor(1);
+ fCanvas->ToggleEventStatus();
+ fPad = new TPad("view3D", "3DView", 0.0, 0.1, 1.0, 1.0, 1, 0, 0);
+ fCanvas->cd();
+ fPad->Draw();
+ }
+ if (!fButton) {
+ fCanvas->cd();
+ fButton = new TButton("Continue",
+ "DrawHitPatterns::fgInstance->Continue()",
+ 0, 0, .5, .05);
+ fButton->Draw();
+ fZoom = new TButton("Zoom",
+ "DrawHitPatterns::fgInstance->Zoom()",
+ .5, 0, .75, .05);
+ fZoom->Draw();
+ fPick = new TButton("Pick",
+ "DrawHitPatterns::fgInstance->Pick()",
+ .75, 0, 1, .05);
+ fPick->Draw();
+ }
+ Info("Begin", "Clearing canvas");
+ // fCanvas->Clear();
+ if (!fGeoManager) {
+ Warning("End", "No geometry manager");
+ return kFALSE;
+ }
+ Info("Begin", "Drawing geometry");
+ fPad->cd();
+ fGeoManager->GetTopVolume()->Draw();
+ Info("Begin", "Adjusting view");
+ Int_t irep;
+ if (fPad->GetView()) {
+ fPad->GetView()->SetView(-200, -40, 80, irep);
+ fPad->GetView()->Zoom();
+ fPad->Modified();
+ fPad->cd();
+ }
+ Info("Begin", "Drawing button");
+
+ return AliFMDInputHits::Begin(event);
+ }
+
+ Bool_t ProcessHit(AliFMDHit* hit, TParticle* p)
+ {
+ if (!hit) {
+ std::cout << "No hit" << std::endl;
+ return kFALSE;
+ }
+
+ if (!p) {
+ std::cout << "No track" << std::endl;
+ return kFALSE;
+ }
+ if (fPrimary && !p->IsPrimary()) return kTRUE;
+ if (hit->IsStop()) return kTRUE;
+ if (fPdg != 0) {
+ if (TMath::Abs(hit->Pdg()) != fPdg) return kTRUE;
+ }
+ fHits->Add(hit);
+ Float_t size = TMath::Min(TMath::Max(hit->Edep() * .1, .1), 1.);
+ Float_t pt = TMath::Sqrt(hit->Py()*hit->Py()+hit->Px()*hit->Px());
+ Float_t theta = TMath::ATan2(pt, hit->Pz());
+ Float_t phi = TMath::ATan2(hit->Py(), hit->Px());
+ TMarker3DBox* marker = new TMarker3DBox(hit->X(), hit->Y(), hit->Z(),
+ size, size, size,
+ theta, phi);
+ marker->SetLineColor((p->IsPrimary() ? 3 : 4));
+ marker->SetRefObject(hit);
+ fMarkers->Add(marker);
+ return kTRUE;
+ }
+ Bool_t End()
+ {
+ // fGeoManager->GetTopVolume()->Draw();
+ fPad->cd();
+ fMarkers->Draw();
+ AppendPad();
+ fPad->Modified(kTRUE);
+ fPad->Update();
+ fPad->cd();
+ fCanvas->Modified(kTRUE);
+ fCanvas->Update();
+ fCanvas->cd();
+
+ fWait = kTRUE;
+ while (fWait) {
+ gApplication->StartIdleing();
+ gSystem->InnerLoop();
+ gApplication->StopIdleing();
+ }
+ Info("End", "After idle loop");
+ fMarkers->Delete();
+ fHits->Clear();
+ Info("End", "After clearing caches");
+ return AliFMDInputHits::End();
+ }
+ ClassDef(DrawHitPatterns,0);
+};
+
+DrawHitPatterns* DrawHitPatterns::fgInstance = 0;
+
+//____________________________________________________________________
+//
+// EOF
+//
#include <iostream>
#include <TStyle.h>
#include <TArrayF.h>
+#include <TParticle.h>
class DrawHits : public AliFMDInputHits
{
TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
{
TArrayF bins(n+1);
+ bins[0] = min;
+ if (n <= 20) {
+ for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
+ return bins;
+ }
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);
Int_t pdg=211)
: fPdg(pdg)
{
+ AddLoad(kKinematics);
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);
+ TString name(Form("elossVsP%s", (fPdg == 0 ? "MQ" : "")));
+ TString title(Form("#Delta E/#Delta x vs. p%s",
+ (fPdg == 0 ? "/(mq^{2})" : "")));
+ fElossVsPMQ = new TH2D(name.Data(), title.Data(),
+ tkine.fN-1, tkine.fArray,
+ eloss.fN-1, eloss.fArray);
fElossVsPMQ->SetXTitle(Form("p%s", (fPdg == 0 ? "/(mq^{2})" : " [GeV]")));
fElossVsPMQ->SetYTitle("#Delta E/#Delta x [MeV/cm]");
}
- Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
+ Bool_t ProcessHit(AliFMDHit* hit, TParticle* p)
{
if (!hit) {
std::cout << "No hit" << std::endl;
return kFALSE;
}
+ if (!p) {
+ std::cout << "No track" << std::endl;
+ return kFALSE;
+ }
+ if (!p->IsPrimary()) return kTRUE;
if (hit->IsStop()) return kTRUE;
Float_t x = hit->P();
Float_t y = hit->Edep()/hit->Length();
gStyle->SetPadColor(0);
gStyle->SetPadBorderSize(0);
fElossVsPMQ->SetStats(kFALSE);
- fElossVsPMQ->Draw("COLZ");
+ fElossVsPMQ->Draw("COLZ box");
return kTRUE;
}
+ ClassDef(DrawHits,0);
};
//____________________________________________________________________
Char_t rng = hit->Ring();
UShort_t sec = hit->Sector();
UShort_t str = hit->Strip();
+ if (str > 511) {
+ AliWarning(Form("Bad strip number %d in hit", str));
+ return kTRUE;
+ }
fMap(det, rng, sec, str).fEdep += hit->Edep();
fMap(det, rng, sec, str).fN++;
return kTRUE;
Char_t rng = digit->Ring();
UShort_t sec = digit->Sector();
UShort_t str = digit->Strip();
+ if (str > 511) {
+ AliWarning(Form("Bad strip number %d in digit", str));
+ continue;
+ }
fElossVsAdc->Fill(fMap(det, rng, sec, str).fEdep, digit->Counts());
}
}
fElossVsAdc->Draw("COLZ");
return kTRUE;
}
+
+ ClassDef(DrawHitsDigits,0);
};
//____________________________________________________________________
--- /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 <AliFMDMultStrip.h>
+#include <AliFMDMultRegion.h>
+#include <AliFMDDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDEdepMap.h>
+#include <AliFMDFloatMap.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TCanvas.h>
+#include <TParticle.h>
+#include <TClonesArray.h>
+#include <TTree.h>
+#include <AliStack.h>
+#include <AliLog.h>
+
+class DrawHitsRecs : public AliFMDInputHits
+{
+private:
+ TH2D* fElossVsMult; // Histogram
+ TH2D* fHitsVsStrip; // Histogram
+ TH2D* fHitsVsRegion; // Histogram
+ AliFMDEdepMap fMap;
+ AliFMDFloatMap fEta;
+ AliFMDFloatMap fPhi;
+ AliFMDFloatMap fMult;
+ Bool_t fPrimary;
+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;
+ }
+ DrawHitsRecs(Bool_t primary=kFALSE,
+ Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
+ Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
+ {
+ fPrimary = primary;
+ AddLoad(kRecPoints);
+ if (fPrimary) AddLoad(kKinematics);
+ TArrayF eloss(MakeLogScale(n, emin, emax));
+ TArrayF mults(m+1);
+ mults[0] = mmin;
+ for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
+ fElossVsMult = new TH2D("elossVsMult",
+ "#Delta E vs. Multiplicity (single)",
+ eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
+ fElossVsMult->SetXTitle("#Delta E/#Delta x [MeV/cm]");
+ fElossVsMult->SetYTitle("Strip Multiplicity");
+
+ Double_t omin = -.5;
+ Double_t omax = 7.5;
+ Int_t o = 8;
+ fHitsVsStrip = new TH2D("hitsVsStrip",
+ "# of Hits vs. Multiplicity (strip)",
+ o, omin, omax, m, mmin, mmax);
+ fHitsVsStrip->SetXTitle("# of Hits");
+ fHitsVsStrip->SetYTitle("Strip Multiplicity");
+ }
+ 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();
+ if (str > 511) {
+ AliWarning(Form("Bad strip number %d in hit", str));
+ return kTRUE;
+ }
+ if (fPrimary) {
+ TParticle* kine = fStack->Particle(hit->Track());
+ if (!kine) return kTRUE;
+ if (!kine->IsPrimary()) return kTRUE;
+ }
+
+ 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 = fTreeR->GetEntries();
+ for (Int_t i = 0; i < nEv; i++) {
+ Int_t singleRead = fTreeR->GetEntry(i);
+ if (singleRead <= 0) continue;
+ Int_t nSingle = fArrayN->GetEntries();
+ if (nSingle <= 0) continue;
+ for (Int_t j = 0; j < nSingle; j++) {
+ AliFMDMultStrip* single =
+ static_cast<AliFMDMultStrip*>(fArrayN->At(j));
+ if (!single) continue;
+ UShort_t det = single->Detector();
+ Char_t rng = single->Ring();
+ UShort_t sec = single->Sector();
+ UShort_t str = single->Strip();
+ if (str > 511) {
+ AliWarning(Form("Bad strip number %d in single", str));
+ continue;
+ }
+ if (fMap(det, rng, sec, str).fEdep > 0)
+ fElossVsMult->Fill(fMap(det, rng, sec, str).fEdep,
+ single->Particles());
+ if (fMap(det, rng, sec, str).fN > 0)
+ fHitsVsStrip->Fill(fMap(det, rng, sec, str).fN,
+ single->Particles());
+ fEta(det, rng, sec, str) = single->Eta();
+ fPhi(det, rng, sec, str) = single->Phi() * 180 / TMath::Pi();
+ }
+ }
+ return kTRUE;
+ }
+ Bool_t Finish()
+ {
+ gStyle->SetPalette(1);
+ gStyle->SetOptTitle(0);
+ gStyle->SetCanvasColor(0);
+ gStyle->SetCanvasBorderSize(0);
+ gStyle->SetPadColor(0);
+ gStyle->SetPadBorderSize(0);
+
+ new TCanvas("c1", "Energy loss vs. Strip Multiplicity");
+ fElossVsMult->SetStats(kFALSE);
+ fElossVsMult->Draw("COLZ");
+
+ new TCanvas("c2", "# of Hits vs. Strip Multiplicity");
+ fHitsVsStrip->SetStats(kFALSE);
+ fHitsVsStrip->Draw("COLZ");
+
+ return kTRUE;
+ }
+
+ ClassDef(DrawHitsRecs,0);
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+void
+WriteTree()
+{
+ TFile* file = TFile::Open("map.root", "RECREATE");
+ TTree* tree = new TTree("T", "T");
+ AliFMDFloatMap* m = new AliFMDFloatMap(1, 1, 1, 3);
+ tree->Branch("map", "AliFMDFloatMap", &m);
+ for (int i = 0; i < 3; i++) m->operator()(1,'I',0,i) = i + 1;
+ tree->Fill();
+ file->Write();
+ file->Close();
+}
+
+void
+ReadTree()
+{
+ TFile* file = TFile::Open("map.root", "READ");
+ TTree* tree = static_cast<TTree*>(file->Get("T"));
+ AliFMDFloatMap* m = 0;
+ tree->SetBranchAddress("map", &m);
+ tree->GetEntry(0);
+ for (int i = 0; i < 3; i++) {
+ std::cout << "Map(1,'I',0," << i << "): " << m->operator()(1,'I',0,i)
+ << std::endl;
+ }
+ file->Close();
+}
+
+
+void
+WriteMap()
+{
+ TFile* file = TFile::Open("map.root", "RECREATE");
+ AliFMDFloatMap* m = new AliFMDFloatMap(1, 1, 1, 3);
+ for (int i = 0; i < 3; i++) m->operator()(1,'I',0,i) = i + 1;
+ m.Write("map");
+ file->Close();
+}
+
+ReadMap()
+{
+ TFile* file = TFile::Open("map.root", "READ");
+ AliFMDFloatMap* m = static_cast<AliFMDFloatMap*>(file->Get("map"));
+ std::cout << "Got map " << map << std::endl;
+ for (int i = 0; i < 3; i++) {
+ std::cout << "Map(1,'I',0," << i << "): " << m->operator()(1,'I',0,i)
+ << std::endl;
+ }
+ file->Close();
+}
+
+
+void
+TestMapIO()
+{
+ WriteMap();
+ ReadMap();
+ WriteTree();
+ ReadTree();
+}
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff