Minor changes, mostly for debugging
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 2 Nov 2005 15:42:46 +0000 (15:42 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 2 Nov 2005 15:42:46 +0000 (15:42 +0000)
FMD/AliFMD.cxx
FMD/AliFMD.h
FMD/AliFMDG3OldSimulator.cxx
FMD/AliFMDG3Simulator.cxx
FMD/AliFMDGeoOldSimulator.cxx
FMD/AliFMDGeoSimulator.cxx
FMD/AliFMDHit.cxx
FMD/AliFMDHit.h
FMD/AliFMDSimulator.cxx
FMD/AliFMDSimulator.h

index b80c7ee..f07b1f9 100644 (file)
@@ -180,7 +180,7 @@ AliFMD::AliFMD(const char *name, const char *title)
   fUseDivided  = kFALSE;
   fUseAssembly = kFALSE;
   fUseGeo      = kTRUE;
-
+  
   // Initialise Hit array
   HitsArray();
   gAlice->GetMCApp()->AddHitList(fHits);
@@ -327,6 +327,14 @@ AliFMD::Init()
   //
 }
 
+//____________________________________________________________________
+void
+AliFMD::FinishEvent()
+{
+  if (fSimulator) fSimulator->EndEvent();
+}
+
+
 //====================================================================
 //
 // Graphics and event display
@@ -606,7 +614,7 @@ AliFMD::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 }
 
 //____________________________________________________________________
-void 
+AliFMDHit*
 AliFMD::AddHitByFields(Int_t    track, 
                       UShort_t detector, 
                       Char_t   ring, 
@@ -662,13 +670,14 @@ AliFMD::AddHitByFields(Int_t    track,
              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;
 }
 
 //____________________________________________________________________
index 12df636..99786a7 100644 (file)
@@ -23,6 +23,7 @@ class TClonesArray;
 class TBrowser;
 class AliDigitizer;
 class AliFMDSimulator;
+class AliFMDHit;
 
 //____________________________________________________________________
 class AliFMD : public AliDetector 
@@ -38,14 +39,15 @@ public:
   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();
@@ -58,7 +60,7 @@ public:
   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, 
index 57b2c6a..f9d1196 100644 (file)
@@ -140,25 +140,26 @@ AliFMDG3OldSimulator::RingGeometry(AliFMDRing* r)
   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;
@@ -167,7 +168,7 @@ AliFMDG3OldSimulator::RingGeometry(AliFMDRing* r)
   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);
 
@@ -185,7 +186,7 @@ AliFMDG3OldSimulator::RingGeometry(AliFMDRing* r)
   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());
@@ -210,12 +211,12 @@ AliFMDG3OldSimulator::RingGeometry(AliFMDRing* r)
 
   // 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);
@@ -229,16 +230,14 @@ AliFMDG3OldSimulator::RingGeometry(AliFMDRing* r)
   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);
 
index aa0c43f..2778cdd 100644 (file)
@@ -167,8 +167,8 @@ AliFMDG3Simulator::RingGeometry(AliFMDRing* r)
   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;
index 02820fa..398bea0 100644 (file)
@@ -143,37 +143,29 @@ AliFMDGeoOldSimulator::RingGeometry(AliFMDRing* r)
   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. 
@@ -209,15 +201,15 @@ AliFMDGeoOldSimulator::RingGeometry(AliFMDRing* r)
 
   // 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);
@@ -226,14 +218,14 @@ AliFMDGeoOldSimulator::RingGeometry(AliFMDRing* r)
   // 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));
 
index c64d6ff..f31fcd9 100644 (file)
@@ -251,9 +251,9 @@ AliFMDGeoSimulator::RingGeometry(AliFMDRing* r)
     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");
@@ -293,7 +293,7 @@ AliFMDGeoSimulator::RingGeometry(AliFMDRing* r)
   
   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), 
@@ -361,16 +361,16 @@ AliFMDGeoSimulator::RingGeometry(AliFMDRing* r)
   // 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), 
index 09ab2b1..41bfc7b 100644 (file)
@@ -24,6 +24,9 @@
 #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)
@@ -103,8 +106,38 @@ AliFMDHit::AliFMDHit(Int_t    shunt,
 }
 
 //____________________________________________________________________
+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" 
@@ -112,6 +145,15 @@ AliFMDHit::Print(Option_t* /* option */) const
        << 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;
+  }
 }
 
 //____________________________________________________________________
index 9050690..8cddebe 100644 (file)
@@ -46,6 +46,9 @@ public:
   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;
index b441c2c..4ef549c 100644 (file)
@@ -80,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
@@ -135,7 +137,8 @@ AliFMDSimulator::AliFMDSimulator()
     fOuterId(-1), 
     fActiveId(4), 
     fUseDivided(kFALSE),
-    fUseAssembly(kTRUE)
+    fUseAssembly(kTRUE), 
+    fBad(0)
 {
   // Default constructor
 }
@@ -149,7 +152,8 @@ AliFMDSimulator::AliFMDSimulator(AliFMD* fmd, Bool_t detailed)
     fOuterId(-1),
     fActiveId(4),
     fUseDivided(kFALSE),
-    fUseAssembly(kTRUE)
+    fUseAssembly(kTRUE),
+    fBad(0)
 {
   // Normal constructor
   // 
@@ -158,6 +162,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");
 }
 
 
@@ -428,7 +433,7 @@ AliFMDSimulator::VMC2FMD(Int_t copy, TLorentzVector& v,
   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) {
@@ -445,7 +450,7 @@ AliFMDSimulator::VMC2FMD(Int_t copy, TLorentzVector& v,
   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;
@@ -493,8 +498,7 @@ AliFMDSimulator::Exec(Option_t* /* option */)
   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;
@@ -565,35 +569,58 @@ AliFMDSimulator::Exec(Option_t* /* option */)
   // 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();
+  }
+}
 
 
 //____________________________________________________________________
index a2b4e73..9e31d94 100644 (file)
@@ -23,6 +23,7 @@ class AliFMDDetector;
 class AliFMD1;
 class AliFMD2;
 class AliFMD3;
+class TObjArray;
 
 /** Simulation of the FMD. 
     This class builds the geometry, and processes hits in the FMD */ 
@@ -39,6 +40,7 @@ public:
   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:  
@@ -96,6 +98,7 @@ 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
 };