From 9ec9143ff8373e8bd268b0533c081f25ae3c0572 Mon Sep 17 00:00:00 2001
From: Zafar <ahmed24z@uregina.ca>
Date: Wed, 29 Jan 2014 19:37:04 -0500
Subject: [PATCH] Updated Cherenkov detector code We have added a new parameter
 hcer_tot_pmts = 2 in hcana.param file. It is equal to total number of PMTs in
 Cherenkov detector. With this update, the THcCherenkov coarse process is
 complete. We have Raw adc, pedestal subtracted adc, number of photo electrons
 and hits histograms in the output root file for gas Cherenkov detector.

---
 examples/PARAM/hcana.param |   4 +
 src/THcCherenkov.cxx       | 194 ++++++++++++++++---------------------
 src/THcCherenkov.h         |  81 ++++------------
 3 files changed, 109 insertions(+), 170 deletions(-)

diff --git a/examples/PARAM/hcana.param b/examples/PARAM/hcana.param
index 47b4766..e437cdf 100644
--- a/examples/PARAM/hcana.param
+++ b/examples/PARAM/hcana.param
@@ -77,3 +77,7 @@ sdc_fix_lr = 0
 # If 1, don't do the the propagation along the wire each time the hit
 # appears in a space point.  (Which means the correction accumulates)
 sdc_fix_propcorr = 0
+
+# Total number of PMTs in Gas Cherenkov detector.
+hcer_tot_pmts = 2
+
diff --git a/src/THcCherenkov.cxx b/src/THcCherenkov.cxx
index 5587a54..fe32f7e 100644
--- a/src/THcCherenkov.cxx
+++ b/src/THcCherenkov.cxx
@@ -4,14 +4,7 @@
 //                                                                                   //
 // Class for an Cherenkov detector consisting of onw pair of PMT's                   //
 //                                                                                   //
-// Zafar Ahmed. Updated on December 24 2013.                                         //
-// Four more variables are added.                                                    //
-//                                                                                   //
-// npe               Total Number of photo electrons                                 //
-// hit_1             Total hits in adc 1                                             //
-// hit_2             Total hits in adc 2                                             //
-// hit               Total hits in adc 1 and 2                                       //
-//                                                                                   //
+// Zafar Ahmed. Second attempt. November 14 2013.                                    //
 // Comment:No need to cahnge the Map file but may need to change the parameter file  //
 //                                                                                   //
 // This code is written for following variables:                                     //
@@ -65,9 +58,8 @@ THcCherenkov::THcCherenkov( const char* name, const char* description,
   THaNonTrackingDetector(name,description,apparatus)
 {
   // Normal constructor with name and description
-  fPosADCHits = new TClonesArray("THcSignalHit",16);
+  fADCHits = new TClonesArray("THcSignalHit",16);
 
-//  fTrackProj = new TClonesArray( "THaTrackProj", 5 );
 }
 
 //_____________________________________________________________________________
@@ -81,9 +73,15 @@ THcCherenkov::THcCherenkov( ) :
 THcCherenkov::~THcCherenkov()
 {
   // Destructor
-  delete [] fPosGain;
-  delete [] fPosPedLimit;
-  delete [] fPosPedMean;
+  delete [] fNPMT;
+  delete [] fADC;
+  delete [] fADC_P;
+  delete [] fNPE;
+
+  delete [] fCerWidth;
+  delete [] fGain;
+  delete [] fPedLimit;
+  delete [] fPedMean;
 }
 
 //_____________________________________________________________________________
@@ -121,21 +119,34 @@ Int_t THcCherenkov::ReadDatabase( const TDatime& date )
   cout << "THcCherenkov::ReadDatabase " << GetName() << endl; // Ahmed
 
   char prefix[2];
+  char parname[100];
 
   prefix[0]=tolower(GetApparatus()->GetName()[0]);
   prefix[1]='\0';
 
-  fNelem = 2;      // Default if not defined                                                                    
-  
-  fPosGain = new Double_t[fNelem];
-  fPosPedLimit = new Int_t[fNelem];
-  fPosPedMean = new Double_t[fNelem];
+  strcpy(parname,prefix);                              // This is taken from 
+  strcat(parname,"cer_tot_pmts");                      // THcScintillatorPlane
+  fNelem = (Int_t)gHcParms->Find(parname)->GetValue(); // class.
+
+  //    fNelem = 2;      // Default if not defined                                                                    
+
+  fNPMT = new Int_t[fNelem];
+  fADC = new Double_t[fNelem];
+  fADC_P = new Double_t[fNelem];
+  fNPE = new Double_t[fNelem];
+
+  fCerWidth = new Double_t[fNelem];
+  fGain = new Double_t[fNelem];
+  fPedLimit = new Int_t[fNelem];
+  fPedMean = new Double_t[fNelem];
   
   DBRequest list[]={
-    {"cer_adc_to_npe", fPosGain, kDouble, fNelem},              // Ahmed
-    {"cer_ped_limit", fPosPedLimit, kInt, fNelem},              // Ahmed
+    {"cer_adc_to_npe", fGain,     kDouble, fNelem},              // Ahmed
+    {"cer_ped_limit",  fPedLimit, kInt,    fNelem},              // Ahmed
+    {"cer_width",      fCerWidth, kDouble, fNelem},              // Ahmed
     {0}
   };
+
   gHcParms->LoadParmValues((DBRequest*)&list,prefix);
 
   fIsInit = true;
@@ -162,17 +173,12 @@ Int_t THcCherenkov::DefineVariables( EMode mode )
   // No.  They show up in tree as Ndata.H.aero.postdchits for example
 
   RVarDef vars[] = {
-    {"adc_1",    "Raw First ADC Amplitude",                  "fA_1"},
-    {"adc_2",    "Raw Second ADC Amplitude",                 "fA_2"},
-    {"adc_p_1",  "Pedestal Subtracted First ADC Amplitude",  "fA_p_1"},
-    {"adc_p_2",  "Pedestal Subtracted Second ADC Amplitude", "fA_p_2"},
-    {"npe_1",    "PEs of First Tube",                        "fNpe_1"},
-    {"npe_2",    "PEs of Second Tube",                       "fNpe_2"},
-    {"npe",      "Total number of PEs",                      "fNpe"},
-    {"hit_1",    "ADC hits First Tube",                      "fNHits_1"},
-    {"hit_2",    "ADC hits Second Tube",                     "fNHits_2"},
-    {"hit",      "Total ADC hits",                           "fNHits"},
-    {"posadchits", "List of Positive ADC hits","fPosADCHits.THcSignalHit.GetPaddleNumber()"},
+    {"PhotoTubes",  "Nuber of Cherenkov photo tubes",            "fNPMT"},
+    {"ADC",         "Raw ADC values",                            "fADC"},
+    {"ADC_P",       "Pedestal Subtracted ADC values",            "fADC_P"},
+    {"NPE",         "Number of Photo electrons",                 "fNPE"},
+    {"NPEsum",      "Sum of Number of Photo electrons",          "fNPEsum"},
+    {"NCherHit",    "Number of Hits(Cherenkov)",                 "fNCherHit"},
     { 0 }
   };
 
@@ -183,22 +189,20 @@ inline
 void THcCherenkov::Clear(Option_t* opt)
 {
   // Clear the hit lists
-  fPosADCHits->Clear();
+  fADCHits->Clear();
 
   // Clear Cherenkov variables  from h_trans_cer.f
   
   fNhits = 0;	     // Don't really need to do this.  (Be sure this called before Decode)
-
-  fA_1 = 0;
-  fA_2 = 0;
-  fA_p_1 = 0;
-  fA_p_2 = 0;
-  fNpe_1 = 0;
-  fNpe_2 = 0;
-  fNpe = 0;
-  fNHits_1 = 0;
-  fNHits_2 = 0;
-  fNHits = 0;
+  fNPEsum = 0;
+  fNCherHit = 0;
+
+  for(Int_t itube = 0;itube < fNelem;itube++) {
+    fNPMT[itube] = 0;
+    fADC[itube] = 0;
+    fADC_P[itube] = 0;
+    fNPE[itube] = 0;
+  }
 
 }
 
@@ -209,29 +213,25 @@ Int_t THcCherenkov::Decode( const THaEvData& evdata )
   fNhits = DecodeToHitList(evdata);
 
   if(gHaCuts->Result("Pedestal_event")) {
-
     AccumulatePedestals(fRawHitList);
-
     fAnalyzePedestals = 1;	// Analyze pedestals first normal events
     return(0);
   }
 
   if(fAnalyzePedestals) {
-     
     CalculatePedestals();
-   
     fAnalyzePedestals = 0;	// Don't analyze pedestals next event
   }
 
 
   Int_t ihit = 0;
-  Int_t nPosADCHits=0;
+  Int_t nADCHits=0;
   while(ihit < fNhits) {
     THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit);
     
-    // ADC positive hit
+    // ADC hit
     if(hit->fADC_pos >  0) {
-      THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+      THcSignalHit *sighit = (THcSignalHit*) fADCHits->ConstructedAt(nADCHits++);
       sighit->Set(hit->fCounter, hit->fADC_pos);
     }
 
@@ -274,7 +274,6 @@ Int_t THcCherenkov::CoarseProcess( TClonesArray&  ) //tracks
       enddo
      ------------------------------------------------------------------------------------------------------------------
     */
-
   for(Int_t ihit=0; ihit < fNhits; ihit++) {
     THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit); // nhit = 1, hcer_tot_hits
 
@@ -291,42 +290,21 @@ Int_t THcCherenkov::CoarseProcess( TClonesArray&  ) //tracks
     if ( ihit != npmt )
       cout << "ihit != npmt " << endl;
 
-    if ( npmt == 0 ) {
-      fA_1 = hit->fADC_pos;
-      fA_p_1 = hit->fADC_pos - fPosPedMean[npmt];
-      if ( ( fA_p_1 > 50 ) && ( hit->fADC_pos < 8000 ) ) {
-	fNpe_1 = fPosGain[npmt]*fA_p_1;
-	fNHits_1 ++;
-      } else if (  hit->fADC_pos > 8000 ) {
-	fNpe_1 = 100.0;
-      } else {
-	fNpe_1 = 0.0;
-      }
-    }
-    
-    if ( npmt == 1 ) {
-      fA_2 = hit->fADC_pos;
-      fA_p_2 = hit->fADC_pos - fPosPedMean[npmt];
-      if ( ( fA_p_2 > 50 ) && ( hit->fADC_pos < 8000 ) ) {
-	fNpe_2 = fPosGain[npmt]*fA_p_2;
-	fNHits_2 ++;
-      } else if (  hit->fADC_pos > 8000 ) {
-	fNpe_2 = 100.0;
-      } else {
-	fNpe_2 = 0.0;
-      }
-    }
-    
-    if ( npmt == 0 ) {
-      fNpe += fNpe_1;
-      fNHits += fNHits_1;
-    }
+    fNPMT[npmt] = hit->fCounter;
+    fADC[npmt] = hit->fADC_pos;
+    fADC_P[npmt] = hit->fADC_pos - fPedMean[npmt];
     
-    if ( npmt == 1 ) {
-      fNpe += fNpe_2;
-      fNHits += fNHits_2;
+    if ( ( fADC_P[npmt] > fCerWidth[npmt] ) && ( hit->fADC_pos < 8000 ) ) {
+      fNPE[npmt] = fGain[npmt]*fADC_P[npmt];
+      fNCherHit ++;
+    } else if (  hit->fADC_pos > 8000 ) {
+      fNPE[npmt] = 100.0;
+    } else {
+      fNPE[npmt] = 0.0;
     }
     
+    fNPEsum += fNPE[npmt];
+
   }
 
   ApplyCorrections();
@@ -346,18 +324,16 @@ void THcCherenkov::InitializePedestals( )
 {
   fNPedestalEvents = 0;
   fMinPeds = 500; 		// In engine, this is set in parameter file
-  fPosPedSum = new Int_t [fNelem];
-  fPosPedSum2 = new Int_t [fNelem];
-  fPosPedLimit = new Int_t [fNelem];
-  fPosPedCount = new Int_t [fNelem];
+  fPedSum = new Int_t [fNelem];
+  fPedSum2 = new Int_t [fNelem];
+  fPedCount = new Int_t [fNelem];
 
-  fPosPed = new Double_t [fNelem];
-  fPosThresh = new Double_t [fNelem];
+  fPed = new Double_t [fNelem];
+  fThresh = new Double_t [fNelem];
   for(Int_t i=0;i<fNelem;i++) {
-    fPosPedSum[i] = 0;
-    fPosPedSum2[i] = 0;
-    fPosPedLimit[i] = 1000;	// In engine, this are set in parameter file
-    fPosPedCount[i] = 0;
+    fPedSum[i] = 0;
+    fPedSum2[i] = 0;
+    fPedCount[i] = 0;
   }
 
 }
@@ -375,13 +351,13 @@ void THcCherenkov::AccumulatePedestals(TClonesArray* rawhits)
     THcCherenkovHit* hit = (THcCherenkovHit *) rawhits->At(ihit);
 
     Int_t element = hit->fCounter - 1;
-    Int_t adcpos = hit->fADC_pos;
-    if(adcpos <= fPosPedLimit[element]) {
-      fPosPedSum[element] += adcpos;
-      fPosPedSum2[element] += adcpos*adcpos;
-      fPosPedCount[element]++;
-      if(fPosPedCount[element] == fMinPeds/5) {
-	fPosPedLimit[element] = 100 + fPosPedSum[element]/fPosPedCount[element];
+    Int_t nadc = hit->fADC_pos;
+    if(nadc <= fPedLimit[element]) {
+      fPedSum[element] += nadc;
+      fPedSum2[element] += nadc*nadc;
+      fPedCount[element]++;
+      if(fPedCount[element] == fMinPeds/5) {
+	fPedLimit[element] = 100 + fPedSum[element]/fPedCount[element];
       }
     }
     ihit++;
@@ -400,19 +376,17 @@ void THcCherenkov::CalculatePedestals( )
   //  cout << "Plane: " << fPlaneNum << endl;
   for(Int_t i=0; i<fNelem;i++) {
     
-    // Positive tubes
-    fPosPed[i] = ((Double_t) fPosPedSum[i]) / TMath::Max(1, fPosPedCount[i]);
-    fPosThresh[i] = fPosPed[i] + 15;
-
-    //    cout << i+1 << " " << fPosPed[i] << " " << fNegPed[i] << endl;
+    // PMT tubes
+    fPed[i] = ((Double_t) fPedSum[i]) / TMath::Max(1, fPedCount[i]);
+    fThresh[i] = fPed[i] + 15;
 
     // Just a copy for now, but allow the possibility that fXXXPedMean is set
     // in a parameter file and only overwritten if there is a sufficient number of
     // pedestal events.  (So that pedestals are sensible even if the pedestal events were
     // not acquired.)
     if(fMinPeds > 0) {
-      if(fPosPedCount[i] > fMinPeds) {
-	fPosPedMean[i] = fPosPed[i];
+      if(fPedCount[i] > fMinPeds) {
+	fPedMean[i] = fPed[i];
       }
     }
   }
@@ -428,9 +402,9 @@ void THcCherenkov::Print( const Option_t* opt) const {
   cout << "Cherenkov Pedestals" << endl;
 
   // Ahmed
-  cout << "No.   Pos" << endl;
+  cout << "No.   ADC" << endl;
   for(Int_t i=0; i<fNelem; i++){
-    cout << " " << i << "    " << fPosPed[i] << endl;
+    cout << " " << i << "    " << fPed[i] << endl;
   }
 
   cout << endl;
diff --git a/src/THcCherenkov.h b/src/THcCherenkov.h
index 4040a56..3196321 100644
--- a/src/THcCherenkov.h
+++ b/src/THcCherenkov.h
@@ -36,78 +36,39 @@ class THcCherenkov : public THaNonTrackingDetector, public THcHitList {
 
   THcCherenkov();  // for ROOT I/O		
  protected:
-  Int_t fAnalyzePedestals;
+  Int_t         fAnalyzePedestals;
 
   // Parameters
-  Double_t* fPosGain;
-  Double_t* fNegGain;
+  Double_t*     fGain;
+  Double_t*     fCerWidth;
 
   // Event information
-  Int_t fNhits;
-
-  Float_t*   fA_Pos;         // [fNelem] Array of ADC amplitudes
-  Float_t*   fA_Neg;         // [fNelem] Array of ADC amplitudes
-  Float_t*   fA_Pos_p;	     // [fNelem] Array of ped-subtracted ADC amplitudes
-  Float_t*   fA_Neg_p;	     // [fNelem] Array of ped-subtracted ADC amplitudes
-  Float_t*   fT_Pos;         // [fNelem] Array of TDCs
-  Float_t*   fT_Neg;         // [fNelem] Array of TDCs
-
-  Double_t  fA_1;         // Ahmed
-  Double_t  fA_2;         // Ahmed
-  Double_t  fNHits_1;     // Ahmed
-  Double_t  fNHits_2;     // Ahmed
-  Double_t  fNHits;       // Ahmed
-  Double_t  fA_p_1;       // Ahmed
-  Double_t  fA_p_2;       // Ahmed
-  Double_t  fNpe_1;       // Ahmed
-  Double_t  fNpe_2;       // Ahmed
-  Double_t  fNpe;         // Ahmed
-
-  Double_t fPosNpeSum;
-  Double_t fNegNpeSum;
-  Double_t fNpeSum;
-  Int_t fNGoodHits;
-  Int_t fNADCPosHits;
-  Int_t fNADCNegHits;
-  Int_t fNTDCPosHits;
-  Int_t fNTDCNegHits;
-
-  Double_t* fPosNpe;		// [fNelem] # Photoelectrons per positive tube
-  Double_t* fNegNpe;		// [fNelem] # Photoelectrons per negative tube
-
+  Int_t         fNhits;
+  Int_t*        fNPMT;            // [fNelem] Array of ADC amplitudes
+  Double_t*     fADC;             // [fNelem] Array of ADC amplitudes
+  Double_t*     fADC_P;           // [fNelem] Array of ADC amplitudes
+  Double_t*     fNPE;             // [fNelem] Array of ADC amplitudes
+  Double_t      fNPEsum;          // [fNelem] Array of ADC amplitudes
+  Double_t      fNCherHit;        // [fNelem] Array of ADC amplitudes
 
   // Hits
-  TClonesArray* fPosTDCHits;
-  TClonesArray* fNegTDCHits;
-  TClonesArray* fPosADCHits;
-  TClonesArray* fNegADCHits;
+  TClonesArray* fADCHits;
 
   // Pedestals
-  Int_t fNPedestalEvents;
-  Int_t fMinPeds;
-  Int_t *fPosPedSum;		/* Accumulators for pedestals */
-  Int_t *fPosPedSum2;
-  Int_t *fPosPedLimit;
-  Int_t *fPosPedCount;
-  Int_t *fNegPedSum;
-  Int_t *fNegPedSum2;
-  Int_t *fNegPedLimit;
-  Int_t *fNegPedCount;
-
-  Double_t *fPosPed;
-  Double_t *fPosSig;
-  Double_t *fPosThresh;
-  Double_t *fNegPed;
-  Double_t *fNegSig;
-  Double_t *fNegThresh;
-
-  Double_t *fPosPedMean; 	/* Can be supplied in parameters and then */
-  Double_t *fNegPedMean;	/* be overwritten from ped analysis */
+  Int_t         fNPedestalEvents;
+  Int_t         fMinPeds;
+  Int_t*        fPedSum;	  /* Accumulators for pedestals */
+  Int_t*        fPedSum2;
+  Int_t*        fPedLimit; 
+  Double_t*     fPedMean; 	  /* Can be supplied in parameters and then */ 
+  Int_t*        fPedCount; 
+  Double_t*     fPed;
+  Double_t*     fThresh;
   
   void Setup(const char* name, const char* description);
   virtual void  InitializePedestals( );
 
-  ClassDef(THcCherenkov,0)   // Generic cherenkov class
+  ClassDef(THcCherenkov,0)        // Generic cherenkov class
 };
 
 #endif
-- 
GitLab