Example: CalRecon Analysis

Note:  The following code was taken directly from RootTreeAnalysis in the RootAnalysis package stored in Offline's CVS repository. In this context, the variable "rec" is a ReconEvent object which was filled in RootTreeAnalysis::Go( ) method.  All file handling is done within the initialization of the RootTreeAnalysis object and is not shown here.

The complete documentation of the reconRootData library which defines ReconEvent and all reconstruction classes can be found at:

http://confluence.slac.stanford.edu/display/WB/reconRootData

Code Snippet
Description
void RootTreeAnalysis::ReconCal() {
    // Purpose and Method:  Process on CalRecon event
   
    CalRecon *calRec = rec->getCalRecon();
    if (!calRec) return;

 
 
 		 Retrieve the CalRecon pointer from the ReconEvent object.  If it is non-null, continue.
    float totXE = 0.;
   
    TObjArray *xtalRecCol = calRec->getCalXtalRecCol();
    TIter xtalIter(xtalRecCol);
    CalXtalRecData *xtal = 0;
    while (xtal = (CalXtalRecData*)xtalIter.Next()) {
        Double_t xtalEnergy = xtal->getEnergy();
        if (xtalEnergy > 2000) {
                   const CalXtalId id = xtal->getPackedId();
    int lyr = id.getLayer();
    int twr = id.getTower();
    int col = id.getColumn();
    CalRangeRecData* rData = xtal->getRangeRecData(0);
    int range = rData->getRange(0);
    double ph0 = xtal->getEnergySelectedRange(range,0);
      continue;
         }
      totXE += xtalEnergy;

    }


    ((TH1F*)GetObjectPtr("CALXTALTOTE"))->Fill(totXE);

    		Retrieve the full collection of CalXtalRecData

Create a ROOT TIter to iterate over the CalXtalRecData.  The while loop uses this iterator and fills the xtal variable with the current CalXtalRecData pointer.

Retrieve the reconstructed energy associated with this crystal.  If that energy is greater than 2000 MeV, retrieve the CalXtalId and use the CalXtalId member functions to find the layer, tower, and column.

Get the first CalRangeRecData which corresponds to "Best Range", retrieve the range and energy associated with that range.

Maintain running total of the energy.

Fill the CALXTALTOTE histogram.
     TObjArray*  clusCol = calRec->getCalClusterCol();
    Long64_t numClus = clusCol->GetEntries();
    ((TH1F*)GetObjectPtr("CALNUMCLUS"))->Fill(numClus);

    float totE = 0.;
    for (int jc=0;jc<numClus; jc++) {
      CalCluster* c1 = (CalCluster*)clusCol->At(jc);
      float clusterEnergy = c1->getEnergySum();
      totE += clusterEnergy;
      ((TH1F*)GetObjectPtr("CALRECESUM"))->Fill(c1->getEnergySum());
      ((TH1F*)GetObjectPtr("CALRECELEAK"))->Fill(c1->getEnergyLeak());
      ((TH1F*)GetObjectPtr("CALRECECORR"))->Fill(c1->getEnergyCorrected());
    }
    ((TH1F*)GetObjectPtr("CALTOTE"))->Fill(totE);

    if (totE > 0.) {
      for (int ic=0;ic<numClus; ic++) {
CalCluster* c2 = (CalCluster*)clusCol->At(ic);
float clusterEnergy = c2->getEnergySum();
((TH1F*)GetObjectPtr("CALECLUS"))->Fill(clusterEnergy);
        float eFraction = clusterEnergy/totE;
        ((TH1F*)GetObjectPtr("CALEFRAC"))->Fill(eFraction);
      }
    }

 		 Retrieve the full collection of CalCluster objects.

Use a simple for loop to iterate over the CalCluster objects.
     Float_t recArr[2] = {clusterCol->GetEntries(), xtalRecCol->GetEntries()};
    ((TNtuple*)GetObjectPtr("calReconTup"))->Fill(recArr);
 		 Fill a simple two element ntuple.

 

If you desire to download this piece of code - here is a full copy:

void RootTreeAnalysis::ReconCal() {
    // Purpose and Method:  Process on CalRecon event
   
    CalRecon *calRec = rec->getCalRecon();
    if (!calRec) return;

    float totXE = 0.;
   
    TObjArray *xtalRecCol = calRec->getCalXtalRecCol();
    TIter xtalIter(xtalRecCol);
    CalXtalRecData *xtal = 0;
    while (xtal = (CalXtalRecData*)xtalIter.Next()) {
        Double_t xtalEnergy = xtal->getEnergy();
        if (xtalEnergy > 2000) {
            const CalXtalId id = xtal->getPackedId();
    int lyr = id.getLayer();
    int twr = id.getTower();
    int col = id.getColumn();
    CalRangeRecData* rData = xtal->getRangeRecData(0);
    int range = rData->getRange(0);
    double ph0 = xtal->getEnergySelectedRange(range,0);
      continue;
         }
      totXE += xtalEnergy;

    }


    ((TH1F*)GetObjectPtr("CALXTALTOTE"))->Fill(totXE);

   
    TObjArray*  clusCol = calRec->getCalClusterCol();
    Long64_t numClus = clusCol->GetEntries();
    ((TH1F*)GetObjectPtr("CALNUMCLUS"))->Fill(numClus);

    float totE = 0.;
    for (int jc=0;jc<numClus; jc++) {
      CalCluster* c1 = (CalCluster*)clusCol->At(jc);
      float clusterEnergy = c1->getEnergySum();
      totE += clusterEnergy;
      ((TH1F*)GetObjectPtr("CALRECESUM"))->Fill(c1->getEnergySum());
      ((TH1F*)GetObjectPtr("CALRECELEAK"))->Fill(c1->getEnergyLeak());
      ((TH1F*)GetObjectPtr("CALRECECORR"))->Fill(c1->getEnergyCorrected());
    }
    ((TH1F*)GetObjectPtr("CALTOTE"))->Fill(totE);

    if (totE > 0.) {
      for (int ic=0;ic<numClus; ic++) {
CalCluster* c2 = (CalCluster*)clusCol->At(ic);
float clusterEnergy = c2->getEnergySum();
((TH1F*)GetObjectPtr("CALECLUS"))->Fill(clusterEnergy);
        float eFraction = clusterEnergy/totE;
        ((TH1F*)GetObjectPtr("CALEFRAC"))->Fill(eFraction);
      }
    }

    Float_t recArr[2] = {clusCol->GetEntries(), xtalRecCol->GetEntries()};
    ((TNtuple*)GetObjectPtr("calReconTup"))->Fill(recArr);

    return;
}

 

Owned by: Heather Kelly

Last updated by: Chuck Patterson 03/28/2005