Commiting before moving code from `AliFMDSimulator' to `AliFMD' (and
[u/mrichter/AliRoot.git] / FMD / AliFMDSimulator.cxx
index 19378586c1b2834f31b1167d178424f2b3243ca0..2cacd68d07040303fb1c7d50c5c866da676f1ecc 100644 (file)
 //                 |                           |             
 //        +--------------------+   +-------------------+
 //        | AliFMDGeoSimulator |   | AliFMDG3Simulator | 
-//        +--------------------+   +---------+---------+
+//        +--------------------+   +-------------------+
+//                                           ^
+//                                           |
+//                                 +----------------------+
+//                                | AliFMDG3OldSimulator |
+//                                +----------------------+
 //      
-//
 // *  AliFMD 
 //    This defines the interface for the various parts of AliROOT that
 //    uses the FMD, like AliFMDSimulator, AliFMDDigitizer, 
 //    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>            // 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
@@ -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,88 +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->AddHit(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<AliFMDHit*>(next()))) 
+      hit->Print("D");
+    fBad->Clear();
+  }
+}
 
 
 //____________________________________________________________________