fUseDivided = kFALSE;
fUseAssembly = kFALSE;
fUseGeo = kTRUE;
-
+
// Initialise Hit array
HitsArray();
gAlice->GetMCApp()->AddHitList(fHits);
//
}
+//____________________________________________________________________
+void
+AliFMD::FinishEvent()
+{
+ if (fSimulator) fSimulator->EndEvent();
+}
+
+
//====================================================================
//
// Graphics and event display
}
//____________________________________________________________________
-void
+AliFMDHit*
AliFMD::AddHitByFields(Int_t track,
UShort_t detector,
Char_t ring,
detector, ring, sector, strip, track, edep, hit->Edep(),
hit->Edep() + edep);
hit->SetEdep(hit->Edep() + edep);
- return;
+ return hit;
}
}
// If hit wasn't already registered, do so know.
hit = new (a[fNhits]) AliFMDHit(fIshunt, track, detector, ring, sector,
strip, x, y, z, px, py, pz, edep, pdg, t);
fNhits++;
+ return hit;
}
//____________________________________________________________________
class TBrowser;
class AliDigitizer;
class AliFMDSimulator;
+class AliFMDHit;
//____________________________________________________________________
class AliFMD : public AliDetector
void UseDivided(Bool_t use=kTRUE) { fUseDivided = use; }
void UseAssembly(Bool_t use=kTRUE) { fUseAssembly = use; }
void UseGeo(Bool_t use=kTRUE) { fUseGeo = use; }
-
+ void UseDetailed(Bool_t use=kTRUE) { fDetailed = use; }
// GEometry ANd Tracking (GEANT :-)
virtual void CreateGeometry();
virtual void CreateMaterials();
virtual void Init();
virtual void StepManager() = 0;
-
+ virtual void FinishEvent();
+
// Graphics and event display
virtual void BuildGeometry();
virtual void DrawDetector();
virtual TClonesArray* SDigits() { return fSDigits; }
virtual void ResetSDigits();
virtual void AddHit(Int_t track, Int_t *vol, Float_t *hits);
- virtual void AddHitByFields(Int_t track,
+ virtual AliFMDHit* AddHitByFields(Int_t track,
UShort_t detector,
Char_t ring,
UShort_t sector,
Char_t id = r->GetId();
Double_t siThick = r->GetSiThickness();
// const Int_t nv = r->GetNVerticies();
- TVector2* a = r->GetVertex(5);
+ //TVector2* a = r->GetVertex(5);
TVector2* b = r->GetVertex(3);
- TVector2* c = r->GetVertex(4);
+ //TVector2* c = r->GetVertex(4);
Double_t theta = r->GetTheta();
- Double_t off = (TMath::Tan(TMath::Pi() * theta / 180)
- * r->GetBondingWidth());
+ //Double_t off = (TMath::Tan(TMath::Pi() * theta / 180)
+ // * r->GetBondingWidth());
Double_t rmax = b->Mod();
Double_t rmin = r->GetLowR();
Double_t pcbThick = r->GetPrintboardThickness();
Double_t copperThick = r->GetCopperThickness(); // .01;
Double_t chipThick = r->GetChipThickness(); // .01;
- Double_t modSpace = r->GetModuleSpacing();
- Double_t legr = r->GetLegRadius();
- Double_t legl = r->GetLegLength();
- Double_t legoff = r->GetLegOffset();
+ //Double_t modSpace = r->GetModuleSpacing();
+ //Double_t legr = r->GetLegRadius();
+ //Double_t legl = r->GetLegLength();
+ //Double_t legoff = r->GetLegOffset();
Int_t ns = r->GetNStrips();
+ Double_t space = r->GetSpacing();
Int_t nsec = Int_t(360 / theta);
- Double_t stripoff = a->Mod();
- Double_t dstrip = (rmax - stripoff) / ns;
+ //Double_t stripoff = a->Mod();
+ //Double_t dstrip = (rmax - stripoff) / ns;
Double_t par[10];
TString name;
TString name2;
Int_t siId = fFMD->GetIdtmed()->At(kSiId);
Int_t airId = fFMD->GetIdtmed()->At(kAirId);
Int_t pcbId = fFMD->GetIdtmed()->At(kPcbId);
- Int_t plaId = fFMD->GetIdtmed()->At(kPlasticId);
+ //Int_t plaId = fFMD->GetIdtmed()->At(kPlasticId);
Int_t copId = fFMD->GetIdtmed()->At(kCopperId);
Int_t chiId = fFMD->GetIdtmed()->At(kSiChipId);
name2 = name;
name = Form(fgkActiveName, id);
Double_t z = - ringWidth / 2 + siThick / 2;
- mc->Gsvolu(name.Data(), "TUBE", (fDetailed ? airId : siId), par, 3);
+ mc->Gsvolu(name.Data(), "TUBE", siId, par, 3);
mc->Gspos(name.Data(), 1, name2.Data(), 0, 0, z, 0);
Int_t sid = mc->VolId(name.Data());
// Shape of Printed circuit Board
Double_t boardThick = (pcbThick + copperThick + chipThick);
- par[0] = rmin - .1;
+ par[0] = rmin + .1;
par[1] = rmax - .1;
par[2] = boardThick / 2;
name2 = Form(fgkRingName, id);
name = Form(fgkPCBName, id, 'B');
- z += siThick / 2 + boardThick / 2;
+ z += siThick / 2 + space + boardThick / 2;
mc->Gsvolu(name.Data(), "TUBE", pcbId, par, 3);
mc->Gspos(name.Data(), 1, name2.Data(), 0, 0, z, 0);
mc->Gspos(name.Data(), 2, name2.Data(), 0, 0, z + boardThick, 0);
mc->Gspos(name.Data(), 1, name2.Data(), 0, 0, z, 0);
// Copper
par[2] = copperThick / 2;
- name2 = name;
name = Form("F%cCO", id);
z += pcbThick / 2 + copperThick / 2;
mc->Gsvolu(name.Data(), "TUBE", copId, par, 3);
mc->Gspos(name.Data(), 1, name2.Data(), 0, 0, z, 0);
// Chip
par[2] = chipThick / 2;
- name2 = name;
name = Form("F%cCH", id);
- z = boardThick / 2 - chipThick / 2;
+ z += copperThick / 2 + chipThick / 2;
mc->Gsvolu(name.Data(), "TUBE", chiId, par, 3);
mc->Gspos(name.Data(), 1, name2.Data(), 0, 0, z, 0);
Int_t airId = fFMD->GetIdtmed()->At(kAirId);
Int_t pcbId = fFMD->GetIdtmed()->At(kPcbId);
Int_t plaId = fFMD->GetIdtmed()->At(kPlasticId);
- Int_t copId = fFMD->GetIdtmed()->At(kCopperId);
- Int_t chiId = fFMD->GetIdtmed()->At(kSiChipId);
+ // Int_t copId = fFMD->GetIdtmed()->At(kCopperId);
+ // Int_t chiId = fFMD->GetIdtmed()->At(kSiChipId);
Double_t ringWidth = r->GetRingDepth();
Double_t x = 0;
Char_t id = r->GetId();
Double_t siThick = r->GetSiThickness();
// const Int_t nv = r->GetNVerticies();
- TVector2* a = r->GetVertex(5);
+ //TVector2* a = r->GetVertex(5);
TVector2* b = r->GetVertex(3);
- TVector2* c = r->GetVertex(4);
+ //TVector2* c = r->GetVertex(4);
Double_t theta = r->GetTheta();
- Double_t off = (TMath::Tan(TMath::Pi() * theta / 180)
- * r->GetBondingWidth());
+ //Double_t off = (TMath::Tan(TMath::Pi() * theta / 180)
+ // * r->GetBondingWidth());
Double_t rmax = b->Mod();
Double_t rmin = r->GetLowR();
Double_t pcbThick = r->GetPrintboardThickness();
Double_t copperThick = r->GetCopperThickness(); // .01;
Double_t chipThick = r->GetChipThickness(); // .01;
- Double_t modSpace = r->GetModuleSpacing();
- Double_t legr = r->GetLegRadius();
- Double_t legl = r->GetLegLength();
- Double_t legoff = r->GetLegOffset();
+ //Double_t modSpace = r->GetModuleSpacing();
+ //Double_t legr = r->GetLegRadius();
+ //Double_t legl = r->GetLegLength();
+ //Double_t legoff = r->GetLegOffset();
Int_t ns = r->GetNStrips();
Int_t nsec = Int_t(360 / theta);
- Double_t stripoff = a->Mod();
- Double_t dstrip = (rmax - stripoff) / ns;
- Double_t par[10];
+ Double_t space = r->GetSpacing();
+ //Double_t stripoff = a->Mod();
+ //Double_t dstrip = (rmax - stripoff) / ns;
TString name;
TString name2;
- TVirtualMC* mc = TVirtualMC::GetMC();
- Int_t siId = fFMD->GetIdtmed()->At(kSiId);
- Int_t airId = fFMD->GetIdtmed()->At(kAirId);
- Int_t pcbId = fFMD->GetIdtmed()->At(kPcbId);
- Int_t plaId = fFMD->GetIdtmed()->At(kPlasticId);
- Int_t copId = fFMD->GetIdtmed()->At(kCopperId);
- Int_t chiId = fFMD->GetIdtmed()->At(kSiChipId);
-
Double_t ringWidth = (siThick + 2 * (pcbThick + copperThick + chipThick));
// Virtual volume shape to divide - This volume is only defined if
// the geometry is set to be detailed.
// Shape of Printed circuit Board
Double_t boardThick = (pcbThick + copperThick + chipThick);
- TGeoShape* boardShape = new TGeoTube(rmin + .1, rmax - .1, boardThick/ 2);
+ TGeoShape* boardShape = new TGeoTube(rmin+.1, rmax-.1, boardThick/ 2);
name = Form(fgkPCBName, id, 'B');
TGeoVolume* boardVolume = new TGeoVolume(name.Data(), boardShape, fAir);
- z += siThick / 2 + boardThick / 2;
+ z += siThick / 2 + space + boardThick / 2;
ringVolume->AddNode(boardVolume, 0, new TGeoTranslation(0, 0, z));
ringVolume->AddNode(boardVolume, 1, new TGeoTranslation(0,0,z+boardThick));
// PCB
- TGeoShape* pcbShape = new TGeoTube(rmin+.1, rmax-.1, pcbThick / 2);
+ TGeoShape* pcbShape = new TGeoTube(rmin+.1,rmax-.1, pcbThick / 2);
name = Form("F%cPC", id);
z = -boardThick / 2 + pcbThick / 2;
TGeoVolume* pcbVolume = new TGeoVolume(name.Data(), pcbShape, fPCB);
// Copper
TGeoShape* cuShape = new TGeoTube(rmin+.1, rmax-.1, copperThick / 2);
name = Form("F%cCO", id);
- z += -pcbThick / 2 + copperThick / 2;
+ z += pcbThick / 2 + copperThick / 2;
TGeoVolume* cuVolume = new TGeoVolume(name.Data(), cuShape, fCopper);
boardVolume->AddNode(cuVolume, 0, new TGeoTranslation(0, 0, z));
// Chip
TGeoShape* chipShape = new TGeoTube(rmin+.1, rmax-.1, chipThick / 2);
name = Form("F%cCH", id);
- z = -copperThick / 2 + chipThick / 2;
+ z += copperThick / 2 + chipThick / 2;
TGeoVolume* chipVolume = new TGeoVolume(name.Data(), chipShape, fChip);
boardVolume->AddNode(chipVolume, 0, new TGeoTranslation(0, 0, z));
if (fDetailed) {
TGeoTubeSeg* activeShape =
new TGeoTubeSeg(rmin, rmax, siThick/2, - theta, theta);
- activeVolume = new TGeoVolume(Form(fgkActiveName, id), activeShape, fSi);
- TGeoVolume* sectorVolume = activeVolume->Divide(Form(fgkSectorName, id),
- 2, 2, -theta, 0, 0, "N");
+ activeVolume = new TGeoVolume(Form(fgkActiveName, id),activeShape,fSi);
+ TGeoVolume* sectorVolume = activeVolume->Divide(Form(fgkSectorName,id),
+ 2, 2, -theta,0,0,"N");
TGeoVolume* stripVolume = sectorVolume->Divide(Form(fgkStripName, id),
1, ns, stripoff, dstrip,
0, "SX");
TGeoMatrix* matrix = 0;
// Back container volume
- Double_t contThick = siThick + pcbThick + legl;
+ Double_t contThick = siThick + pcbThick + legl + space;
TGeoTubeSeg* backShape = new TGeoTubeSeg(rmin, rmax, contThick/2,
- theta, theta);
TGeoVolume* backVolume = new TGeoVolume(Form(fgkBackVName, id),
// frontVolume->VisibleDaughters(kTRUE);
// Ring mother volume
- TGeoTube* ringShape = new TGeoTube(rmin, rmax, contThick / 2);
- TGeoVolume* ringVolume = new TGeoVolume(Form(fgkRingName, id), ringShape,
- fAir);
+ TGeoTube* ringShape = new TGeoTube(rmin, rmax, contThick / 2);
+ TGeoVolume* ringVolume = new TGeoVolume(Form(fgkRingName,id),
+ ringShape,fAir);
Int_t nmod = r->GetNModules();
AliDebug(10, Form("making %d modules in ring %c", nmod, id));
for (Int_t i = 0; i < nmod; i++) {
Bool_t isFront = (i % 2 == 0);
TGeoVolume* vol = (isFront ? frontVolume : backVolume);
- TGeoRotation* rot = new TGeoRotation(Form("FMD Ring %c rotation %d",id,i));
+ TGeoRotation* rot =new TGeoRotation(Form("FMD Ring %c rotation %d",id,i));
rot->RotateZ((i + .5) * 2 * theta);
Double_t z = (isFront ? 0 : modSpace) / 2;
matrix = new TGeoCombiTrans(Form("FMD Ring %c transform %d", id, i),
#include "AliFMDHit.h" // ALIFMDHIT_H
#include "AliLog.h" // ALILOG_H
#include "Riostream.h" // ROOT_Riostream
+#include <TDatabasePDG.h>
+#include <TMath.h>
+#include <TString.h>
//____________________________________________________________________
ClassImp(AliFMDHit)
fZ = z;
}
+//____________________________________________________________________
+Float_t
+AliFMDHit::P() const
+{
+ // Get the momentum of the particle of the particle that made this
+ // hit.
+ return TMath::Sqrt(fPx * fPx + fPy * fPy + fPz * fPz);
+}
+
+//____________________________________________________________________
+Float_t
+AliFMDHit::M() const
+{
+ // Get the mass of the particle that made this hit.
+ TDatabasePDG* pdgDB = TDatabasePDG::Instance();
+ TParticlePDG* pdg = pdgDB->GetParticle(fPdg);
+ return (pdg ? pdg->Mass() : -1);
+}
+
+//____________________________________________________________________
+Float_t
+AliFMDHit::Q() const
+{
+ // Get the charge of the particle that made this hit.
+ TDatabasePDG* pdgDB = TDatabasePDG::Instance();
+ TParticlePDG* pdg = pdgDB->GetParticle(fPdg);
+ return (pdg ? pdg->Charge() : 0);
+}
+
+
//____________________________________________________________________
void
-AliFMDHit::Print(Option_t* /* option */) const
+AliFMDHit::Print(Option_t* option) const
{
// Print Hit to standard out
cout << "AliFMDHit: FMD"
<< setw(3) << fSector << ","
<< setw(3) << fStrip << "] = "
<< fEdep << endl;
+ TString opt(option);
+ if (opt.Contains("D", TString::kIgnoreCase)) {
+ TDatabasePDG* pdgDB = TDatabasePDG::Instance();
+ TParticlePDG* pdg = pdgDB->GetParticle(fPdg);
+ cout << "\tPDG:\t" << fPdg << " " << (pdg ? pdg->GetName() : "?") << "\n"
+ << "\tP:\t(" << fPx << "," << fPy << "," << fPz << ") "<<P() << "\n"
+ << "\tX:\t" << fX << "," << fY << "," << fZ << "\n"
+ << "\tTrack #:\t" << fTrack << std::endl;
+ }
}
//____________________________________________________________________
Float_t Px() const { return fPx; }
Float_t Py() const { return fPy; }
Float_t Pz() const { return fPz; }
+ Float_t P() const;
+ Float_t M() const;
+ Float_t Q() const;
Int_t Pdg() const { return fPdg; }
Float_t Time() const { return fTime; }
void Print(Option_t* opt="") const;
#include "AliFMD2.h" // ALIFMD2_H
#include "AliFMD3.h" // ALIFMD3_H
#include "AliFMD.h" // ALIFMD_H
+#include "AliFMDHit.h" // ALIFMDHIT_H
#include <AliRun.h> // ALIRUN_H
#include <AliMC.h> // ALIMC_H
#include <AliMagF.h> // ALIMAGF_H
#include <TVirtualMC.h> // ROOT_TVirtualMC
#include <TArrayD.h> // ROOT_TArrayD
+
//====================================================================
ClassImp(AliFMDSimulator)
#if 0
fOuterId(-1),
fActiveId(4),
fUseDivided(kFALSE),
- fUseAssembly(kTRUE)
+ fUseAssembly(kTRUE),
+ fBad(0)
{
// Default constructor
}
fOuterId(-1),
fActiveId(4),
fUseDivided(kFALSE),
- fUseAssembly(kTRUE)
+ fUseAssembly(kTRUE),
+ fBad(0)
{
// Normal constructor
//
// fmd Pointer to AliFMD object
// detailed Whether to make a detailed simulation or not
//
+ fBad = new TClonesArray("AliFMDHit");
}
Int_t det; mc->CurrentVolOffID(fDetectorOff, det); detector = det;
AliFMDGeometry* fmd = AliFMDGeometry::Instance();
- Double_t rz = fmd->GetDetector(detector)->GetRingZ(ring);
+ //Double_t rz = fmd->GetDetector(detector)->GetRingZ(ring);
Int_t n = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors();
#if 0
if (rz < 0) {
sector--;
// Get track position
mc->TrackPosition(v);
- AliDebug(40, Form("<2> Inside an active FMD volume FMD%d%c[%2d,%3d] %s",
+ AliDebug(15, Form("<2> Inside an active FMD volume FMD%d%c[%2d,%3d] %s",
detector, ring, sector, strip, mc->CurrentVolPath()));
return kTRUE;
Int_t copy;
Int_t vol = mc->CurrentVolID(copy);
if (!IsActive(vol)) {
- AliDebug(50, Form("Not an FMD volume %d '%s' (%d or %d)",
- vol, mc->CurrentVolName(), fInnerId, fOuterId));
+ AliDebug(50, Form("Not an FMD volume %d '%s'",vol,mc->CurrentVolName()));
return;
}
TLorentzVector v;
// our parameters.
if (entering) {
AliDebug(15, Form("Track # %8d entering active FMD volume %s: "
- "Edep=%f", trackno, mc->CurrentVolPath(), edep));
+ "Edep=%f (%f,%f,%f)", trackno, mc->CurrentVolPath(),
+ edep, v.X(), v.Y(), v.Z()));
fCurrentP = p;
fCurrentV = v;
fCurrentDeltaE = edep;
- fCurrentPdg = mc->IdFromPDG(pdg);
+ fCurrentPdg = pdg; // mc->IdFromPDG(pdg);
}
// If the track is inside, then update the energy deposition
if (inside && fCurrentDeltaE >= 0) {
fCurrentDeltaE += edep;
AliDebug(15, Form("Track # %8d inside active FMD volume %s: Edep=%f, "
- "Accumulated Edep=%f", trackno, mc->CurrentVolPath(),
- edep, fCurrentDeltaE));
+ "Accumulated Edep=%f (%f,%f,%f)", trackno,
+ mc->CurrentVolPath(), edep, fCurrentDeltaE,
+ v.X(), v.Y(), v.Z()));
}
// The track exits the volume, or it disappeared in the volume, or
// the track is stopped because it no longer fulfills the cuts
// defined, then we create a hit.
- if (out && fCurrentDeltaE >= 0) {
- fCurrentDeltaE += edep;
- AliDebug(15, Form("Track # %8d exiting active FMD volume %s: Edep=%f, "
- "Accumulated Edep=%f", trackno, mc->CurrentVolPath(),
- edep, fCurrentDeltaE));
- fFMD->AddHitByFields(trackno, detector, ring, sector, strip,
- fCurrentV.X(), fCurrentV.Y(), fCurrentV.Z(),
- fCurrentP.X(), fCurrentP.Y(), fCurrentP.Z(),
- fCurrentDeltaE, fCurrentPdg, fCurrentV.T());
+ if (out) {
+ if (fCurrentDeltaE >= 0) {
+ fCurrentDeltaE += edep;
+ AliDebug(15, Form("Track # %8d exiting active FMD volume %s: Edep=%g, "
+ "Accumulated Edep=%g (%f,%f,%f)", trackno,
+ mc->CurrentVolPath(), edep, fCurrentDeltaE,
+ v.X(), v.Y(), v.Z()));
+ AliFMDHit* h =
+ fFMD->AddHitByFields(trackno, detector, ring, sector, strip,
+ fCurrentV.X(), fCurrentV.Y(), fCurrentV.Z(),
+ fCurrentP.X(), fCurrentP.Y(), fCurrentP.Z(),
+ fCurrentDeltaE, fCurrentPdg, fCurrentV.T());
+ // Add a copy
+ if (isBad && fBad) {
+ new ((*fBad)[fBad->GetEntries()]) AliFMDHit(*h);
+ }
+ }
fCurrentDeltaE = -1;
}
}
+//____________________________________________________________________
+void
+AliFMDSimulator::EndEvent()
+{
+ if (fBad && fBad->GetEntries() > 0) {
+ Warning("EndEvent", "got %d 'bad' hits", fBad->GetEntries());
+ TIter next(fBad);
+ AliFMDHit* hit;
+ while ((hit = static_cast<AliFMDHit*>(next())))
+ hit->Print("D");
+ fBad->Clear();
+ }
+}
//____________________________________________________________________
class AliFMD1;
class AliFMD2;
class AliFMD3;
+class TObjArray;
/** Simulation of the FMD.
This class builds the geometry, and processes hits in the FMD */
virtual void DefineGeometry() = 0;
/** Deal with a hit in the FMD */
virtual void Exec(Option_t* option="");
+ virtual void EndEvent();
virtual void UseDivided(Bool_t use=kTRUE) { fUseDivided = use; }
virtual void UseAssembly(Bool_t use=kTRUE) { fUseAssembly = use; }
protected:
Int_t fModuleOff; // Module offset in volume tree
Int_t fRingOff; // Ring offset in the volume tree
Int_t fDetectorOff; // Detector offfset in the volume tree
+ TObjArray* fBad; //! List of bad hits
ClassDef(AliFMDSimulator,0) // Simulation class for the FMD
};
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff