---

EDAnalysis.cxx

Go to the documentation of this file.

#include "TSystem.h"
#include "TStyle.h"

#include "bfemDisplay/EDControl.h"
#include "bfemDisplay/EDAnalysis.h"
#include "bfemDisplay/Enum.h"
#include "bfemDisplay/EDHistos.h"


TH1 *hACD[13];
TH1 *hXGT[4];
TH1 *hTKR[26];
TH1 *hTKRtotL[26];
TH1 *hTKRtotR[26];

TH1F *LOGID;
TH1F *NLOGS;
TH1F *BESTRANGE;
TH1F *ELOG;
TH1F *ETOTADCCAL;
TH1F *ETOTCAL;
TH1F *ADC0;
TH1F *ADC1;
TH1F *ADC2;
TH1F *ADC3;
TH2F *ETOTVSNLOGS;


extern EDControl  *m_evtCtrl;   // Pointer to event control.

EDHistos *h;             // Pointer to histos.

ClassImp(EDAnalysis)


EDAnalysis::EDAnalysis()

{

  fTree = m_evtCtrl->getRawTree();

  event = m_evtCtrl->getEvent();

  h = m_evtCtrl->getHistos();

  m_ped = m_evtCtrl->getCALPeds();

  m_rawdata = m_evtCtrl->getRawDataType(); // Get type of input data: real or mc.


  m_StartEvent = 0;

}


void EDAnalysis::StartWithEvent(Int_t event) {

// set a starting event for Go.

  m_StartEvent = event;

}




void EDAnalysis::Rewind()

{

  //   Start input file at beginning

       m_StartEvent = 0;

       return;

} 




void EDAnalysis::Go(Int_t numEvents)

{
  
  TH1F *hist;
  
  // Event loop_______________________________________ 
  

  fTree = m_evtCtrl->getRawTree();

  event = m_evtCtrl->getEvent();

  m_rawdata = m_evtCtrl->getRawDataType(); 

  if ( fTree == 0 ) return;

  if ( event ) event->Clean();
  
  
  Int_t nentries;
    
  nentries = m_evtCtrl->getEntries();
  
  
  printf("\nNum Events in File is: %i\n", nentries);
  
  Int_t curI;
  
  Int_t nMax = TMath::Min(numEvents+m_StartEvent,nentries);
  
  if (m_StartEvent == nentries) {
    
    printf(" all events in file read\n");
    
    return;

  }



  // Refresh histogram pointers_______________________________________



  TFile *histFile = (TFile*)gROOT->GetFile("Histograms.root");


  ADC0 = (TH1F*)histFile->Get("ADC0");

  ADC1 = (TH1F*)histFile->Get("ADC1");

  ADC2 = (TH1F*)histFile->Get("ADC2");

  ADC3 = (TH1F*)histFile->Get("ADC3");

  LOGID = (TH1F*)histFile->Get("LOGID");

  BESTRANGE = (TH1F*)histFile->Get("BESTRANGE");

  ELOG = (TH1F*)histFile->Get("ELOG");

  ETOTVSNLOGS = (TH2F*)histFile->Get("ETOTVSNLOGS");

  NLOGS = (TH1F*)histFile->Get("NLOGS");

  ETOTCAL = (TH1F*)histFile->Get("ETOTCAL");

  ETOTADCCAL = (TH1F*)histFile->Get("ETOTADCCAL");

  for (Int_t iacd = 0; iacd < NUM_ACD_HISTS; ++iacd) {

    hACD[iacd] = (TH1*)h->m_histList->At(iacd);

  }

  for (Int_t ixgt = 0; ixgt < NUM_XGT_HISTS; ++ixgt) {

    hXGT[ixgt] = (TH1*)h->m_histList->At(ixgt+13);

  }

  for (Int_t ilayer = 0; ilayer < NUM_TKR_HISTS; ++ilayer) {

    hTKR[ilayer] = (TH1*)h->m_histList->At(ilayer+28);
    
    hTKRtotL[ilayer] = (TH1*)h->m_histList->At(ilayer+54);
    
    hTKRtotR[ilayer] = (TH1*)h->m_histList->At(ilayer+80);
  }

  // Start analysis code _______________________________________

  m_Stop = kFALSE;

  for (Int_t ievent=m_StartEvent; ievent<nMax; ievent++, curI=ievent) {


    gSystem->ProcessEvents();

    if ( m_Stop ) break;

    if (event) event->Clean();

    if ( !( ievent%10 ) ) m_evtCtrl->UpdateGuiEvents( ievent );  // Update slider and event count.

    if ( ievent==100 ) m_evtCtrl->UpdateCanvas();  // Show a few preliminary events in histos.

    if ( !( ievent%1000 ) ) m_evtCtrl->UpdateCanvas();   // Update histos every 1000 events.

    if (ievent %100 == 0) printf(" Event: %i\n ", ievent);

    fTree->GetEvent(ievent);

   
    // ACD

    AcdTile *iTile;

    Int_t kADC = 1;

    if(kADC) {

      Int_t nACD = event->getAcdDigi()->GetEntries();

      Float_t eTotalACD;        
          
      for (Int_t ihit=0; ihit < nACD; ihit++) {   
 
        event->getAcdDigi()->Sort();

        AcdTile *hit = (AcdTile*)event->getAcdDigi()->At(ihit);

        AcdId *tile = (AcdId*) hit->getId();

        UInt_t id = tile->getId();

        if (id >1000) continue;

        iTile = ( AcdTile* )event->getAcdTile(id);

        eTotalACD = (Float_t) iTile->getPulseHeight();

        hist = (TH1F*)hACD[ihit];

        hist->Fill(eTotalACD);

      }

    }


    // XGT


    Int_t kXGT = 1;

    if(kXGT) {

      Int_t nXGT = event->getXgtDigi()->GetEntries();
          
      for (Int_t ihit=0; ihit < nXGT; ihit++) {   
 
        event->getXgtDigi()->Sort();

        AcdTile *hit = (AcdTile*)event->getXgtDigi()->At(ihit);

        AcdId *xgt = (AcdId*) hit->getId();

        Int_t id = xgt->getId();

        Float_t eTotalXGT;

        iTile = ( AcdTile* )event->getXgt(id);

        eTotalXGT = iTile->getPulseHeight();

        hist = (TH1F*)hXGT[ihit];

        hist->Fill(eTotalXGT);

      }

    }




    // Tracker...




    Int_t kTKR = 1;

    if(kTKR) {  

      event->getTkrDigi()->Sort();

      int nTKR = event->getTkrDigi()->GetEntries();

      for (Int_t ilay = 0; ilay < nTKR; ilay++) {

        const TkrLayer* layer = event->getTkrLayer(ilay); 

        // getXY = 0 measures Y ==> strips along X
        
        //      Int_t isY = (Int_t)layer->getXY();

        Int_t LayerNum= layer->getLayerNum();
       
 
        hist = (TH1F*)hTKR[LayerNum];       
            
        TObjArray  *stripList = 0;
            
        stripList = ( ( TkrLayer* )layer )->getStrips();

        Int_t nStrips = stripList->GetEntries();

        for (Int_t iStrip=0; iStrip < nStrips; iStrip++) {

          StripId* strip = (StripId*)(stripList->At(iStrip));

          Int_t stripID = strip->getId();

          if (stripID <=1600) hist->Fill(stripID);

        }

        Int_t totR =  ( ( TkrLayer* )layer )->getToT(0);
         
        Int_t totL =  ( ( TkrLayer* )layer )->getToT(1);

        hist = (TH1F*)hTKRtotL[LayerNum];
            
        if (totL) hist->Fill(totL);

        hist = (TH1F*)hTKRtotR[LayerNum];
            
        if (totR) hist->Fill(totR);

      }
      
    }




    // Calorimeter...




    Int_t kCAL = 1;

    if(kCAL) {  

      Int_t nCAL = event->getCalDigi()->GetEntries();

      Int_t nLOGS = 0;

      Float_t eLog = 0;

      Float_t eTotcal = 0, eTotadc = 0;

      Int_t index;

      //Float_t gain[4] = {0.000065, 0.00022, 0.0037, 0.022};
      Float_t gain[4] = {0.000075, 0.0003, 0.0043, 0.025};       //Eric's gain Sept. 11 2001  - in MeV
         
      
      for (Int_t ihit=0; ihit < nCAL; ihit++) {

        Int_t En = 1;

        CalLog *hitcal = (CalLog*) event->getCalDigi()->At(ihit);       

        LogId *log = (LogId*)hitcal->getLogId();

        Int_t id = log->getId(); // ==> id differs from the old id! range now from 0 to 39

        Int_t lay = log->getLayer();

        Int_t column = log->getColumn();

        Int_t iped = lay*10+column;

        Int_t hitADC;


        ((TH1*)h->m_histList->At(17))->Fill(id); // LOGID
           
           

        for (Int_t end=0; end <= 1; end++) {
             
                     
          Float_t eMax = -1000.; Int_t bestRange = -1;

          for (Int_t adc = 0; adc <= 3; adc++) {

            if (!end) hitADC = 4095 - hitcal->getAdc(hitcal->POS,adc);
            else hitADC = 4095 - hitcal->getAdc(hitcal->NEG,adc);

            if (hitADC < 1) continue;

            if(adc == 0) ((TH1*)h->m_histList->At(23))->Fill(hitADC); // ADC0

            if(adc == 1) ((TH1*)h->m_histList->At(24))->Fill(hitADC); // ADC1

            if(adc == 2) ((TH1*)h->m_histList->At(25))->Fill(hitADC); // ADC2

            if(adc == 3) ((TH1*)h->m_histList->At(26))->Fill(hitADC); // ADC3

            if (hitADC > 3900) continue;

            if (hitADC > eMax) {

              eMax = hitADC;

              bestRange = adc;

            }

          }



          if(!end) index = 2 * bestRange;

          else index = 2 * bestRange + 1;

          // Correct for unexpected ordering in ROOT ped file.
          
          /*      if(index==1) index=2;
          
                  if(index==2) index=1;*/



          ((TH1*)h->m_histList->At(19))->Fill(bestRange); // BESTRANGE

          Float_t eHit;

          if (bestRange > -1) {

            if ( m_rawdata ) {  // real raw data 


            if ( m_ped ) {
               
              TVector *vv = (TVector *)m_ped->At(index);
                               
              Float_t calped = vv->operator()( iped );

              eHit = (eMax-calped)*gain[bestRange];

              if (eHit < 0) En = 0;

            }

            else { eHit = 0; En = 0; }
            }

            else { // mc data

              eHit = eMax*gain[bestRange];
              if (eHit <= 0) En = 0;

            }

            eLog += eHit/2.;

            eTotadc += eMax/2.;       


          }

        }

        if (En) {

          ((TH1*)h->m_histList->At(20))->Fill(eLog); // ELOG

          eTotcal += eLog;

          nLOGS++;

        }

        eLog = 0;

      } // End CAL logs loop

         

      if (eTotcal > 0) ((TH1*)h->m_histList->At(22))->Fill(eTotcal); // ETOTCAL
      
      ((TH1*)h->m_histList->At(21))->Fill(eTotadc); // ETOTADCCAL

      ((TH1*)h->m_histList->At(18))->Fill(nLOGS); // NLOGS

      ((TH1*)h->m_histList->At(27))->Fill(eTotadc,(Float_t)nLOGS); // ETOTVSNLOGS

              

    }


  } // End event loop

  // end analysis code in event loop ___________________________


  histFile->Write();

  if ( !( curI%10 ) )m_evtCtrl->UpdateGuiEvents( curI );

  m_StartEvent = curI;

      

  cout << "\nAnalysis finished at event " << curI << ".\n";


} // End GO()

---