X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=FMD%2FAliFMDSimulator.cxx;h=2cacd68d07040303fb1c7d50c5c866da676f1ecc;hp=d70ac4591c6a9c064ce99b9e7db664eb8e3ad23a;hb=f4ee3d577a361fdf930eaac0ce55d7ed8b8d6444;hpb=69b696b9972ea1724c192ec1d88dae1f431f6ef9 diff --git a/FMD/AliFMDSimulator.cxx b/FMD/AliFMDSimulator.cxx index d70ac4591c6..2cacd68d070 100644 --- a/FMD/AliFMDSimulator.cxx +++ b/FMD/AliFMDSimulator.cxx @@ -37,9 +37,13 @@ // | | // +--------------------+ +-------------------+ // | AliFMDGeoSimulator | | AliFMDG3Simulator | -// +--------------------+ +---------+---------+ +// +--------------------+ +-------------------+ +// ^ +// | +// +----------------------+ +// | AliFMDG3OldSimulator | +// +----------------------+ // -// // * AliFMD // This defines the interface for the various parts of AliROOT that // uses the FMD, like AliFMDSimulator, AliFMDDigitizer, @@ -64,6 +68,10 @@ // in the corners should be cut away at run time (but currently // isn't). // +// * AliFMDG3OldSimulator +// This is a concrete implementation of AliFMDSimulator. It +// approximates the of the rings as segmented disks. +// #include "AliFMDSimulator.h" // ALIFMDSIMULATOR_H #include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H #include "AliFMDDetector.h" // ALIFMDDETECTOR_H @@ -72,6 +80,7 @@ #include "AliFMD2.h" // ALIFMD2_H #include "AliFMD3.h" // ALIFMD3_H #include "AliFMD.h" // ALIFMD_H +#include "AliFMDHit.h" // ALIFMDHIT_H #include // ALIRUN_H #include // ALIMC_H #include // ALIMAGF_H @@ -93,6 +102,7 @@ #include // ROOT_TVirtualMC #include // ROOT_TArrayD + //==================================================================== ClassImp(AliFMDSimulator) #if 0 @@ -123,19 +133,23 @@ const Char_t* AliFMDSimulator::fgkFlangeName = "F3SF"; AliFMDSimulator::AliFMDSimulator() : fFMD(0), fDetailed(kFALSE), - fInnerId(-1), - fOuterId(-1) + fActiveId(4), + fUseDivided(kFALSE), + fUseAssembly(kTRUE), + fBad(0) { // Default constructor } //____________________________________________________________________ AliFMDSimulator::AliFMDSimulator(AliFMD* fmd, Bool_t detailed) - : TTask("FMDsimulator", "Forward Multiplicity Detector Simulator"), + : TTask("FMDSimulator", "Forward Multiplicity Detector Simulator"), fFMD(fmd), fDetailed(detailed), - fInnerId(-1), - fOuterId(-1) + fActiveId(4), + fUseDivided(kFALSE), + fUseAssembly(kTRUE), + fBad(0) { // Normal constructor // @@ -144,6 +158,7 @@ AliFMDSimulator::AliFMDSimulator(AliFMD* fmd, Bool_t detailed) // fmd Pointer to AliFMD object // detailed Whether to make a detailed simulation or not // + fBad = new TClonesArray("AliFMDHit"); } @@ -159,7 +174,7 @@ AliFMDSimulator::DefineMaterials() // FMD Si$ Silicon (active medium in sensors) // FMD C$ Carbon fibre (support cone for FMD3 and vacuum pipe) // FMD Al$ Aluminium (honeycomb support plates) - // FMD PCB$ Printed Circuit Board (FEE board with VA1_ALICE) + // FMD PCB$ Printed Circuit Board (FEE board with VA1_3) // FMD Chip$ Electronics chips (currently not used) // FMD Air$ Air (Air in the FMD) // FMD Plastic$ Plastic (Support legs for the hybrid cards) @@ -168,6 +183,10 @@ AliFMDSimulator::DefineMaterials() // singleton. These pointers are later used when setting up the // geometry AliDebug(10, "\tCreating materials"); + AliDebug(1, Form("\tGeometry options: %s, %s, %s", + (fDetailed ? "detailed" : "coarse"), + (fUseDivided ? "divided into strips" : "one volume"), + (fUseAssembly ? "within assemblies" : "in real volumes"))); // Get pointer to geometry singleton object. AliFMDGeometry* geometry = AliFMDGeometry::Instance(); geometry->Init(); @@ -196,9 +215,9 @@ AliFMDSimulator::DefineMaterials() precision = .001; minStepSize = .001; id = kSiId; - fFMD->AliMaterial(id, "FMD Si$", + fFMD->AliMaterial(id, "Si$", a, z, density, radiationLength, absorbtionLength); - fFMD->AliMedium(kSiId, "FMD Si$", + fFMD->AliMedium(kSiId, "Si$", id,1,fieldType,maxField,maxBending, maxStepSize,maxEnergyLoss,precision,minStepSize); @@ -213,9 +232,9 @@ AliFMDSimulator::DefineMaterials() precision = .003; minStepSize = .003; id = kCarbonId; - fFMD->AliMaterial(id, "FMD Carbon$", + fFMD->AliMaterial(id, "Carbon$", a, z, density, radiationLength, absorbtionLength); - fFMD->AliMedium(kCarbonId, "FMD Carbon$", + fFMD->AliMedium(kCarbonId, "Carbon$", id,0,fieldType,maxField,maxBending, maxStepSize,maxEnergyLoss,precision,minStepSize); @@ -225,34 +244,46 @@ AliFMDSimulator::DefineMaterials() density = 2.7; radiationLength = 8.9; id = kAlId; - fFMD->AliMaterial(id, "FMD Aluminum$", + fFMD->AliMaterial(id, "Aluminum$", a, z, density, radiationLength, absorbtionLength); - fFMD->AliMedium(kAlId, "FMD Aluminum$", + fFMD->AliMedium(kAlId, "Aluminum$", id, 0, fieldType, maxField, maxBending, maxStepSize, maxEnergyLoss, precision, minStepSize); + // Copper + a = 63.546; + z = 29; + density = 8.96; + radiationLength = 1.43; + id = kCopperId; + fFMD->AliMaterial(id, "Copper$", + a, z, density, radiationLength, absorbtionLength); + fFMD->AliMedium(kCopperId, "Copper$", + id, 0, fieldType, maxField, maxBending, + maxStepSize, maxEnergyLoss, precision, minStepSize); + + // Silicon chip { Float_t as[] = { 12.0107, 14.0067, 15.9994, - 1.00794, 28.0855, 107.8682 }; + 1.00794, 28.0855, 107.8682 }; Float_t zs[] = { 6., 7., 8., 1., 14., 47. }; Float_t ws[] = { 0.039730642, 0.001396798, 0.01169634, 0.004367771, 0.844665, 0.09814344903 }; - density = 2.36436; + density = 2.36436; maxBending = 10; maxStepSize = .01; precision = .003; minStepSize = .003; - id = kSiChipId; - fFMD->AliMixture(id, "FMD Si Chip$", as, zs, density, 6, ws); - fFMD->AliMedium(kSiChipId, "FMD Si Chip$", + id = kSiChipId; + fFMD->AliMixture(id, "Si Chip$", as, zs, density, 6, ws); + fFMD->AliMedium(kSiChipId, "Si Chip$", id, 0, fieldType, maxField, maxBending, maxStepSize, maxEnergyLoss, precision, minStepSize); } -#if 0 // Kaption { Float_t as[] = { 1.00794, 12.0107, 14.010, 15.9994}; @@ -263,13 +294,12 @@ AliFMDSimulator::DefineMaterials() maxStepSize = .001; precision = .001; minStepSize = .001; - id = KaptionId; - fFMD->AliMixture(id, "FMD Kaption$", as, zs, density, 4, ws); - fFMD->AliMedium(kAlId, "FMD Kaption$", + id = kKaptonId; + fFMD->AliMixture(id, "Kaption$", as, zs, density, 4, ws); + fFMD->AliMedium(kKaptonId, "Kaption$", id,0,fieldType,maxField,maxBending, maxStepSize,maxEnergyLoss,precision,minStepSize); } -#endif // Air { @@ -282,8 +312,8 @@ AliFMDSimulator::DefineMaterials() precision = .001; minStepSize = .001; id = kAirId; - fFMD->AliMixture(id, "FMD Air$", as, zs, density, 4, ws); - fFMD->AliMedium(kAirId, "FMD Air$", + fFMD->AliMixture(id, "Air$", as, zs, density, 4, ws); + fFMD->AliMedium(kAirId, "Air$", id,0,fieldType,maxField,maxBending, maxStepSize,maxEnergyLoss,precision,minStepSize); } @@ -308,8 +338,8 @@ AliFMDSimulator::DefineMaterials() precision = .001; minStepSize = .001; id = kPcbId; - fFMD->AliMixture(id, "FMD PCB$", as, zs, density, 14, ws); - fFMD->AliMedium(kPcbId, "FMD PCB$", + fFMD->AliMixture(id, "PCB$", as, zs, density, 14, ws); + fFMD->AliMedium(kPcbId, "PCB$", id,0,fieldType,maxField,maxBending, maxStepSize,maxEnergyLoss,precision,minStepSize); } @@ -325,13 +355,107 @@ AliFMDSimulator::DefineMaterials() precision = .003; minStepSize = .003; id = kPlasticId; - fFMD->AliMixture(id, "FMD Plastic$", as, zs, density, -2, ws); - fFMD->AliMedium(kPlasticId, "FMD Plastic$", + fFMD->AliMixture(id, "Plastic$", as, zs, density, -2, ws); + fFMD->AliMedium(kPlasticId, "Plastic$", id,0,fieldType,maxField,maxBending, maxStepSize,maxEnergyLoss,precision,minStepSize); } } +//____________________________________________________________________ +Bool_t +AliFMDSimulator::IsActive(Int_t volId) const +{ + for (Int_t i = 0; i < fActiveId.fN; i++) + if (volId == fActiveId[i]) return kTRUE; + return kFALSE; +} + +//____________________________________________________________________ +Bool_t +AliFMDSimulator::VMC2FMD(TLorentzVector& v, UShort_t& detector, + Char_t& ring, UShort_t& sector, UShort_t& strip) +{ + TVirtualMC* mc = TVirtualMC::GetMC(); + + // Get track position + mc->TrackPosition(v); + Int_t moduleno; mc->CurrentVolOffID(fModuleOff, moduleno); + Int_t iring; mc->CurrentVolOffID(fRingOff, iring); ring = Char_t(iring); + Int_t det; mc->CurrentVolOffID(fDetectorOff, det); detector = det; + + + // Get the ring geometry + AliFMDGeometry* fmd = AliFMDGeometry::Instance(); + //Int_t nsec = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors(); + Int_t nstr = fmd->GetDetector(detector)->GetRing(ring)->GetNStrips(); + Double_t lowr = fmd->GetDetector(detector)->GetRing(ring)->GetLowR(); + Double_t highr= fmd->GetDetector(detector)->GetRing(ring)->GetHighR(); + Double_t theta= fmd->GetDetector(detector)->GetRing(ring)->GetTheta(); + + // Figure out the strip number + Double_t r = TMath::Sqrt(v.X() * v.X() + v.Y() * v.Y()); + Double_t pitch = (highr - lowr) / nstr; + Int_t str = Int_t((r - lowr) / pitch); + if (str < 0 || str >= nstr) return kFALSE; + strip = str; + + // Figure out the sector number + Double_t phi = TMath::ATan2(v.Y(), v.X()) * 180. / TMath::Pi(); + if (phi < 0) phi = 360. + phi; + Double_t t = phi - 2 * moduleno * theta; + sector = 2 * moduleno; + if (t < 0 || t > 2 * theta) return kFALSE; + else if (t > theta) sector += 1; + + AliDebug(40, Form("<1> Inside an active FMD volume FMD%d%c[%2d,%3d] %s", + detector, ring, sector, strip, mc->CurrentVolPath())); + return kTRUE; +} + +//____________________________________________________________________ +Bool_t +AliFMDSimulator::VMC2FMD(Int_t copy, TLorentzVector& v, + UShort_t& detector, Char_t& ring, + UShort_t& sector, UShort_t& strip) +{ + TVirtualMC* mc = TVirtualMC::GetMC(); + + strip = copy - 1; + Int_t sectordiv; mc->CurrentVolOffID(fSectorOff, sectordiv); + if (fModuleOff >= 0) { + Int_t module; mc->CurrentVolOffID(fModuleOff, module); + sector = 2 * module + sectordiv; + } + else + sector = sectordiv; + Int_t iring; mc->CurrentVolOffID(fRingOff, iring); ring = Char_t(iring); + Int_t det; mc->CurrentVolOffID(fDetectorOff, det); detector = det; + + AliFMDGeometry* fmd = AliFMDGeometry::Instance(); + //Double_t rz = fmd->GetDetector(detector)->GetRingZ(ring); + Int_t n = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors(); +#if 0 + if (rz < 0) { + Int_t s = ((n - sector + n / 2) % n) + 1; + AliDebug(1, Form("Recalculating sector to %d (=%d-%d+%d/2%%%d+1 z=%f)", + s, n, sector, n, n, rz)); + sector = s; + } +#endif + if (sector < 1 || sector > n) { + Warning("Step", "sector # %d out of range (0-%d)", sector-1, n-1); + return kFALSE; + } + sector--; + // Get track position + mc->TrackPosition(v); + AliDebug(15, Form("<2> Inside an active FMD volume FMD%d%c[%2d,%3d] %s", + detector, ring, sector, strip, mc->CurrentVolPath())); + + return kTRUE; +} + //____________________________________________________________________ void AliFMDSimulator::Exec(Option_t* /* option */) @@ -367,89 +491,136 @@ AliFMDSimulator::Exec(Option_t* /* option */) // - ENDIF // TVirtualMC* mc = TVirtualMC::GetMC(); - if (!mc->IsTrackAlive()) return; - if (TMath::Abs(mc->TrackCharge()) <= 0) return; + Double_t absQ = TMath::Abs(mc->TrackCharge()); + if (absQ <= 0) return; Int_t copy; Int_t vol = mc->CurrentVolID(copy); - if (vol != fInnerId && vol != fOuterId) { - AliDebug(15, Form("Not an FMD volume %d '%s' (%d or %d)", - vol, mc->CurrentVolName(), fInnerId, fOuterId)); + if (!IsActive(vol)) { + AliDebug(50, Form("Not an FMD volume %d '%s'",vol,mc->CurrentVolName())); return; } - + TLorentzVector v; + UShort_t detector; + Char_t ring; + UShort_t sector; + UShort_t strip; + + if (fUseDivided) { + if (!VMC2FMD(copy, v, detector, ring, sector, strip)) return; + } else { + if (!VMC2FMD(v, detector, ring, sector, strip)) return; + } + TLorentzVector p; + mc->TrackMomentum(p); + Int_t trackno = gAlice->GetMCApp()->GetCurrentTrackNumber(); + Int_t pdg = mc->TrackPid(); + Double_t mass = mc->TrackMass(); + Double_t edep = mc->Edep() * 1000; // keV + Double_t poverm = (mass == 0 ? 0 : p.P() / mass); + Bool_t isBad = kFALSE; + + // This `if' is to debug abnormal energy depositions. We trigger on + // p/m approx larger than or equal to a MIP, and a large edep - more + // than 1 keV - a MIP is 100 eV. + if (edep > absQ * absQ && poverm > 1) { + isBad = kTRUE; + TArrayI procs; + mc->StepProcesses(procs); + TString processes; + for (Int_t ip = 0; ip < procs.fN; ip++) { + if (ip != 0) processes.Append(","); + processes.Append(TMCProcessName[procs.fArray[ip]]); + } + TDatabasePDG* pdgDB = TDatabasePDG::Instance(); + TParticlePDG* particleType = pdgDB->GetParticle(pdg); + TString pname(particleType ? particleType->GetName() : "???"); + TString what; + if (mc->IsTrackEntering()) what.Append("entering "); + if (mc->IsTrackExiting()) what.Append("exiting "); + if (mc->IsTrackInside()) what.Append("inside "); + if (mc->IsTrackDisappeared()) what.Append("disappeared "); + if (mc->IsTrackStop()) what.Append("stopped "); + if (mc->IsNewTrack()) what.Append("new "); + if (mc->IsTrackAlive()) what.Append("alive "); + if (mc->IsTrackOut()) what.Append("out "); + + Int_t mother = gAlice->GetMCApp()->GetPrimary(trackno); + Warning("Step", "Track # %5d deposits a lot of energy\n" + " Volume: %s\n" + " Momentum: (%7.4f,%7.4f,%7.4f)\n" + " PDG: %d (%s)\n" + " Edep: %-14.7f keV (mother %d)\n" + " p/m: %-7.4f/%-7.4f = %-14.7f\n" + " Processes: %s\n" + " What: %s\n", + trackno, mc->CurrentVolPath(), p.X(), p.Y(), p.Z(), + pdg, pname.Data(), edep, mother, p.P(), mass, + poverm, processes.Data(), what.Data()); + } + // Check that the track is actually within the active area Bool_t entering = mc->IsTrackEntering(); Bool_t inside = mc->IsTrackInside(); Bool_t out = (mc->IsTrackExiting()|| mc->IsTrackDisappeared()|| mc->IsTrackStop()); - // Reset the energy deposition for this track, and update some of // our parameters. if (entering) { - AliDebug(15, "Entering active FMD volume"); - fCurrentDeltaE = 0; - - // Get production vertex and momentum of the track - mc->TrackMomentum(fCurrentP); - mc->TrackPosition(fCurrentV); - fCurrentPdg = mc->IdFromPDG(mc->TrackPid()); + AliDebug(15, Form("Track # %8d entering active FMD volume %s: " + "Edep=%f (%f,%f,%f)", trackno, mc->CurrentVolPath(), + edep, v.X(), v.Y(), v.Z())); + fCurrentP = p; + fCurrentV = v; + fCurrentDeltaE = edep; + fCurrentPdg = pdg; // mc->IdFromPDG(pdg); } - // If the track is inside, then update the energy deposition - if (inside && fCurrentDeltaE >= 0) - AliDebug(15, "Inside active FMD volume"); - fCurrentDeltaE += 1000 * mc->Edep(); - + if (inside && fCurrentDeltaE >= 0) { + fCurrentDeltaE += edep; + AliDebug(15, Form("Track # %8d inside active FMD volume %s: Edep=%f, " + "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) { - AliDebug(15, Form("Leaving active FMD volume %s", mc->CurrentVolPath())); - - Int_t strip = copy - 1; - Int_t sectordiv; - mc->CurrentVolOffID(fSectorOff, sectordiv); - Int_t module; - mc->CurrentVolOffID(fModuleOff, module); - Int_t sector = 2 * module + sectordiv; - Int_t iring; - mc->CurrentVolOffID(fRingOff, iring); - Char_t ring = Char_t(iring); - Int_t detector; - mc->CurrentVolOffID(fDetectorOff, detector); - - - AliFMDGeometry* fmd = AliFMDGeometry::Instance(); - Double_t rz = fmd->GetDetector(detector)->GetRingZ(ring); - Int_t n = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors(); - if (rz < 0) { - Int_t s = ((n - sector + n / 2) % n) + 1; - AliDebug(40, Form("Recalculating sector to %d (=%d-%d+%d/2%%%d+1 z=%f)", - s, n, sector, n, n, rz)); - sector = s; - } - if (sector < 1 || sector > n) { - Warning("Step", "sector # %d out of range (0-%d)", sector-1, n-1); - return; + 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); + } } - sector--; - fCurrentDeltaE += 1000 * mc->Edep(); - - AliDebug(20, Form("Processing hit in FMD%d%c[%2d,%3d]: %f", - detector, ring, sector, strip, fCurrentDeltaE)); - - fFMD->AddHitByFields(gAlice->GetMCApp()->GetCurrentTrackNumber(), - UShort_t(detector), ring, UShort_t(sector), - UShort_t(strip), - fCurrentV.X(), fCurrentV.Y(), fCurrentV.Z(), - fCurrentP.X(), fCurrentP.Y(), fCurrentP.Z(), - fCurrentDeltaE, fCurrentPdg, fCurrentV.T()); 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(next()))) + hit->Print("D"); + fBad->Clear(); + } +} //____________________________________________________________________