Removing the tasks from the digitization (Ruben)
[u/mrichter/AliRoot.git] / FMD / AliFMDHitDigitizer.cxx
index 659e6b2..ad7d7f3 100644 (file)
 #include <AliLoader.h>         // ALILOADER_H
 #include <AliRunLoader.h>      // ALIRUNLOADER_H
 #include <AliFMDHit.h>
+#include <AliStack.h>
 #include <TFile.h>
+#include <TParticle.h>
 
 //====================================================================
 ClassImp(AliFMDHitDigitizer)    
@@ -209,17 +211,36 @@ ClassImp(AliFMDHitDigitizer)
 
 //____________________________________________________________________
 AliFMDHitDigitizer::AliFMDHitDigitizer(AliFMD* fmd, Output_t  output)
-  : AliFMDBaseDigitizer("FMD", (fOutput == kDigits ? 
+  : AliFMDBaseDigitizer("FMD", (output == kDigits ? 
                                "FMD Hit->Digit digitizer" :
                                "FMD Hit->SDigit digitizer")),
-    fOutput(output)
+    fOutput(output), 
+    fHoldTime(2e-6),
+    fStack(0)
 {
   fFMD = fmd;
 }
 
 //____________________________________________________________________
+AliFMDHitDigitizer& 
+AliFMDHitDigitizer::operator=(const AliFMDHitDigitizer& o) 
+{
+  /** 
+   * Assignment operator
+   *
+   * @param o Object to assign from 
+   * @return Reference to this 
+   */
+  AliFMDBaseDigitizer::operator=(o);
+  fHoldTime    = o.fHoldTime;
+  fOutput      = o.fOutput;
+  fStack       = o.fStack;
+  return *this;
+}
+
+//____________________________________________________________________
 void
-AliFMDHitDigitizer::Exec(Option_t* /*option*/)
+AliFMDHitDigitizer::Digitize(Option_t* /*option*/)
 {
   // Run this digitizer 
   // Get an inititialize parameter manager
@@ -262,10 +283,21 @@ AliFMDHitDigitizer::Exec(Option_t* /*option*/)
     }
     
     // Read in the event
-    AliFMDDebug(5, ("Now digitizing (Hits->%s) event # %d", 
+    AliFMDDebug(1, ("Now digitizing (Hits->%s) event # %d", 
                    (fOutput == kDigits ? "digits" : "sdigits"), event));
     thisLoader->GetEvent(event);
     
+    // Load kinematics to get primary information for SDigits
+    fStack = 0;
+    if (fOutput == kSDigits) {
+      if (thisLoader->LoadKinematics("READ")) {
+       AliError("Failed to get kinematics from event loader");
+       return;
+      }
+      AliFMDDebug(5, ("Loading stack of kinematics"));
+      fStack = thisLoader->Stack();
+    }
+
     // Check that we have the hits 
     if (!loader->TreeH() && loader->LoadHits()) {
       AliError("Failed to load hits");
@@ -310,7 +342,7 @@ AliFMDHitDigitizer::Exec(Option_t* /*option*/)
     
     // Write digits to tree
     Int_t write = outTree->Fill();
-    AliFMDDebug(1, ("Wrote %d bytes to digit tree", write));
+    AliFMDDebug(5, ("Wrote %d bytes to digit tree", write));
 
     // Store the digits
     StoreDigits(loader);
@@ -322,6 +354,12 @@ AliFMDHitDigitizer::Exec(Option_t* /*option*/)
 TTree*
 AliFMDHitDigitizer::MakeOutputTree(AliLoader* loader)
 {
+  /** 
+   * Make the output tree using the passed loader 
+   *
+   * @param loader 
+   * @return The generated tree. 
+   */
   if (fOutput == kDigits) 
     return AliFMDBaseDigitizer::MakeOutputTree(loader);
   
@@ -363,24 +401,73 @@ AliFMDHitDigitizer::SumContributions(TBranch* hitsBranch)
   for (Int_t track = 0; track < ntracks; track++)  {
     // Read in entry number `track' 
     read += hitsBranch->GetEntry(track);
-    
+
     // Get the number of hits 
     Int_t nhits = fmdHits->GetEntries ();
     for (Int_t hit = 0; hit < nhits; hit++) {
       // Get the hit number `hit'
       AliFMDHit* fmdHit = 
        static_cast<AliFMDHit*>(fmdHits->UncheckedAt(hit));
+
+      // Ignore hits that arrive too late
+      if (fmdHit->Time() > fHoldTime) continue;
       
+
+      // Check if this is a primary particle
+      Bool_t isPrimary = kTRUE;
+      Int_t  trackno   = -1;
+      if (fStack) {
+       trackno = fmdHit->Track();
+       AliFMDDebug(10, ("Will get track # %d/%d from entry # %d", 
+                       trackno, fStack->GetNtrack(), track));
+       if (fStack->GetNtrack() < trackno) {
+         AliError(Form("Track number %d/%d out of bounds", 
+                       trackno, fStack->GetNtrack()));
+         continue;
+       }
+#if 1
+       isPrimary = fStack->IsPhysicalPrimary(trackno);
+#else // This is our hand-crafted code.  We use the ALICE definition
+       TParticle* part    = fStack->Particle(trackno);
+       isPrimary          = part->IsPrimary();
+       if (!isPrimary) { 
+         // Extended testing of mother status - this is for Pythia6.
+         Int_t      mother1   = part->GetFirstMother();
+         TParticle* mother    = fStack->Particle(mother1);
+         if (!mother || mother->GetStatusCode() > 1)
+           isPrimary = kTRUE;
+         AliFMDDebug(15,
+                     ("Track %d secondary, mother: %d - %s - status %d: %s", 
+                      trackno, mother1, 
+                      (mother ? "found"                 : "not found"), 
+                      (mother ? mother->GetStatusCode() : -1),
+                      (isPrimary ? "primary" : "secondary")));
+       }
+#endif
+      }
+    
       // Extract parameters 
+      AliFMDDebug(15,("Adding contribution %7.5f for FMD%d%c[%2d,%3d] "
+                     " for trackno %6d (%s)", 
+                     fmdHit->Edep(),
+                     fmdHit->Detector(), 
+                     fmdHit->Ring(),
+                     fmdHit->Sector(), 
+                     fmdHit->Strip(), 
+                     trackno, 
+                     (isPrimary ? "primary" : "secondary")));
       AddContribution(fmdHit->Detector(),
                      fmdHit->Ring(),
                      fmdHit->Sector(),
                      fmdHit->Strip(),
-                     fmdHit->Edep());
+                     fmdHit->Edep(), 
+                     isPrimary, 
+                     1, 
+                     &trackno);
     }  // hit loop
   } // track loop
   AliFMDDebug(5, ("Size of cache: %d bytes, read %d bytes", 
-                  sizeof(fEdep), read));
+                 int(sizeof(fEdep)), read));
 }
 
 
@@ -400,26 +487,48 @@ AliFMDHitDigitizer::MakePedestal(UShort_t  detector,
 
 //____________________________________________________________________
 void
-AliFMDHitDigitizer::AddDigit(UShort_t  detector, 
-                            Char_t    ring,
-                            UShort_t  sector, 
-                            UShort_t  strip, 
-                            Float_t   edep, 
-                            UShort_t  count1, 
-                            Short_t   count2, 
-                            Short_t   count3,
-                            Short_t   count4) const
+AliFMDHitDigitizer::AddDigit(UShort_t        detector, 
+                            Char_t          ring,
+                            UShort_t        sector, 
+                            UShort_t        strip, 
+                            Float_t         edep, 
+                            UShort_t        count1, 
+                            Short_t         count2, 
+                            Short_t         count3,
+                            Short_t         count4, 
+                            UShort_t        ntotal,
+                            UShort_t        nprim, 
+                            const TArrayI&  refs) const
 {
   // Add a digit or summable digit
   if (fOutput == kDigits) { 
+    AliFMDDebug(15,("Adding digit for FMD%d%c[%2d,%3d] = (%x,%x,%x,%x)",
+                   detector, ring, sector, strip, 
+                   count1, count2, count3, count4));
     AliFMDBaseDigitizer::AddDigit(detector, ring, sector, strip, 0,
-                                 count1, count2, count3, count4);
+                                 count1, count2, count3, count4, 
+                                 ntotal, nprim, refs);
+    return;
+  }
+  if (edep <= 0) { 
+    AliFMDDebug(15, ("Digit edep = %f <= 0 for FMD%d%c[%2d,%3d]", 
+                   edep, detector, ring, sector, strip));
+    return;
+  }
+  if (count1 == 0 && count2 <= 0 && count3 <= 0 && count4 <= 0) {
+    AliFMDDebug(15, ("Digit counts = (%x,%x,%x,%x) <= 0 for FMD%d%c[%2d,%3d]", 
+                   count1, count2, count3, count4, 
+                   detector, ring, sector, strip));
     return;
   }
-  if (edep <= 0) return;
-  if (count1 == 0 && count2 <= 0 && count3 <= 0 && count4 <= 0) return;
+  AliFMDDebug(15, ("Adding sdigit for FMD%d%c[%2d,%3d] = "
+                  "(%x,%x,%x,%x) [%d/%d] %d",
+                  detector, ring, sector, strip, 
+                  count1, count2, count3, count4, nprim, ntotal, refs.fN));
   fFMD->AddSDigitByFields(detector, ring, sector, strip, edep,
-                         count1, count2, count3, count4);
+                         count1, count2, count3, count4, 
+                         ntotal, nprim, fStoreTrackRefs ? refs.fArray : 0);
+  if (fStoreTrackRefs && nprim > 3) fIgnoredLabels += nprim - 3;
 }
 
 //____________________________________________________________________
@@ -454,8 +563,13 @@ AliFMDHitDigitizer::CheckDigit(AliFMDDigit*    digit,
 
 //____________________________________________________________________
 void
-AliFMDHitDigitizer::StoreDigits(AliLoader* loader)
+AliFMDHitDigitizer::StoreDigits(const AliLoader* loader)
 {
+  /** 
+   * Store the data using the loader 
+   *
+   * @param loader The loader 
+   */
   if (fOutput == kDigits) { 
     AliFMDBaseDigitizer::StoreDigits(loader);
     return;