PhotonStatSummarizer.cxx

Go to the documentation of this file.

#include "PhotonStatSummarizer.h"
#include "PulserCalibration/PulserConventions.h"
#include "Calibrator/ValueErr.h"
#include "Calibrator/CalPmtDrift.h"
#include "Calibrator/Calibrator.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "DatabaseInterface/DbiSqlContext.h"
#include "DatabaseInterface/DbiWriter.h"
#include "DatabaseInterface/DbiStatement.h"
#include "DatabaseInterface/DbiTableProxyRegistry.h"
#include "DatabaseInterface/DbiCascader.h"
#include <TFile.h>
#include <TTree.h>
#include <iostream>
#include <TH1.h>
#include "MessageService/MsgService.h"
#include <cmath>

CVSID( " $Id: PhotonStatSummarizer.cxx,v 1.25 2006/05/27 10:05:44 rhatcher Exp $ " );

using namespace std;

ClassImp(PhotonStatSummarizer)
ClassImp(SpotDrift)
ClassImp(PmtDrift)
ClassImp(DetectorDrift)


bool operator<(const SpotDrift& lhs, const SpotDrift& rhs)
{
  // return true if lhs is better than rhs.
  if((lhs.bad>0) != (rhs.bad>0)){
    if(rhs.bad>0) return true;
    else          return false;
  }
  return (lhs.gain.GetError() < rhs.gain.GetError());
}

void PmtDrift::Print(const char*) const
{
  cout << "PmtDrift " << pmt.AsString() << "(" << encoded << ") "
       << rawGain << "  " << meanGain << "  " << nspot << "  " << ngood
       << endl;
}


FloatErr GetWeightedMean(std::vector<SpotDrift>::const_iterator begin,
                         std::vector<SpotDrift>::const_iterator end)
{
  double n=0;
  double b=0; // sum of x_i/sigma_i^2
  double c=0; // sum of 1/sigma_i^2

  std::vector<SpotDrift>::const_iterator d = begin;;
  for(; d!=end; ++d) {
    n+=1.0;
    //double s2 = d->gain.GetError2(); // Does it properly.
    double s2 = 9.92e3/(d->n); // Roughly converts to 4.5% on 70 ADC counts.
    b+=d->gain.GetValue() / s2;
    c+=1.0/s2;
  }

  if(n<1) {
    return FloatErr(0,0);
  }

  return FloatErr((float)(b/c),
                  (float)sqrt(1/c));
                    
}


  
PhotonStatSummarizer::PhotonStatSummarizer()
{
  CreateTables();

  fTask = CalPmtDrift::kPhotonStat;
  if(fWhichEnd==Pulser::kFarEnd) fTask = CalPmtDrift::kPhotonStatFarEnd;

  Registry r;
  r.Set("WriteDB",1);
  r.Set("OverwriteDB",0);
  r.Set("WriteDetTree",1);
  r.Set("WritePmtTree",1);
  r.Set("WriteSpotTree",0);
  r.Set("FileName","drifts.root");
  r.Set("Prescale",0);
  r.Set("WhichEnd",Pulser::kNearEnd); // near 1 far 2

  InitializeConfig(r);

  fDetTree = fPmtTree = fSpotTree = 0;
  fFile = 0;
}

void PhotonStatSummarizer::ConfigModified()
{
  fWriteDB = GetConfig().GetInt("WriteDB");
  fOverwriteDB = GetConfig().GetInt("OverwriteDB");
  fWriteDetTree = GetConfig().GetInt("WriteDetTree");
  fWritePmtTree = GetConfig().GetInt("WritePmtTree");
  fWriteSpotTree = GetConfig().GetInt("WriteSpotTree");
  fPrescale = GetConfig().GetInt("Prescale");
  fWhichEnd = GetConfig().GetInt("WhichEnd");

}


PhotonStatSummarizer::~PhotonStatSummarizer()
{
  if(fDetTree) {
    MSG("PmtDrift",Msg::kInfo) << "Writing Detector tree data." << endl;
    fDetTree->Write("",TObject::kOverwrite);
  }

  if(fPmtTree) {
    MSG("PmtDrift",Msg::kInfo) << "Writing PMT tree data." << endl;
    fPmtTree->Write("",TObject::kOverwrite);
  }

  if(fSpotTree) {
    MSG("PmtDrift",Msg::kInfo) << "Writing PixelSpot tree data." << endl;
    fSpotTree->Write("",TObject::kOverwrite);
  }

  if(fFile) {    
    MSG("PmtDrift",Msg::kInfo) << "Closing file." << endl;  
    fFile->Write();
    delete fFile;
  };
}




void PhotonStatSummarizer::SummarizePoint(
                                          VldContext vldPointStart, // Time at which this data became valid
                                          VldContext vldPoint       // Time at which the tables are all valid.
                                          )
{
  MSG("PmtDrift",Msg::kInfo) << endl;
  MSG("PmtDrift",Msg::kInfo) << "Looking at drift point starting " << vldPointStart.AsString() << endl;

  const int kGoodVersion = 2; // If data exists with this version, don't overwrite
  const int kCurrentVersion = 4; // Write new data with this version.

  if(fPrescale) {
    static int iteration = 0;
    iteration++;
    if((iteration % fPrescale)!=0) {
      MSG("PmtDrift",Msg::kInfo) << "--Prescaling away " << vldPointStart.AsString() << endl;
      return; // skip.
    }
  }

  if(fWriteDB) {
    if(!fOverwriteDB) {
      for(int iVersion=kGoodVersion; iVersion<=kCurrentVersion; iVersion++){
        if(DoesTableExist(iVersion, vldPointStart, fTask)) {
          MSG("PmtDrift",Msg::kInfo) << "Skipping drift point at " << vldPointStart.AsString() 
                                     << " because version " << iVersion << " data already exists." << endl;
          return;
        }
      }
    }
  }

  // Get the Pixel data:
  LoadDataFromDB(vldPoint);

    // Build the PMT info:
  ComputePmtDrifts();

  ComputeDetectorDrift();

  if(fWriteDB) {
    // Fill the database:
    FillPmtDatabase(kCurrentVersion, vldPointStart, fTask);
  };
  
  if(fWritePmtTree||fWriteSpotTree||fWriteDetTree) 
    FillTrees(vldPoint);
  
}


void PhotonStatSummarizer::CreateTables()
{
  // create tables.

  DbiStatement* s = DbiTableProxyRegistry::Instance()
    .GetCascader()
    .CreateStatement(0);
  
  //s->ExecuteUpdate("drop table CALPMTDRIFTVLD");
  //s->ExecuteUpdate("drop table CALPMTDRIFT");


  s->ExecuteUpdate("create table if not exists CALPMTDRIFTVLD ("
                   " SEQNO         int not null primary key auto_increment,"
                   " TIMESTART     datetime not null,"
                   " TIMEEND       datetime not null,"
                   " DETECTORMASK  tinyint,"
                   " SIMMASK       tinyint,"
                   " TASK          int,"
                   " AGGREGATENO   int,"
                   " CREATIONDATE  datetime not null,"
                   " INSERTDATE    datetime not null )");

  s->ExecuteUpdate("create table if not exists CALPMTDRIFT ("
                   " SEQNO integer,"
                   " ROW_COUNTER int, "
                   " PMT bigint, "
                   " CHANNEL bigint,"
                   " DRIFT float, "
                   " STATERROR float,"
                   " SYSERROR float,"
                   " SPOTSUSED float,"
                   " SPOTSAVAIL float,"
                   " CRUDEDRIFT float,"
                   " MINDRIFT float,"
                   " MAXDRIFT float,"
                   " index (SEQNO))");
}



void PhotonStatSummarizer::LoadDataFromDB( VldContext vldPoint )
{
  fDriftTable.NewQuery(vldPoint,fWhichEnd);
  MSG("PmtDrift",Msg::kInfo) << "Loaded " << vldPoint.AsString() << " end = " << fWhichEnd << " rows = " << fDriftTable.GetNumRows() << endl;

  fSpots.clear();

  PlexHandle plex(vldPoint);

  for(UInt_t i=0;i<fDriftTable.GetNumRows();i++) {
    const PulserDrift* row = fDriftTable.GetRow(i);
    if(row) {
      PlexStripEndId seid = 
        PlexStripEndId::UnbuildPlnStripEndKey(vldPoint.GetDetector(),
                                              row->GetStripEnd()
                                              );

      PlexPixelSpotId psid = plex.GetPixelSpotId(seid);
      //cout << "Looking at " << row->GetStripEnd() << " " << seid.AsString() << " " << psid.AsString() << endl;
      //cout << row->GetStripEnd() << " " << row->GetAggregateNo() << " " << row->GetMean() << " " << row->GetError() << endl;

      if(psid.IsValid()) {
        SpotDrift& s = fSpots[psid];
        if(s.n != 0) MSG("PmtDrift",Msg::kWarning) << "Conflict! Two strips with psid " << psid.AsString() << " " << seid.AsString() << endl;
        s.seid = seid;
        s.psid = psid;
        s.rcid = plex.GetRawChannelId(seid);
        s.crate  =  s.rcid.GetCrate();
        s.geo      = s.rcid.GetGeographicAddress();
        s.masterch = s.rcid.GetMasterChannel();
        s.minderch = s.rcid.GetMinderChannel();
        s.varc = s.rcid.GetVarcId();
        s.vmm = s.rcid.GetVmm();
        s.vfb = s.rcid.GetVaAdcSel();
        s.vachip = s.rcid.GetVaChip();
        s.vachan = s.rcid.GetVaChannel();
          
        // Notes on the math:
        // Assuming a gaussian distribution:
        //
        // The variance on the mean is Var(x) = sigma^2
        // Error on the mean is = sigma/sqrt(N)
        // 
        // The variance on sigma^2 is:
        //   Var(sigma^2) = 2 sigma^4
        
        s.n        = row->GetNumEntries();
        s.nflash   = row->GetNumTriggers();
        s.mean     = FloatErr(row->GetMean(),row->GetError());
        float err  = row->GetError();
        s.rms2     = err*err*s.n;                // rms squared
        if(s.n>0)
          s.rms2.SetError(sqrt(2.0)*s.rms2.GetValue()/sqrt(s.n)); // Statisical error on sigma squared
        else 
          s.rms2.SetError(s.rms2.GetValue() * 1e9); // big
                
        const FloatErr k1PeWidthCorr(1.0/1.2,0);
        
        // Includes error calculation:
        s.gain = (s.rms2/s.mean)*k1PeWidthCorr; // Rough correction for 1-pe width.
        
        // Apply zero-correction if the gain is obviously low.
        if(s.mean<150) {          
          FloatErr corrMean = s.mean*s.n/s.nflash;
          FloatErr corrRms2 = s.rms2; //((s.n-1)*s.rms2 + (s.nflash-s.n)*s.mean*s.mean) / FloatErr(s.nflash);
          s.gain = corrRms2/corrMean * k1PeWidthCorr;
        }

        if(psid.GetElecType()==ElecType::kVA) {
          // Appropriate for FD:        
          s.bad = 0; // good by default
          
          //if((s.n/s.nflash)<0.950) s.bad = 1; // Missed more than 5% of flashes           
          if(s.mean<150.)          s.bad = 1; // mean is < 2.5 pe
          if(s.mean>9000.)         s.bad = 2; // Mean is > 150pe
        } else {
          // ND:
          s.bad = 0; // good by default
          
          //if((s.n/s.nflash)<0.950) s.bad = 1; // Missed more than 5% of flashes           
          if(s.mean<250.)           s.bad = 1; // mean is < 2.5 pe
          if(s.mean>15000.)         s.bad = 2; // Mean is > 150pe
          if(s.gain>400.)           s.bad = 3;
        }       
        

      }
    }
  }
  MSG("PmtDrift",Msg::kInfo) << "Found " << fSpots.size() << " spots." << endl;
}




void PhotonStatSummarizer::ComputePmtDrifts()
{
  // Fill the spot data into the PMT objects:
  fPmts.clear();

  SpotList_t::iterator it;
  SpotList_t::iterator end= fSpots.end();

  for(it = fSpots.begin(); it!=end; it++) {
    PlexPixelSpotId pmt = it->first;
    FloatErr spotGain = it->second.gain;

    
    PmtDrift& d = fPmts[pmt.GetUniquePmtEncodedValue()];
    d.pmt = pmt;

    if(d.nspot>127) {
      MSG("PmtDrift",Msg::kError) << "Found too many spots. " << it->first.AsString() << endl;
    } else {
     d.spots[d.nspot] = it->second;
     d.nspot    += 1;
    }
    if(spotGain > d.highestGain) d.highestGain = spotGain;
    if(spotGain < d.lowestGain) d.lowestGain = spotGain;    
    d.rcid = it->second.rcid;
    //cout << "Inserting spot " << pmt.AsString() << " into pmt " << pmt.GetUniquePmtEncodedValue() << " gain " << spotGain << endl;
    //d.Print();
    //fPmts.begin()->second.Print();
  }

  // Compute the PMT statistics.

  PmtList_t::iterator pit;
  PmtList_t::iterator pend = fPmts.end();

  for(pit = fPmts.begin(); pit!=pend; pit++) {
    PmtDrift& d = pit->second;

    PlexPixelSpotId pmt = d.pmt;
    d.pmt.SetPixel(127);
    d.pmt.SetSpot(15);
    d.encoded = pmt.GetEncoded();
    d.side = pmt.GetEastWest();
    d.level = pmt.GetRackLevel();
    d.bay = pmt.GetRackBay();
    d.muxbox = pmt.GetInRack();
    d.tube   = pmt.GetTube();

    // Compute raw mean:
    d.rawGain = 0;
    for(int i=0;i<d.nspot;i++) d.rawGain += d.spots[i].gain/FloatErr(d.nspot);
    // Find the biggest gain and throw it out.
    int   big = 0;
    for(int i=0;i<d.nspot;i++) if(d.spots[i].gain > d.spots[big].gain) { big=i; };
    d.spots[big].bad=3; // make it bad.

    // sort the spots by goodness:
    std::sort(d.spots.begin(),d.spots.begin() + d.nspot);
    
    // Start throwing them out.
    d.ngood = 0;
    d.nbadLow = d.nbadHigh = 0;

    for(Int_t i=0;i<d.nspot;i++) {
      if(d.spots[i].bad ==0) d.ngood++;
      else if(d.spots[i].bad ==1) d.nbadLow++;
      else if(d.spots[i].bad ==2) d.nbadHigh++;
    }
    
    // Because we sorted, the first d.ngood spots are the good ones.
    
    if(d.ngood>0) {

      // Work out the mean gain, using only good spots.
      d.meanGain = 0;
      for(int i=0;i<d.ngood;i++) d.meanGain+=d.spots[i].gain / FloatErr(d.ngood,0);
      
      // The final gain is created using weights.
      d.finalGain = GetWeightedMean(d.spots.begin(), d.spots.begin() + d.ngood);

      // Work out the systematic error.
      int worst = d.ngood-1; // Worst is the last in list.
      if( (d.ngood<d.nspot)               // There are bad ones in the list.
          //&& (d.spots[d.ngood-1].bad<2) // The first bad one isn't VERY bad.
          ) {
          worst = d.ngood;  // Use the first bad one for the error estimate.
      }

      // Without the worst one:
      FloatErr meanWithout(0,0);
      FloatErr meanWith(0,0);

      //meanWithout = GetWeightedMean(d.spots.begin(), d.spots.begin() + worst);
      for(int i=0;i<worst;i++) meanWithout += d.spots[i].gain / FloatErr(worst);

      // With the worst one:
      //meanWith = GetWeightedMean(d.spots.begin(), d.spots.begin() + worst + 1);
      meanWith = (meanWithout*FloatErr(worst) + d.spots[worst].gain) / FloatErr(worst+1);
 
      //cout << meanWithout << "\t" << meanWith << "\t" << meanWithout-meanWith << endl;

      d.sysError = fabs(meanWith - meanWithout);

      MSG("PmtDrift",Msg::kDebug) << "Pmt " << d.pmt.AsString() 
                                  << Form("%3d %3d", d.nspot, d.ngood)
                                  << " \t" << d.finalGain 
                                  << " \t" << d.rawGain 
                                  << " \t" << d.sysError
                                  << " \t" << d.lowestGain 
                                  << " \t" << d.highestGain 
                                  << endl;

    }
  }
}

void PhotonStatSummarizer::ComputeDetectorDrift()
{
  fDetectorDrift.Reset();

  double tot = 0;
  double tot2 = 0;
  double n = 0;

  // Loop spots.
  SpotList_t::iterator it;
  SpotList_t::iterator end= fSpots.end();
  for(it = fSpots.begin(); it!=end; it++) {
    double g = it->second.gain;
    tot += g;
    tot2 += g*g;
    n+=1.0;
  }
  if(n>2) {
    fDetectorDrift.meanSpotGain = tot/n;
    fDetectorDrift. rmsSpotGain = sqrt(fabs(tot2 - tot*tot/n)/(n-1));
  }
  fDetectorDrift.numSpots = n;


  // Loop PMTs.
  tot=tot2=n=0;
  double stots[8];
  double sn[8];
  for(UInt_t i=0;i<8;i++) { stots[i] = sn[i] = 0; };

  PmtList_t::iterator pit;
  PmtList_t::iterator pend = fPmts.end();
  for(pit = fPmts.begin(); pit!=pend; pit++) {
    PlexPixelSpotId psid = pit->second.pmt;
    int section = 0;
    if(psid.GetEastWest()!='E') section+=1;
    if(psid.GetRackLevel()!='L') section+=2;
    if(psid.GetRackBay()>8) section+=4;

    double g = pit->second.finalGain;
    tot += g;
    tot2+= g*g;
    n+=1;
    if(g<fDetectorDrift.lowestPmt) fDetectorDrift.lowestPmt = g;
    if(g>fDetectorDrift.highestPmt) fDetectorDrift.highestPmt = g;    
    stots[section] +=g;
    sn[section] += 1.0;
  }

  if(n>2) {
    fDetectorDrift.meanPmtGain = tot/n;
    fDetectorDrift. rmsPmtGain = sqrt(fabs(tot2 - tot*tot/n)/(n-1));
  }
  fDetectorDrift.numPmts = n;
  for(UInt_t i=0;i<8;i++) {
    if(sn[i]>0)  fDetectorDrift.section[i] = stots[i]/sn[i];
  };
  
}


void PhotonStatSummarizer::FillPmtDatabase(Int_t version,
                                           VldContext vldPointStart,
                                           Int_t task)
{
  if(fPmts.size() == 0) {
    MSG("PmtDrift",Msg::kWarning) << "No PMTs to fill into DB. " << vldPointStart.AsString() << endl;
  }
  VldTimeStamp creationTime = vldPointStart.GetTimeStamp() + VldTimeStamp(version,0);
  VldTimeStamp beginTime = vldPointStart.GetTimeStamp();
  VldTimeStamp endTime = beginTime + VldTimeStamp(3600*24,0); // 1 day later
  VldRange range(vldPointStart.GetDetector(),
                 vldPointStart.GetSimFlag(),
                 beginTime,
                 endTime,
                 "PhotonStatSummarizer::Summarize" );

  // Each crate in the detector is an aggregate, so loop through them all.
  int ncrates = 16;
  if(vldPointStart.GetDetector()==Detector::kNear) ncrates = 8;

  int crates_with_data = 0;
  int rows_total = 0;
  
  for(int crate = 0; crate< ncrates; crate++) {
    // Set up the DB writer.
    DbiWriter<CalPmtDrift> writer(range,
                                  crate, //aggregate,
                                  task, // task no
                                  creationTime,
                                  0,
                                  Form("Code version %d. Created by %s",
                                       version, 
                                       "$Id: PhotonStatSummarizer.cxx,v 1.25 2006/05/27 10:05:44 rhatcher Exp $")
                                  );

    PmtList_t::iterator pit;
    PmtList_t::iterator pend = fPmts.end();
    Int_t nrows = 0;
    
    for(pit = fPmts.begin(); pit!=pend; pit++) {
      PmtDrift& d         = pit->second;
      
      if(d.rcid.GetCrate()==crate) {
        CalPmtDrift cpd(d.pmt,
                        d.rcid,
                        d.finalGain.GetValue(),
                        d.finalGain.GetError(),
                        d.sysError,
                        d.ngood,
                        d.nspot,
                        d.rawGain,
                        d.lowestGain,
                        d.highestGain);
        
        writer << cpd;    
        nrows++;
      }
    }
    
    // Write only if data exists. Otherwise, trash it.
    if(nrows>0){
      writer.Close();
      crates_with_data++;
      rows_total += nrows;
    }
  }

  MSG("PmtDrift",Msg::kInfo) << "Wrote " << rows_total << " rows to DB from " 
                             << fPmts.size() << " PMTs, in " << crates_with_data << " aggregates" << endl;
}


bool PhotonStatSummarizer::DoesTableExist( Int_t version,
                                           VldContext vldPointStart,
                                           Int_t task )
{
  // Question: Does a table already exist with this data?
  // Return true if it does.

  // Compose the VldTimeStamp for the creation date.
  // This is the time the data was taken, plus an integer number of seconds
  // representing the version number of the constants generator.
  VldTimeStamp creationTime = vldPointStart.GetTimeStamp() + VldTimeStamp(version,0);

  const char* sqltxt = Form("(CREATIONDATE='%s') and (TASK=%d) and (DETECTORMASK & %d) and (SIMMASK & %d) order by TIMESTART asc limit 1",
                            creationTime.AsString("s"),
                            task,
                            vldPointStart.GetDetector(),
                            vldPointStart.GetSimFlag() );
  
  DbiSqlContext myContext(sqltxt);
  
  DbiResultPtr<CalPmtDrift> resPtr("CALPMTDRIFT",myContext);
  VldTimeStamp ts= resPtr.GetValidityRec()->GetVldRange().GetTimeStart();
  
  if(ts == vldPointStart.GetTimeStamp()) {
    // This table already exists.
    if(resPtr.GetNumRows()==0) {
      MSG("PmtDrift",Msg::kWarning) << "..Found existing table version " << version << " but no data. Treating as empty." << endl;
      return false;
    }

    return true;
  } 
  return false;
}


void PhotonStatSummarizer::FillTrees( VldContext vldPoint )
{
  static DetectorDrift* det_ptr = 0;
  static PmtDrift* pmt_ptr =0;
  static SpotDrift* spot_ptr =0;
  static VldContext* cx_ptr = 0;
  static Float_t* temp_ptr = 0;
  static Float_t temp;

  if(fFile==0) {
    fFile = new TFile(GetConfig().GetCharString("FileName"),"RECREATE"); 
    fDetTree = (TTree*) fFile->Get("det");
    fPmtTree = (TTree*) fFile->Get("pmts");
    fSpotTree = (TTree*) fFile->Get("spots");
  }
  
  cx_ptr = &vldPoint;
  temp_ptr = &temp;
  Calibrator::Instance().Reset(vldPoint);
  temp = Calibrator::Instance().GetTemperature();
  cout << vldPoint.AsString() << " Temp = " << temp << endl;

  if(fWriteDetTree) {
    if(fPmtTree==0) {
      fDetTree = new TTree("det","det");
      fDetTree->Branch("det_br","DetectorDrift",&det_ptr,32000,3);
      fDetTree->Branch("cx_br","VldContext",&cx_ptr,32000,3);
      fDetTree->Branch("temp_br",temp_ptr,"T/F",32000);
    } else {
      fDetTree->SetBranchAddress("det_br",&det_ptr);
      fDetTree->SetBranchAddress("cx_br",&cx_ptr);
      fDetTree->SetBranchAddress("temp_br",temp_ptr);
    }    
    det_ptr = &fDetectorDrift;
    fDetTree->Fill();
    fDetTree->AutoSave();
  }

  
  if(fWritePmtTree) {
    if(fPmtTree==0) {
      fPmtTree = new TTree("pmts","pmts");
      fPmtTree->Branch("pmt_br","PmtDrift",&pmt_ptr,32000,3);
      fPmtTree->Branch("cx_br","VldContext",&cx_ptr,32000,3);
      fPmtTree->Branch("temp_br",temp_ptr,"T/F",32000);
    } else {
      fPmtTree->SetBranchAddress("pmt_br",&pmt_ptr);
      fPmtTree->SetBranchAddress("cx_br",&cx_ptr);
      fPmtTree->SetBranchAddress("temp_br",temp_ptr);
    }
    
    PmtList_t::iterator pit;
    PmtList_t::iterator pend = fPmts.end();
    for(pit = fPmts.begin(); pit!=pend; pit++) {
      pmt_ptr = &(pit->second);
      fPmtTree->Fill();
    }
    fPmtTree->AutoSave();
  }

  if(fWriteSpotTree) {
    if(fSpotTree==0) {
      fSpotTree = new TTree("spots","spots");
      fSpotTree->Branch("spot_br","SpotDrift",&spot_ptr,32000,3);
      fSpotTree->Branch("cx_br","VldContext",&cx_ptr,32000,3);
      fSpotTree->Branch("temp_br",temp_ptr,"T/F",32000);
    } else {
      fSpotTree->SetBranchAddress("spot_br",&spot_ptr);
      fSpotTree->SetBranchAddress("cx_br",&cx_ptr);
      fSpotTree->SetBranchAddress("temp_br",temp_ptr);
    }

    SpotList_t::iterator it;
    SpotList_t::iterator end= fSpots.end();
    
    for(it = fSpots.begin(); it!=end; it++) {
    spot_ptr = &(it->second);
    fSpotTree->Fill();
    }
    fSpotTree->AutoSave();
  }
}

---