Skip to content
Snippets Groups Projects
THcCherenkov.cxx 13.4 KiB
Newer Older
  • Learn to ignore specific revisions
  • 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412
    ///////////////////////////////////////////////////////////////////////////////////////
    //                                                                                   //
    // THcCherenkov                                                                      //
    //                                                                                   //
    // Class for an Cherenkov detector consisting of onw pair of PMT's                   //
    //                                                                                   //
    // 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:                                     //
    //                                                                                   //
    // Variable Name     Description                                                     //
    // adc_1             Raw value of adc 1                                              //
    // adc_2             Raw value of adc 2                                              //
    // adc_p_1           Pedestal substracted value of adc 1                             //
    // adc_p_2           Pedestal substracted value of adc 2                             //
    // npe_1             Number of photo electrons of adc 1                              //
    // npe_2             Number of photo electrons of adc 2                              //
    //                                                                                   //
    ///////////////////////////////////////////////////////////////////////////////////////
    
    #include "THcCherenkov.h"
    #include "TClonesArray.h"
    #include "THcSignalHit.h"
    #include "THaEvData.h"
    #include "THaDetMap.h"
    #include "THcDetectorMap.h"
    #include "THcGlobals.h"
    #include "THaCutList.h"
    #include "THcParmList.h"
    #include "THcHitList.h"
    #include "VarDef.h"
    #include "VarType.h"
    #include "THaTrack.h"
    #include "TClonesArray.h"
    #include "TMath.h"
    
    #include "THaTrackProj.h"
    
    #include <cstring>
    #include <cstdio>
    #include <cstdlib>
    #include <iostream>
    
    using namespace std;
    
    using std::cout;
    using std::cin;
    using std::endl;
    
    #include <iomanip>
    using std::setw;
    using std::setprecision;
    
    //_____________________________________________________________________________
    THcCherenkov::THcCherenkov( const char* name, const char* description,
    				  THaApparatus* apparatus ) :
      THaNonTrackingDetector(name,description,apparatus)
    {
      // Normal constructor with name and description
      fPosADCHits = new TClonesArray("THcSignalHit",16);
    
    //  fTrackProj = new TClonesArray( "THaTrackProj", 5 );
    }
    
    //_____________________________________________________________________________
    THcCherenkov::THcCherenkov( ) :
      THaNonTrackingDetector()
    {
      // Constructor
    }
    
    //_____________________________________________________________________________
    THcCherenkov::~THcCherenkov()
    {
      // Destructor
      delete [] fPosGain;
      delete [] fPosPedLimit;
      delete [] fPosPedMean;
    }
    
    //_____________________________________________________________________________
    THaAnalysisObject::EStatus THcCherenkov::Init( const TDatime& date )
    {
      static const char* const here = "Init()";
    
      cout << "THcCherenkov::Init " << GetName() << endl;
    
      // Should probably put this in ReadDatabase as we will know the
      // maximum number of hits after setting up the detector map
      THcHitList::InitHitList(fDetMap, "THcCherenkovHit", 100); // 100 is max hits
    
      EStatus status;
      if( (status = THaNonTrackingDetector::Init( date )) )
        return fStatus=status;
    
      // Will need to determine which apparatus it belongs to and use the
      // appropriate detector ID in the FillMap call
      if( gHcDetectorMap->FillMap(fDetMap, "HCER") < 0 ) {
        Error( Here(here), "Error filling detectormap for %s.", 
    	     "HCER");
          return kInitError;
      }
    
      return fStatus = kOK;
    }
    
    //_____________________________________________________________________________
    Int_t THcCherenkov::ReadDatabase( const TDatime& date )
    {
      // This function is called by THaDetectorBase::Init() once at the beginning
      // of the analysis.
    
      cout << "THcCherenkov::ReadDatabase " << GetName() << endl; // Ahmed
    
      char prefix[2];
    
      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];
      
      DBRequest list[]={
        {"cer_adc_to_npe", fPosGain, kDouble, fNelem},              // Ahmed
        {"cer_ped_limit", fPosPedLimit, kInt, fNelem},              // Ahmed
        {0}
      };
      gHcParms->LoadParmValues((DBRequest*)&list,prefix);
    
      fIsInit = true;
    
      // Create arrays to hold pedestal results
      InitializePedestals();
    
      return kOK;
    }
    
    //_____________________________________________________________________________
    Int_t THcCherenkov::DefineVariables( EMode mode )
    {
      // Initialize global variables for histogramming and tree
    
      cout << "THcCherenkov::DefineVariables called " << GetName() << endl;
    
      if( mode == kDefine && fIsSetup ) return kOK;
      fIsSetup = ( mode == kDefine );
      
      // Register variables in global list
    
      // Do we need to put the number of pos/neg TDC/ADC hits into the variables?
      // No.  They show up in tree as Ndata.H.aero.postdchits for example
    
      RVarDef vars[] = {
        {"adc_1",  "Raw First ADC Amplitudes",    "fA_1"},
        {"adc_2",  "Raw Second ADC Amplitudes",   "fA_2"},
        {"adc_p_1",  "Raw First ADC Amplitudes",  "fA_p_1"},
        {"adc_p_2",  "Raw Second ADC Amplitudes", "fA_p_2"},
        {"npe_1","PEs First Tube", "fNpe_1"},
        {"npe_2","PEs Second Tube","fNpe_2"},
        {"posadchits", "List of Positive ADC hits","fPosADCHits.THcSignalHit.GetPaddleNumber()"},
        { 0 }
      };
    
      return DefineVarsFromList( vars, mode );
    }
    //_____________________________________________________________________________
    inline 
    void THcCherenkov::Clear(Option_t* opt)
    {
      // Clear the hit lists
      fPosADCHits->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;
    
    }
    
    //_____________________________________________________________________________
    Int_t THcCherenkov::Decode( const THaEvData& evdata )
    {
      // Get the Hall C style hitlist (fRawHitList) for this event
      fNhits = THcHitList::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;
      while(ihit < fNhits) {
        THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit);
        
        // ADC positive hit
        if(hit->fADC_pos >  0) {
          THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
          sighit->Set(hit->fCounter, hit->fADC_pos);
        }
    
        ihit++;
      }
      return ihit;
    }
    
    //_____________________________________________________________________________
    Int_t THcCherenkov::ApplyCorrections( void )
    {
      return(0);
    }
    
    //_____________________________________________________________________________
    Int_t THcCherenkov::CoarseProcess( TClonesArray&  ) //tracks
    {
      /* 
         ------------------------------------------------------------------------------------------------------------------
         h_trans_cer.f code:
    
          hcer_num_hits = 0                                      <--- clear event
          do tube=1,hcer_num_mirrors                             
            hcer_npe(tube) = 0.                                  <--- clear event
            hcer_adc(tube) = 0.                                  <--- clear event
          enddo
          hcer_npe_sum = 0.                                      <--- clear event
          do nhit = 1, hcer_tot_hits                             <--- loop over total hits. Very first line of this method
            tube = hcer_tube_num(nhit)                           <--- tube is number of PMT on either side and it is this 
    	                                                          line:  Int_t npmt = hit->fCounter - 1
            hcer_adc(tube) = hcer_raw_adc(nhit) - hcer_ped(tube) <--- This is done above:
                                                                      fA_Pos_p[npmt] = hit->fADC_pos - fPosPedMean[npmt];
                                                                      fA_Neg_p[npmt] = hit->fADC_neg - fNegPedMean[npmt];
            if (hcer_adc(tube) .gt. hcer_width(tube)) then       <--- This needs to convert in hcana
              hcer_num_hits = hcer_num_hits + 1
              hcer_tube_num(hcer_num_hits) = tube
              hcer_npe(tube) = hcer_adc(tube) * hcer_adc_to_npe(tube)
              hcer_npe_sum = hcer_npe_sum + hcer_npe(tube)
            endif
          enddo
         ------------------------------------------------------------------------------------------------------------------
        */
      for(Int_t ihit=0; ihit < fNhits; ihit++) {
        THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit); // nhit = 1, hcer_tot_hits
    
        // Pedestal subtraction and gain adjustment
    
        // An ADC value of less than zero occurs when that particular
        // channel has been sparsified away and has not been read. 
        // The NPE for that tube will be assigned zero by this code.
        // An ADC value of greater than 8192 occurs when the ADC overflows on
        // an input that is too large. Tubes with this characteristic will
        // be assigned NPE = 100.0.
        Int_t npmt = hit->fCounter - 1;                             // tube = hcer_tube_num(nhit)
        // Should probably check that npmt is in range
        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;
          } 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;
          } else if (  hit->fADC_pos > 8000 ) {
    	fNpe_2 = 100.0;
          } else {
    	fNpe_2 = 0.0;
          }
        }
      }
          
      ApplyCorrections();
    
      return 0;
    }
    
    //_____________________________________________________________________________
    Int_t THcCherenkov::FineProcess( TClonesArray& tracks )
    {
    
      return 0;
    }
    
    //_____________________________________________________________________________
    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];
    
      fPosPed = new Double_t [fNelem];
      fPosThresh = 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;
      }
    
    }
    
    //_____________________________________________________________________________
    void THcCherenkov::AccumulatePedestals(TClonesArray* rawhits)
    {
      // Extract data from the hit list, accumulating into arrays for
      // calculating pedestals
    
      Int_t nrawhits = rawhits->GetLast()+1;
    
      Int_t ihit = 0;
      while(ihit < nrawhits) {
        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];
          }
        }
        ihit++;
      }
    
      fNPedestalEvents++;
    
      return;
    }
    
    //_____________________________________________________________________________
    void THcCherenkov::CalculatePedestals( )
    {
      // Use the accumulated pedestal data to calculate pedestals
      // Later add check to see if pedestals have drifted ("Danger Will Robinson!")
      //  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;
    
        // 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];
          }
        }
      }
      //  cout << " " << endl;
      
    }
    void THcCherenkov::Print( const Option_t* opt) const {
      THaNonTrackingDetector::Print(opt);
    
      // Print out the pedestals
    
      cout << endl;
      cout << "Cherenkov Pedestals" << endl;
    
      // Ahmed
      cout << "No.   Pos" << endl;
      for(Int_t i=0; i<fNelem; i++){
        cout << " " << i << "    " << fPosPed[i] << endl;
      }
    
      cout << endl;
    }
    
    ClassImp(THcCherenkov)
    ////////////////////////////////////////////////////////////////////////////////