Added code to put N_primary and N_total into all SDigits.
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 9 Dec 2008 05:45:53 +0000 (05:45 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 9 Dec 2008 05:45:53 +0000 (05:45 +0000)
This is for background studies, so that we don't need to store hits.
In the not too distant future, we'll also add track references.  See also
FMD/scripts/DrawSDigits.C.

Hans, perhaps you could check this code.

25 files changed:
FMD/AliFMD.cxx
FMD/AliFMD.h
FMD/AliFMDBaseDigitizer.cxx
FMD/AliFMDBaseDigitizer.h
FMD/AliFMDDigitizer.cxx
FMD/AliFMDEdepHitPair.h
FMD/AliFMDEdepMap.cxx
FMD/AliFMDHitDigitizer.cxx
FMD/AliFMDHitDigitizer.h
FMD/AliFMDSDigit.cxx
FMD/AliFMDSDigit.h
FMD/Calib/AltroMap/Run0_999999_v0_s0.root
FMD/Calib/Dead/Run0_999999_v0_s0.root
FMD/Calib/Pedestal/Run0_999999_v0_s0.root
FMD/Calib/PulseGain/Run0_999999_v0_s0.root
FMD/Calib/SampleRate/Run0_999999_v0_s0.root
FMD/Calib/StripRange/Run0_999999_v0_s0.root
FMD/Calib/ZeroSuppression/Run0_999999_v0_s0.root
FMD/Config.C
FMD/Simulate.C
FMD/scripts/Compile.C
FMD/scripts/DrawSDigits.C [new file with mode: 0644]
FMD/scripts/Hits2SDigits.C
FMD/scripts/MakeCalibration.C
FMD/scripts/PatternSDigits.C

index 7ec9355de5a40f1d86fe1b931a0a5fdcba9ed3a5..3c493b6972055d1235261e38785e51fd214025a0 100644 (file)
@@ -956,15 +956,18 @@ AliFMD::AddSDigit(Int_t* digits)
   //    digits[6]  [Short_t]  ADC Count, -1 if not used
   //    digits[7]  [Short_t]  ADC Count, -1 if not used 
   // 
-  AddSDigitByFields(UShort_t(digits[0]),  // Detector #
-                   Char_t(digits[1]),    // Ring ID
-                   UShort_t(digits[2]),  // Sector #
-                   UShort_t(digits[3]),  // Strip #
-                   Float_t(digits[4]),   // Edep
-                   UShort_t(digits[5]),  // ADC Count1 
-                   Short_t(digits[6]),   // ADC Count2 
-                   Short_t(digits[7]),   // ADC Count3 
-                   Short_t(digits[8]));
+  AddSDigitByFields(UShort_t(digits[0]),   // Detector #
+                   Char_t(digits[1]),     // Ring ID
+                   UShort_t(digits[2]),   // Sector #
+                   UShort_t(digits[3]),   // Strip #
+                   Float_t(digits[4]),    // Edep
+                   UShort_t(digits[5]),   // ADC Count1 
+                   Short_t(digits[6]),    // ADC Count2 
+                   Short_t(digits[7]),    // ADC Count3 
+                   Short_t(digits[8]),    // ADC Count4
+                   UShort_t(digits[9]),   // N particles
+                   UShort_t(digits[10])); // N primaries
+   
 }
 
 //____________________________________________________________________
@@ -977,7 +980,9 @@ AliFMD::AddSDigitByFields(UShort_t detector,
                          UShort_t count1, 
                          Short_t  count2,
                          Short_t  count3, 
-                         Short_t  count4)
+                         Short_t  count4, 
+                         UShort_t ntot, 
+                         UShort_t nprim)
 {
   // add a summable digit
   // 
@@ -997,7 +1002,7 @@ AliFMD::AddSDigitByFields(UShort_t detector,
   
   new (a[fNsdigits++]) 
     AliFMDSDigit(detector, ring, sector, strip, edep, 
-                count1, count2, count3, count4);
+                count1, count2, count3, count4, ntot, nprim);
 }
 
 //____________________________________________________________________
index d127efa580fb1b08ee6a0b19c341d3f3daf1fa5d..449761f2ab04756772896f3b50ba4c2f6c625f42 100644 (file)
@@ -456,9 +456,14 @@ public:
       - digits[1]  [Char_t]   Ring ID
       - digits[2]  [UShort_t] Sector #
       - digits[3]  [UShort_t] Strip #
-      - digits[4]  [UShort_t] ADC Count 
-      - digits[5]  [Short_t]  ADC Count, -1 if not used
-      - digits[6]  [Short_t]  ADC Count, -1 if not used  */
+      - digits[4]  [Float_t]  Edep
+      - digits[5]  [UShort_t] ADC Count 
+      - digits[6]  [Short_t]  ADC Count, -1 if not used
+      - digits[7]  [Short_t]  ADC Count, -1 if not used  
+      - digits[8]  [Short_t]  ADC Count, -1 if not used
+      - digits[9]  [UShort_t] N total particles
+      - digits[10] [UShort_t] N total primary particles
+  */
   virtual        void   AddSDigit(Int_t *digits);
   /** add a summable digit - as coming from data
       @param detector  Detector # (1, 2, or 3)                      
@@ -477,7 +482,9 @@ public:
                                          UShort_t count1=0, 
                                          Short_t  count2=-1, 
                                          Short_t  count3=-1,
-                                         Short_t  count4=-1);
+                                         Short_t  count4=-1, 
+                                         UShort_t ntot=0, 
+                                         UShort_t nprim=0);
   /** @}*/
 
   /** @{ */
index 8134f8bf5c90609bce41371ac5eb6a52071eb9ae..c063bb378a877242017ce700b35bc0861659ab95 100644 (file)
@@ -309,7 +309,8 @@ AliFMDBaseDigitizer::AddContribution(UShort_t detector,
                                     Char_t   ring, 
                                     UShort_t sector, 
                                     UShort_t strip, 
-                                    Float_t  edep)
+                                    Float_t  edep, 
+                                    Bool_t   isPrimary)
 {
   // Add edep contribution from (detector,ring,sector,strip) to cache
   AliFMDParameters* param = AliFMDParameters::Instance();
@@ -333,8 +334,14 @@ AliFMDBaseDigitizer::AddContribution(UShort_t detector,
                    detector, ring, sector, strip));
       
   // Sum energy deposition
-  fEdep(detector, ring, sector, strip).fEdep  += edep;
-  fEdep(detector, ring, sector, strip).fN     += 1;
+  fEdep(detector, ring, sector, strip).fEdep    += edep;
+  fEdep(detector, ring, sector, strip).fN       += 1;
+  if (isPrimary)
+    fEdep(detector, ring, sector, strip).fNPrim += 1;
+  AliFMDDebug(15, ("Adding contribution %f to FMD%d%c[%2d,%3d] (%f)", 
+                 edep, detector, ring, sector, strip,
+                 fEdep(detector, ring, sector, strip).fEdep));
+  
 }
 
 //____________________________________________________________________
@@ -378,10 +385,16 @@ AliFMDBaseDigitizer::DigitizeHits() const
          // VA1_ALICE channel. 
          if (strip % 128 == 0) last = 0;
          
-         Float_t edep = fEdep(detector, ring, sector, strip).fEdep;
+         Float_t  edep  = fEdep(detector, ring, sector, strip).fEdep;
+         UShort_t ntot  = fEdep(detector, ring, sector, strip).fN;
+         UShort_t nprim = fEdep(detector, ring, sector, strip).fNPrim;
+         if (edep > 0)
+           AliFMDDebug(15, ("Edep = %f for FMD%d%c[%2d,%3d]", 
+                            edep, detector, ring, sector, strip));
          ConvertToCount(edep, last, detector, ring, sector, strip, counts);
          last = edep;
          
+
          // The following line was introduced - wrongly - by Peter
          // Hristov.  It _will_ break the digitisation and the
          // following reconstruction.  The behviour of the
@@ -401,7 +414,8 @@ AliFMDBaseDigitizer::DigitizeHits() const
          //   if (edep<=0) continue;
          AddDigit(detector, ring, sector, strip, edep, 
                   UShort_t(counts[0]), Short_t(counts[1]), 
-                  Short_t(counts[2]), Short_t(counts[3]));
+                  Short_t(counts[2]), Short_t(counts[3]), 
+                  ntot, nprim);
          AliFMDDebug(15, ("   Adding digit in FMD%d%c[%2d,%3d]=%d", 
                          detector,ring,sector,strip,counts[0]));
 #if 0
@@ -466,7 +480,8 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
   //                  = E + (l - E) * ext(-B * t)
   // 
   AliFMDParameters* param = AliFMDParameters::Instance();
-  Float_t  convF          = (param->GetDACPerMIP()*param->GetPulseGain(detector,ring,sector,strip)) / param->GetEdepMip();
+  Float_t  convF          = (param->GetDACPerMIP() / param->GetEdepMip() *
+                            param->GetPulseGain(detector,ring,sector,strip));
   Int_t    ped            = MakePedestal(detector,ring,sector,strip);
   Int_t    maxAdc         = param->GetAltroChannelSize()-1;
   if (maxAdc < 0) {
@@ -475,7 +490,7 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
   }
   UShort_t rate           = param->GetSampleRate(detector,ring,sector,strip);
   AliFMDDebug(15, ("Sample rate for FMD%d%c[%2d,%3d] = %d", 
-                 detector, ring, sector, strip, rate));
+                  detector, ring, sector, strip, rate));
   if (rate < 1 || rate > 4) {
     AliWarning(Form("Invalid sample rate for for FMD%d%c[%2d,%3d] = %d", 
                    detector, ring, sector, strip, rate));
@@ -492,6 +507,7 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
                     detector,ring,sector,strip,edep,counts[0],convF,ped));
     return;
   }
+
   
   // Create a pedestal 
   Float_t b = fShapingTime;
@@ -502,6 +518,12 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
     if (a < 0) a = 0;
     counts[i]  = UShort_t(TMath::Min(a, Float_t(maxAdc)));
   }
+  AliFMDDebug(15, ("Converted edep = %f to ADC (%x,%x,%x,%x) "
+                  "[gain: %f=(%f/%f*%f), pedestal: %d, rate: %d]", 
+                  edep, counts[0], counts[1], counts[2], counts[3], 
+                  convF, param->GetDACPerMIP(),param->GetEdepMip(),
+                  param->GetPulseGain(detector,ring,sector,strip), 
+                  ped, rate));
 }
 
 //____________________________________________________________________
@@ -514,7 +536,9 @@ AliFMDBaseDigitizer::AddDigit(UShort_t  detector,
                              UShort_t  count1, 
                              Short_t   count2, 
                              Short_t   count3,
-                             Short_t   count4) const
+                             Short_t   count4,
+                             UShort_t  /* ntot */, 
+                             UShort_t  /* nprim */) const
 {
   // Add a digit or summable digit
   
index 473fbe769852e814e9f76af967653a39ace30e69..636038f4effe9e7e031da2f583eaf4f12e9fc52f 100644 (file)
@@ -224,7 +224,8 @@ protected:
                               Char_t   ring, 
                               UShort_t sector, 
                               UShort_t strip, 
-                              Float_t  edep);
+                              Float_t  edep, 
+                              Bool_t   isPrimary);
   /** Add a digit to output */
   virtual void     AddDigit(UShort_t detector, 
                            Char_t   ring,
@@ -234,7 +235,9 @@ protected:
                            UShort_t count1, 
                            Short_t  count2, 
                            Short_t  count3,
-                           Short_t  count4) const;
+                           Short_t  count4, 
+                           UShort_t ntot, 
+                           UShort_t nprim) const;
   /** Make the output tree using the passed loader 
       @param loader 
       @return The generated tree. */
index 2be36313fa68519e05234dda44ea02e8970383ed..aa68a17a6aade953fa04828f64917fa68a5d322e 100644 (file)
@@ -365,7 +365,7 @@ AliFMDDigitizer::Exec(Option_t*)
   
   // Fill the 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(outFMD);
@@ -377,7 +377,7 @@ AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch)
 {
   // Sum energy deposited contributions from each sdigits in a cache
   // (fEdep).  
-  AliFMDDebug(1, ("Runnin our version of SumContributions"));
+  AliFMDDebug(3, ("Runnin our version of SumContributions"));
 
   // Get a list of hits from the FMD manager 
   TClonesArray *fmdSDigits = fFMD->SDigits();
@@ -393,7 +393,7 @@ AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch)
     
     // Get the number of sdigits 
     Int_t nsdigits = fmdSDigits->GetEntries ();
-    AliFMDDebug(1, ("Got %5d SDigits", nsdigits));
+    AliFMDDebug(3, ("Got %5d SDigits", nsdigits));
     for (Int_t sdigit = 0; sdigit < nsdigits; sdigit++) {
       // Get the sdigit number `sdigit'
       AliFMDSDigit* fmdSDigit = 
@@ -404,7 +404,8 @@ AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch)
                      fmdSDigit->Ring(),
                      fmdSDigit->Sector(),
                      fmdSDigit->Strip(),
-                     fmdSDigit->Edep());
+                     fmdSDigit->Edep(), 
+                     kTRUE);
     }  // sdigit loop
   } // event loop
 
index 33bf40690fe9e3bf1f225c06d612f67a24943273..56e889e87ac60f54311999eb4b5f23e2061feaaf 100644 (file)
 class AliFMDEdepHitPair 
 {
 public:
-  Float_t  fEdep; // summed energy deposition
-  UShort_t fN;    // Number of hits
+  Float_t  fEdep;  // summed energy deposition
+  UShort_t fN;     // Number of hits
+  UShort_t fNPrim; // Number of primaries;
+  
   /** CTOR  */
-  AliFMDEdepHitPair() : fEdep(0), fN(0) {}
+  AliFMDEdepHitPair() : fEdep(0), fN(0), fNPrim(0) {}
   /** DTOR */
   virtual ~AliFMDEdepHitPair() {}
   /** Assignment operator 
       @param o Object to assign from 
       @return Reference to this object */
   AliFMDEdepHitPair& operator=(const AliFMDEdepHitPair& o) 
-  { fEdep = o.fEdep; fN    = o.fN; return *this; }
+  { 
+    fEdep  = o.fEdep; 
+    fN     = o.fN; 
+    fNPrim = o.fNPrim;
+    return *this; 
+  }
   /** Copy CTOR 
       @param o Object to copy from */
-  AliFMDEdepHitPair(const AliFMDEdepHitPair& o) : fEdep(o.fEdep), fN(o.fN) {}
-  ClassDef(AliFMDEdepHitPair, 1)
+  AliFMDEdepHitPair(const AliFMDEdepHitPair& o) 
+    : fEdep(o.fEdep), fN(o.fN), fNPrim(o.fNPrim)
+  {}
+  ClassDef(AliFMDEdepHitPair, 2)
 };
 
 #endif 
index 843a9ca299bb5ed9faef62fe6fe80c08cf7ced4b..84d869976465d9cf456360555df681b2058935a0 100644 (file)
@@ -93,8 +93,9 @@ AliFMDEdepMap::Reset()
 {
   // Reset to zero
   for (Int_t i = 0; i < fTotal; i++) { 
-    fData[i].fEdep = 0; 
-    fData[i].fN = 0; 
+    fData[i].fEdep  = 0; 
+    fData[i].fN     = 0; 
+    fData[i].fNPrim = 0;
   };
 }
 
@@ -104,14 +105,16 @@ AliFMDEdepMap::Reset(const AliFMDEdepHitPair& val)
 {
   // Reset to val
   for (Int_t i = 0; i < fTotal; i++) { 
-    fData[i].fEdep = val.fEdep; 
-    fData[i].fN = val.fN; 
+    fData[i].fEdep  = val.fEdep; 
+    fData[i].fN     = val.fN; 
+    fData[i].fNPrim = val.fNPrim;
   };
 }
 
 //____________________________________________________________________
 AliFMDEdepHitPair& 
-AliFMDEdepMap::operator()(UShort_t det, Char_t ring, UShort_t sec, UShort_t str) 
+AliFMDEdepMap::operator()(UShort_t det, Char_t ring, 
+                         UShort_t sec, UShort_t str) 
 {
   // Get data 
   // 
@@ -128,7 +131,8 @@ AliFMDEdepMap::operator()(UShort_t det, Char_t ring, UShort_t sec, UShort_t str)
 
 //____________________________________________________________________
 const AliFMDEdepHitPair& 
-AliFMDEdepMap::operator()(UShort_t det, Char_t ring, UShort_t sec, UShort_t str) const
+AliFMDEdepMap::operator()(UShort_t det, Char_t ring, 
+                         UShort_t sec, UShort_t str) const
 {
   // Get data 
   // 
index 659e6b24ffdce350ae450c7691e8b23ff4fc1502..4087677ec19fcc2b703d0cf6fe73104f298bbc2a 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)    
@@ -212,7 +214,8 @@ AliFMDHitDigitizer::AliFMDHitDigitizer(AliFMD* fmd, Output_t  output)
   : AliFMDBaseDigitizer("FMD", (fOutput == kDigits ? 
                                "FMD Hit->Digit digitizer" :
                                "FMD Hit->SDigit digitizer")),
-    fOutput(output)
+    fOutput(output), 
+    fStack(0)
 {
   fFMD = fmd;
 }
@@ -266,6 +269,17 @@ AliFMDHitDigitizer::Exec(Option_t* /*option*/)
                    (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 +324,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);
@@ -363,7 +377,7 @@ 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++) {
@@ -371,12 +385,42 @@ AliFMDHitDigitizer::SumContributions(TBranch* hitsBranch)
       AliFMDHit* fmdHit = 
        static_cast<AliFMDHit*>(fmdHits->UncheckedAt(hit));
       
+      // Check if this is a primary particle
+      Bool_t isPrimary = kTRUE;
+      if (fStack) {
+       Int_t      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;
+       }
+       
+       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")));
+       }
+      }
+    
       // Extract parameters 
       AddContribution(fmdHit->Detector(),
                      fmdHit->Ring(),
                      fmdHit->Sector(),
                      fmdHit->Strip(),
-                     fmdHit->Edep());
+                     fmdHit->Edep(), 
+                     isPrimary);
     }  // hit loop
   } // track loop
   AliFMDDebug(5, ("Size of cache: %d bytes, read %d bytes", 
@@ -408,18 +452,33 @@ AliFMDHitDigitizer::AddDigit(UShort_t  detector,
                             UShort_t  count1, 
                             Short_t   count2, 
                             Short_t   count3,
-                            Short_t   count4) const
+                            Short_t   count4, 
+                            UShort_t  ntotal,
+                            UShort_t  nprim) const
 {
   // Add a digit or summable digit
   if (fOutput == kDigits) { 
     AliFMDBaseDigitizer::AddDigit(detector, ring, sector, strip, 0,
-                                 count1, count2, count3, count4);
+                                 count1, count2, count3, count4, 0, 0);
+    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 digit for FMD%d%c[%2d,%3d] = (%x,%x,%x,%x) [%d/%d]",
+                  detector, ring, sector, strip, 
+                  count1, count2, count3, count4, nprim, ntotal));
   fFMD->AddSDigitByFields(detector, ring, sector, strip, edep,
-                         count1, count2, count3, count4);
+                         count1, count2, count3, count4, 
+                         ntotal, nprim);
 }
 
 //____________________________________________________________________
index 614c36d904048f5aaea315a2640ae6cf8f57fa6d..656f94eba35c8abe69578a8d5bbaa69f198c9c25 100644 (file)
@@ -30,6 +30,7 @@ class AliFMD;
 class AliLoader;
 class AliRunLoader;
 class AliFMDDigit;
+class AliStack;
 
 
 
@@ -50,7 +51,8 @@ public:
   /** CTOR */
   AliFMDHitDigitizer() 
     : AliFMDBaseDigitizer(), 
-      fOutput(kDigits) 
+      fOutput(kDigits), 
+      fStack(0)
   {}
   /** CTOR 
       @param name Name */
@@ -64,7 +66,8 @@ protected:
       @param o Object to copy from */
   AliFMDHitDigitizer(const AliFMDHitDigitizer& o) 
     : AliFMDBaseDigitizer(o),
-      fOutput(o.fOutput)
+      fOutput(o.fOutput), 
+      fStack(o.fStack)
   {}
   /** Assignment operator
       @param o Object to assign from 
@@ -111,7 +114,9 @@ protected:
                UShort_t  count1, 
                Short_t   count2, 
                Short_t   count3,
-               Short_t   count4) const;
+               Short_t   count4, 
+               UShort_t  ntot, 
+               UShort_t  nprim) const;
   /** Check that digit data is consistent
       @param digit   Digit
       @param nhits   Number of hits
@@ -125,6 +130,7 @@ protected:
   
 
   Output_t      fOutput;           // Output mode
+  AliStack*     fStack;            // Kinematics
 
   ClassDef(AliFMDHitDigitizer,1) // Make Digits from Hits
 };
index eb3fa44d3c6251de6e063709c92726002c4a5331..ae98a01adaf19d1e09f7911793ed3fc2024593d1 100644 (file)
@@ -92,13 +92,17 @@ AliFMDSDigit::AliFMDSDigit(UShort_t detector,
                           UShort_t count1,
                           Short_t  count2, 
                           Short_t  count3, 
-                          Short_t  count4)
+                          Short_t  count4,
+                          UShort_t npart,
+                          UShort_t nprim)
   : AliFMDBaseDigit(detector, ring, sector, strip), 
     fEdep(edep),
     fCount1(count1),
     fCount2(count2),
     fCount3(count3),
-    fCount4(count4)
+    fCount4(count4),
+    fNParticles(npart), 
+    fNPrimaries(nprim)
 {
   //
   // Creates a real data digit object
index 87003f429febc111107017839ba1cd2710698bd1..6e6b1387ac34b6a8739b4994d67ec6d44630c882 100644 (file)
@@ -42,7 +42,9 @@ public:
               UShort_t count=0, 
               Short_t  count2=-1, 
               Short_t  count3=-1,
-              Short_t  count4=-1);
+              Short_t  count4=-1,
+              UShort_t npart=0,
+              UShort_t nprim=0);
   /** DTOR */
   virtual ~AliFMDSDigit() {}
   /** @return ADC count (first sample) */
@@ -57,6 +59,11 @@ public:
   UShort_t Counts()                const;
   /** @return Energy deposited */
   Float_t  Edep()                  const { return fEdep;     }
+  /** @return Number of particles that hit this strip */
+  UShort_t NParticles() const { return fNParticles; }
+  /** @return Number of primary particles that hit this strip */
+  UShort_t NPrimaries() const { return fNPrimaries; }
+     
   /** Print info 
       @param opt Not used */
   void     Print(Option_t* opt="") const;
@@ -66,7 +73,9 @@ protected:
   Short_t  fCount2;     // Digital signal (-1 if not used)
   Short_t  fCount3;     // Digital signal (-1 if not used)
   Short_t  fCount4;     // Digital signal (-1 if not used)
-  ClassDef(AliFMDSDigit,2)     // Summable FMD digit
+  UShort_t fNParticles; // Total number of particles that hit this strip
+  UShort_t fNPrimaries; // Number of primary particles that his this strip
+  ClassDef(AliFMDSDigit,3)     // Summable FMD digit
 };
   
 inline UShort_t 
index 7e289571c795d8e82f5cef6817bf28891e8b5314..6cb9204f51da3a9ed794139d147d63ed4dff8c8a 100644 (file)
Binary files a/FMD/Calib/AltroMap/Run0_999999_v0_s0.root and b/FMD/Calib/AltroMap/Run0_999999_v0_s0.root differ
index 665e2223a7535857a67944c4c468dd3e759dbf9b..58732008586641957d977c76b870d1bd5c1e93b8 100644 (file)
Binary files a/FMD/Calib/Dead/Run0_999999_v0_s0.root and b/FMD/Calib/Dead/Run0_999999_v0_s0.root differ
index e55166028f35254ebfef29d9f7fa94b46427cd62..09d6c5bedae77b5b24ba2b1b795afd1900a80984 100644 (file)
Binary files a/FMD/Calib/Pedestal/Run0_999999_v0_s0.root and b/FMD/Calib/Pedestal/Run0_999999_v0_s0.root differ
index f02337ad2b5cbfc74e6b0478c02f33ada0ce791c..e2a4a43d15c6d58f8241c3337f7b6a7fb9cc20be 100644 (file)
Binary files a/FMD/Calib/PulseGain/Run0_999999_v0_s0.root and b/FMD/Calib/PulseGain/Run0_999999_v0_s0.root differ
index 43f10a914cb617f39f36146a4033cdabe88d61df..c1baf976cf7f0e2758a40d22e13dd894332231a2 100644 (file)
Binary files a/FMD/Calib/SampleRate/Run0_999999_v0_s0.root and b/FMD/Calib/SampleRate/Run0_999999_v0_s0.root differ
index 8b68d63b154616d7be9fe7a7451c568ebda42734..598aa051e9187720ca4fd1d5a2af94fe444b3363 100644 (file)
Binary files a/FMD/Calib/StripRange/Run0_999999_v0_s0.root and b/FMD/Calib/StripRange/Run0_999999_v0_s0.root differ
index 5d81a207e6fae6b932cf4223d9d9220bbf77bd01..4bc3f9516199faaa31bb03b2cf899c115b8d5ef5 100644 (file)
Binary files a/FMD/Calib/ZeroSuppression/Run0_999999_v0_s0.root and b/FMD/Calib/ZeroSuppression/Run0_999999_v0_s0.root differ
index f4352c6f2178eb65a50902778e1a98830d031fe5..a5bc57bb915709086b83cb6bf0b24a03e3c208cb 100644 (file)
@@ -184,8 +184,8 @@ Config()
   // EG_t  eg   = kParam_fmd;
   // EG_t  eg   = kFMDFlat; // kParam_2000; // kPythia;
   // EG_t  eg   = kFMDFlat;
-  // EG_t  eg   = kPythia6;
-  EG_t  eg   = kFMD2Flat;
+  EG_t  eg   = kPythia6;
+  // EG_t  eg   = kFMD2Flat;
   Geo_t geo  = kNoHoles;
   Rad_t rad  = kGluonRadiation;
   Mag_t mag  = k5kG;
@@ -345,7 +345,7 @@ Config()
   gMC->SetProcess("MUNU",1);
   gMC->SetProcess("CKOV",1);
   gMC->SetProcess("HADR",1);
-  gMC->SetProcess("LOSS",2); // 0:none 1,3:dray 2:nodray 4:nofluct (def:2)
+  gMC->SetProcess("LOSS",1); // 0:none 1,3:dray 2:nodray 4:nofluct (def:2)
   gMC->SetProcess("MULS",1);
   gMC->SetProcess("RAYL",1);
 
@@ -419,12 +419,12 @@ Config()
   Bool_t useFMD   = kTRUE; 
   Bool_t useFRAME = kFALSE; 
   Bool_t useHALL  = kFALSE; 
-  Bool_t useITS   = kTRUE;
+  Bool_t useITS   = kFALSE;
   Bool_t useMAG   = kFALSE; 
   Bool_t useMUON  = kFALSE; 
   Bool_t usePHOS  = kFALSE; 
   Bool_t usePIPE  = kTRUE; 
-  Bool_t usePMD   = kTRUE; 
+  Bool_t usePMD   = kFALSE; 
   Bool_t useHMPID = kFALSE; 
   Bool_t useSHIL  = kFALSE; 
   Bool_t useT0    = kFALSE; 
@@ -433,7 +433,7 @@ Config()
   Bool_t useTRD   = kFALSE; 
   Bool_t useZDC   = kFALSE; 
   Bool_t useEMCAL = kFALSE; 
-  Bool_t useVZERO = kTRUE;
+  Bool_t useVZERO = kFALSE;
 
   cout << "\t* Creating the detectors ..." << endl;
   // ================= Alice BODY parameters =========================
index 61046af4bc20616598b6a4296a15e756593b5499..da627ef8cdad80a03d72d42f7a18bbb92cf65ec6 100644 (file)
@@ -24,10 +24,10 @@ void
 Simulate(Int_t n=1)
 {
   AliSimulation sim;
-  // AliLog::SetModuleDebugLevel("FMD", 1);
+  AliLog::SetModuleDebugLevel("FMD", 2);
   sim.SetConfigFile("$(ALICE_ROOT)/FMD/Config.C");
-  // sim.SetMakeSDigits("FMD");
-  sim.SetMakeDigits("FMD ITS VZERO PMD"); 
+  sim.SetMakeSDigits("FMD");
+  sim.SetMakeDigits("FMD"); 
   sim.SetWriteRawData("FMD"); 
   // sim.SetMakeDigitsFromHits("FMD"); 
   sim.SetRunQA(":");
index 3b6a5c3627a9def06b87908a5189fa796170bf91..588fde9e966bdcf3fc09dabc181baf4d892838ab 100644 (file)
@@ -24,6 +24,9 @@ Compile(const char* script, Option_t* option="g")
     std::cerr << "No script to compile!" << std::endl;
     return kFALSE;
   }
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libFMDanalysis.so");
   gSystem->Load("libFMDutil.so");
   gSystem->Load("libFMDflow.so");
   TString macroPath(gROOT->GetMacroPath());
diff --git a/FMD/scripts/DrawSDigits.C b/FMD/scripts/DrawSDigits.C
new file mode 100644 (file)
index 0000000..53d404a
--- /dev/null
@@ -0,0 +1,205 @@
+//____________________________________________________________________
+//
+// $Id: DrawSDigits.C 20907 2007-09-25 08:44:03Z cholm $
+//
+// Script that contains a class to draw eloss from hits, versus ADC
+// counts from digits, using the AliFMDInputHits class in the util library. 
+//
+// It draws the energy loss versus the p/(mq^2).  It can be overlayed
+// with the Bethe-Bloc curve to show how the simulation behaves
+// relative to the expected. 
+//
+// Use the script `Compile.C' to compile this class using ACLic. 
+//
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TProfile2D.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <AliFMDSDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDEdepMap.h>
+#include <AliFMDGeometry.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TLatex.h>
+#include <TArrayF.h>
+#include <AliLog.h>
+
+/** @class DrawSDigits
+    @brief Draw hit energy loss versus digit ADC
+    @code 
+    Root> .L Compile.C
+    Root> Compile("DrawSDigits.C")
+    Root> DrawSDigits c
+    Root> c.Run();
+    @endcode
+    @ingroup FMD_script
+ */
+class DrawSDigits : public AliFMDInput
+{
+private:
+  TH1D* fAdc; // Histogram 
+  TProfile2D* fPrimRatio[5]; // Histograms
+public:
+  void Idx2Det(UShort_t idx, UShort_t& d, Char_t& r) const
+  {
+    d = 0;
+    r = '\0';
+    switch (idx) { 
+    case 0: d = 1; r = 'I'; break;
+    case 1: d = 2; r = 'I'; break;
+    case 2: d = 2; r = 'O'; break;
+    case 3: d = 3; r = 'I'; break;
+    case 4: d = 3; r = 'O'; break;
+    }
+  }
+  Short_t Det2Idx(UShort_t d, Char_t r) const
+  { 
+    Short_t idx = -1;
+    switch (d) { 
+    case 1: idx = 0; break;
+    case 2: idx = 1; break;
+    case 3: idx = 3; break;
+    }
+    return (idx + ((r == 'O' || r == 'o') ? 1 : 0));
+  }
+  
+    
+  //__________________________________________________________________
+  DrawSDigits(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1023.5) 
+  { 
+    AddLoad(kSDigits);
+    AddLoad(kGeometry);
+    fAdc = new TH1D("adc", "ADC", m, amin, amax);
+    fAdc->SetXTitle("ADC value");
+    fAdc->Sumw2();
+
+    Int_t   eN   = 50;
+    Float_t eMin = -4;
+    Float_t eMax = 6;
+    Int_t   pN   = 40;
+    Float_t pMin = -.5;
+    Float_t pMax = 359.5;
+    for (UShort_t i = 0; i < 5; i++) { 
+      UShort_t d;
+      Char_t   r;
+      Idx2Det(i, d, r);
+      fPrimRatio[i] = new TProfile2D(Form("primRatio_fmd%d%c", d, r), 
+                                    Form("Primary/Total - FMD%d%c", d, r), 
+                                    eN, eMin, eMax, pN, pMin, pMax);
+      fPrimRatio[i]->SetXTitle("#eta");
+      fPrimRatio[i]->SetYTitle("#varphi");
+      fPrimRatio[i]->SetZTitle("M_{ch,primary}/M_{ch,total}");
+      fPrimRatio[i]->GetXaxis()->SetTitleFont(132);
+      fPrimRatio[i]->GetYaxis()->SetTitleFont(132);
+      fPrimRatio[i]->GetZaxis()->SetTitleFont(132);
+      fPrimRatio[i]->GetXaxis()->SetLabelFont(132);
+      fPrimRatio[i]->GetYaxis()->SetLabelFont(132);
+      fPrimRatio[i]->GetZaxis()->SetLabelFont(132);
+      fPrimRatio[i]->GetXaxis()->SetTitleSize(.06);
+      fPrimRatio[i]->GetYaxis()->SetTitleSize(.06);
+      fPrimRatio[i]->GetZaxis()->SetTitleSize(.06);
+      fPrimRatio[i]->GetXaxis()->SetLabelSize(.06);
+      fPrimRatio[i]->GetYaxis()->SetLabelSize(.06);
+      fPrimRatio[i]->GetZaxis()->SetLabelSize(.06);
+      fPrimRatio[i]->GetXaxis()->SetNdivisions(10);
+      fPrimRatio[i]->GetYaxis()->SetNdivisions(10);
+      fPrimRatio[i]->GetZaxis()->SetNdivisions(10);
+    }
+  }
+  Bool_t Init() 
+  {
+    Bool_t ret = AliFMDInput::Init();
+    AliFMDGeometry::Instance()->Init();
+    AliFMDGeometry::Instance()->InitTransformations();
+    return ret;
+  }
+    
+  //__________________________________________________________________
+  Bool_t ProcessSDigit(AliFMDSDigit* digit)
+  {
+    if (!digit) return kTRUE;
+    fAdc->Fill(digit->Counts());
+    if (digit->NParticles() == 0) return kTRUE;
+    
+
+    Short_t primIdx = Det2Idx(digit->Detector(), digit->Ring());
+    if (primIdx < 0) return kTRUE;
+    
+    AliFMDGeometry* geom = AliFMDGeometry::Instance();
+    Double_t x, y, z;
+    geom->Detector2XYZ(digit->Detector(), 
+                      digit->Ring(), 
+                      digit->Sector(), 
+                      digit->Strip(), 
+                      x, y, z);
+    // We should get the primary vertex and do 
+    // z     += fCurrentVertex;
+    Double_t phi   =  TMath::ATan2(y, x) * 180. / TMath::Pi();
+    Double_t r     =  TMath::Sqrt(y * y + x * x);
+    Double_t theta =  TMath::ATan2(r, z);
+    Double_t eta   = -TMath::Log(TMath::Tan(theta / 2));
+    Double_t ratio = digit->NPrimaries() / digit->NParticles();
+    if (phi < 0) phi += 360;
+    fPrimRatio[primIdx]->Fill(eta, phi, ratio);
+    
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  Bool_t Finish()
+  {
+    gStyle->SetPalette(1);
+    gStyle->SetOptTitle(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
+    
+    TCanvas* c1 = new TCanvas("fmdSdigitSpectra", "FMD SDigit spectra");
+    fAdc->SetStats(kFALSE);
+    if (fAdc->GetEntries() != 0) 
+      fAdc->Scale(1. / fAdc->GetEntries());
+    fAdc->Draw();
+    c1->Modified();
+    c1->Update();
+    c1->cd();
+
+    TCanvas* c2 = new TCanvas("fmdSdigitPrims", "FMD SDigit # prim/# ntotal",
+                             800, 800);
+    c2->SetTopMargin(0);
+    c2->SetRightMargin(0);
+    c2->Divide(2, 3, 0, 0);
+    for (UShort_t i = 0; i < 5; i++) { 
+      Int_t        np = i+1+(i == 0 ? 0 : 1);
+      TVirtualPad* p   = c2->cd(np);
+      p->SetFillColor(0);
+      p->SetTopMargin((np <= 2) ? 0.05 : 0);
+      p->SetLeftMargin((np % 2 == 1) ? 0.20 : 0);
+      p->SetRightMargin((np % 2 == 1) ? 0 : 0.20);
+      p->SetBottomMargin((np >= 5) ? 0.15 : 0);
+      fPrimRatio[i]->SetStats(kFALSE);
+      fPrimRatio[i]->Draw(Form("COL%c", (np % 2 == 1) ? ' ' : 'Z'));
+      UShort_t d;
+      Char_t   r;
+      Idx2Det(i, d, r);
+      TLatex* text = new TLatex(0, 180, Form("FMD%d%c", d, r));
+      text->SetTextFont(132);
+      text->SetTextAlign(22);
+      text->SetTextSize(0.08);
+      text->Draw();
+    }
+    c2->Modified();
+    c2->Update();
+    c2->cd();
+
+    return kTRUE;
+  }
+
+  ClassDef(DrawSDigits,0);
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
index ea37f4ac045030016984d1ce678591013cd72a4e..3a13ae18cec4fc89f666151c650541146fa10e61 100644 (file)
@@ -1,6 +1,9 @@
 void
 Hits2SDigits()
 {
+  AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
+  AliCDBManager::Instance()->SetRun(0);
+  
   AliRunLoader* runLoader = AliRunLoader::Open("galice.root", "Alice", "read");
   if (!runLoader) {
     AliError("Coulnd't read the file galice.root");
@@ -29,8 +32,7 @@ Hits2SDigits()
   }
   TTree* treeE = runLoader->TreeE();
   
-  AliCDBManager::Instance()->SetRun(0);
-  AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
-  
+  AliLog::SetModuleDebugLevel("FMD", 1);
+
   fmd->Hits2SDigits();
 }
index 6611f805616aa3132c89952fc705627169ad358c..c5dba42c5ccbb6927191f8a2ecf75d1e4db15549 100644 (file)
@@ -24,6 +24,9 @@ MakeCalibration(const char* base="local://$ALICE_ROOT")
   AliCDBManager* cdb   = AliCDBManager::Instance();
   cdb->SetDefaultStorage(base);
 
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libFMDanalysis.so");
   gSystem->Load("libFMDutil.so");
   AliFMDCalibFaker f(AliFMDCalibFaker::kAll, 0);
   f.SetRunRange(0,999999);
index 6980bcc0cfe0470954987e996e35b8e8eab5ad5f..8bc2ac2b8b1c59c0d0731e55ee9829da9b45c184 100644 (file)
@@ -15,6 +15,9 @@ PatternSDigits()
   cdb->SetDefaultStorage("local://$ALICE_ROOT");
   cdb->SetRun(0);
   AliFMDParameters::Instance()->Init();
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libFMDanalysis.so");
   gSystem->Load("libFMDutil.so");
   AliFMDPattern* d = new AliFMDPattern;
   d->SetName("sdigit");