#include <EcalLaserAnalyzer.h>

Inheritance diagram for EcalLaserAnalyzer:

Public Types
enum	VarCol { iBlue, iRed, nColor }
Public Member Functions
virtual void	analyze (const edm::Event &e, const edm::EventSetup &c)
virtual void	beginJob ()
	EcalLaserAnalyzer (const edm::ParameterSet &iConfig)
virtual void	endJob ()
void	setGeomEB (int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr)
void	setGeomEE (int etaG, int phiG, int iX, int iY, int iZ, int module, int tower, int ch, int apdRefTT, int channel, int lmr)
	~EcalLaserAnalyzer ()
Private Attributes
double	_alpha
double	_beta
double	_chi2cut
int	_debug
bool	_docorpn
std::string	_ecalPart
int	_fedid
unsigned int	_firstsample
unsigned int	_firstsamplePN
bool	_fitab
unsigned int	_lastsample
unsigned int	_lastsamplePN
unsigned int	_nevtmax
unsigned int	_niter
double	_noise
unsigned int	_nsamples
unsigned int	_nsamplesPN
unsigned int	_presample
double	_presamplecut
unsigned int	_presamplePN
double	_qualpercent
double	_ratiomaxcutlow
double	_ratiomincuthigh
double	_ratiomincutlow
bool	_saveallevents
unsigned int	_timingcuthigh
unsigned int	_timingcutlow
unsigned int	_timingqualhigh
unsigned int	_timingquallow
double	adc [10]
TFile *	ADCFile
std::string	ADCfile
int	adcG [10]
double *	adcNoPed
TTree *	ADCtrees [NCRYSEB]
std::string	alphafile
std::string	alphainitfile
double	APD [6]
double	apdAmpl
double	apdAmplA
double	apdAmplB
TAPD *	APDAnal [NCRYSEB][nColor]
TFile *	APDFile
std::string	APDfile
TAPD *	APDFirstAnal [NCRYSEB][nColor]
double	APDoAPDA [6]
double	APDoAPDB [6]
double	APDoPN [6]
double	APDoPNA [6]
double	APDoPNB [6]
TAPDPulse *	APDPulse
std::map< int, unsigned int >	apdRefMap [2]
double	apdTime
TTree *	APDtrees [NCRYSEB]
int	channelID
int	channelIteratorEE
std::map< unsigned int, unsigned int >	channelMapEE
int	color
int	colorref
std::vector< int >	colors
int	dccID
TMom *	Delta01
TMom *	Delta12
std::string	digiCollection_
std::string	digiPNCollection_
std::string	digiProducer_
bool	doesABTreeExist
int	eta
int	event
std::string	eventHeaderCollection_
std::string	eventHeaderProducer_
int	eventref
int	fedID
unsigned int	firstChanMod [NMODEE]
int	flag
int	flagAB
int	iChannelID [NCRYSEB]
int	idccID [NCRYSEB]
int	iEta [NCRYSEB]
int	iEvent
unsigned int	iModule [NCRYSEB]
int	iPhi [NCRYSEB]
unsigned int	isFirstChanModFilled [NMODEE]
bool	isGainOK
int	iside [NCRYSEB]
int	IsThereDataADC [NCRYSEB][nColor]
bool	isTimingOK
int	iTowerID [NCRYSEB]
int	iZ
int	laserEvents
int	lightside
TMem *	Mem
int	moduleID
std::vector< int >	modules
unsigned int	nCrys
unsigned int	nevtAB [NCRYSEB]
int	nEvtBadGain [NCRYSEB]
int	nEvtBadTiming [NCRYSEB]
int	nEvtTot [NCRYSEB]
unsigned int	nMod
unsigned int	nPNPerMod
unsigned int	nRefChan
unsigned int	nRefTrees
unsigned int	nSides
int	phi
double	pn [50]
double	PN [6]
double	pn0
double	pn1
double	pnAmpl
TPN *	PNAnal [NMODEE][NPNPERMOD][nColor]
std::string	pncorfile_
TPNCor *	pnCorrector
TPN *	PNFirstAnal [NMODEE][NPNPERMOD][nColor]
int	pnG [50]
int	pnID
double *	pnNoPed
double	PNoPN [6]
double	PNoPNA [6]
double	PNoPNB [6]
TPNPulse *	PNPulse
TTree *	RefAPDtrees [NREFCHAN][NMODEE]
std::string	resdir_
TFile *	resFile
std::string	resfile
TTree *	respntrees [nColor]
TTree *	restrees [nColor]
int	runNum
int	runType
TShapeAnalysis *	shapana
int	side
double	Time [6]
int	towerID
bool	wasABCalcOK [NCRYSEB]
bool	wasGainOK [NCRYSEB]
bool	wasTimingOK [NCRYSEB]

Detailed Description

Definition at line 47 of file EcalLaserAnalyzer.h.

Member Enumeration Documentation

enum EcalLaserAnalyzer::VarCol

Enumerator:

iBlue
iRed
nColor

Definition at line 64 of file EcalLaserAnalyzer.h.

{ iBlue, iRed, nColor };

Constructor & Destructor Documentation

EcalLaserAnalyzer::EcalLaserAnalyzer ( const edm::ParameterSet & iConfig ) [explicit]

Definition at line 58 of file EcalLaserAnalyzer.cc.

References _ecalPart, _fedid, _firstsample, _lastsample, _nsamples, _nsamplesPN, _presample, _presamplePN, _ratiomaxcutlow, _ratiomincuthigh, _ratiomincutlow, _timingcuthigh, _timingcutlow, _timingqualhigh, _timingquallow, APDPulse, Delta01, Delta12, digiCollection_, digiPNCollection_, digiProducer_, eventHeaderCollection_, eventHeaderProducer_, firstChanMod, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), iChannelID, idccID, iEta, iModule, iPhi, isFirstChanModFilled, isGainOK, iside, isTimingOK, iTowerID, iZ, j, ME::lmmodFromDcc(), Mem, modules, nCrys, NCRYSEB, NCRYSEE, nMod, nRefChan, NREFCHAN, pncorfile_, pnCorrector, PNPulse, resdir_, wasGainOK, and wasTimingOK.

  :
iEvent(0),

// Framework parameters with default values

_nsamples(        iConfig.getUntrackedParameter< unsigned int >( "nSamples",           10 ) ),
_presample(       iConfig.getUntrackedParameter< unsigned int >( "nPresamples",         2 ) ),
_firstsample(     iConfig.getUntrackedParameter< unsigned int >( "firstSample",         1 ) ),
_lastsample(      iConfig.getUntrackedParameter< unsigned int >( "lastSample",          2 ) ),
_nsamplesPN(      iConfig.getUntrackedParameter< unsigned int >( "nSamplesPN",         50 ) ), 
_presamplePN(     iConfig.getUntrackedParameter< unsigned int >( "nPresamplesPN",       6 ) ),
_firstsamplePN(   iConfig.getUntrackedParameter< unsigned int >( "firstSamplePN",       7 ) ),
_lastsamplePN(    iConfig.getUntrackedParameter< unsigned int >( "lastSamplePN",        8 ) ),
_timingcutlow(    iConfig.getUntrackedParameter< unsigned int >( "timingCutLow",        2 ) ),
_timingcuthigh(   iConfig.getUntrackedParameter< unsigned int >( "timingCutHigh",       9 ) ),
_timingquallow(   iConfig.getUntrackedParameter< unsigned int >( "timingQualLow",       3 ) ),
_timingqualhigh(  iConfig.getUntrackedParameter< unsigned int >( "timingQualHigh",      8 ) ),
_ratiomincutlow(  iConfig.getUntrackedParameter< double       >( "ratioMinCutLow",    0.4 ) ),
_ratiomincuthigh( iConfig.getUntrackedParameter< double       >( "ratioMinCutHigh",  0.95 ) ),
_ratiomaxcutlow(  iConfig.getUntrackedParameter< double       >( "ratioMaxCutLow",    0.8 ) ),
_presamplecut(    iConfig.getUntrackedParameter< double       >( "presampleCut",      5.0 ) ),
_niter(           iConfig.getUntrackedParameter< unsigned int >( "nIter",               3 ) ),
_fitab(           iConfig.getUntrackedParameter< bool         >( "fitAB",           false ) ),
_alpha(           iConfig.getUntrackedParameter< double       >( "alpha",       1.5076494 ) ),
_beta(            iConfig.getUntrackedParameter< double       >( "beta",        1.5136036 ) ),
_nevtmax(         iConfig.getUntrackedParameter< unsigned int >( "nEvtMax",           200 ) ),
_noise(           iConfig.getUntrackedParameter< double       >( "noise",             2.0 ) ),
_chi2cut(         iConfig.getUntrackedParameter< double       >( "chi2cut",          10.0 ) ), 
_ecalPart(        iConfig.getUntrackedParameter< std::string  >( "ecalPart",         "EB" ) ),
_docorpn(         iConfig.getUntrackedParameter< bool         >( "doCorPN",         false ) ),
_fedid(           iConfig.getUntrackedParameter< int          >( "fedID",            -999 ) ),
_saveallevents(   iConfig.getUntrackedParameter< bool         >( "saveAllEvents",   false ) ),
_qualpercent(     iConfig.getUntrackedParameter< double       >( "qualPercent",       0.2 ) ),
_debug(           iConfig.getUntrackedParameter< int          >( "debug",              0  ) ),
nCrys(                                                                               NCRYSEB),
nPNPerMod(                                                                         NPNPERMOD),
nMod(                                                                                 NMODEE),
nSides(                                                                               NSIDES),
runType(-1), runNum(0), fedID(-1), dccID(-1), side(2), lightside(2), iZ(1),
phi(-1), eta(-1), event(0), color(-1),pn0(0), pn1(0), apdAmpl(0), apdAmplA(0), 
apdAmplB(0),apdTime(0),pnAmpl(0),
pnID(-1), moduleID(-1), channelIteratorEE(0)


  //========================================================================


{

  // Initialization from cfg file

  resdir_                 = iConfig.getUntrackedParameter<std::string>("resDir");
  pncorfile_              = iConfig.getUntrackedParameter<std::string>("pnCorFile");

  digiCollection_         = iConfig.getParameter<std::string>("digiCollection");
  digiPNCollection_       = iConfig.getParameter<std::string>("digiPNCollection");
  digiProducer_           = iConfig.getParameter<std::string>("digiProducer");

  eventHeaderCollection_  = iConfig.getParameter<std::string>("eventHeaderCollection");
  eventHeaderProducer_    = iConfig.getParameter<std::string>("eventHeaderProducer");


  // Geometrical constants initialization 

  if (_ecalPart == "EB") {
    nCrys    = NCRYSEB;
  } else {
    nCrys    = NCRYSEE;
  }
  iZ         =  1;
  if(_fedid <= 609 ) 
    iZ       = -1;
  modules    = ME::lmmodFromDcc(_fedid);
  nMod       = modules.size(); 
  nRefChan   = NREFCHAN;
  
  for(unsigned int j=0;j<nCrys;j++){
    iEta[j]=-1;
    iPhi[j]=-1;
    iModule[j]=10;
    iTowerID[j]=-1;
    iChannelID[j]=-1;
    idccID[j]=-1;
    iside[j]=-1;
    wasTimingOK[j]=true;
    wasGainOK[j]=true;
  }
  
  for(unsigned int j=0;j<nMod;j++){
    int ii= modules[j];
    firstChanMod[ii-1]=0;
    isFirstChanModFilled[ii-1]=0;
  }
  

  // Quality check flags
  
  isGainOK=true;
  isTimingOK=true;

  // PN linearity corrector

  pnCorrector = new TPNCor(pncorfile_.c_str());


  // Objects dealing with pulses

  APDPulse = new TAPDPulse(_nsamples, _presample, _firstsample, _lastsample, 
                           _timingcutlow, _timingcuthigh, _timingquallow, _timingqualhigh,
                           _ratiomincutlow,_ratiomincuthigh, _ratiomaxcutlow);
  PNPulse = new TPNPulse(_nsamplesPN, _presamplePN);

  // Object dealing with MEM numbering

  Mem = new TMem(_fedid);

  // Objects needed for npresample calculation

  Delta01=new TMom();
  Delta12=new TMom();

}

EcalLaserAnalyzer::~EcalLaserAnalyzer ( )

Definition at line 184 of file EcalLaserAnalyzer.cc.

                                     {
  //========================================================================


  // do anything here that needs to be done at destruction time
  // (e.g. close files, deallocate resources etc.)

}

Member Function Documentation

void EcalLaserAnalyzer::analyze	(	const edm::Event &	e,
		const edm::EventSetup &	c
	)		`[virtual]`

Implements edm::EDAnalyzer.

Definition at line 310 of file EcalLaserAnalyzer.cc.

References _debug, _docorpn, _ecalPart, _fedid, _firstsamplePN, _fitab, _lastsamplePN, _nevtmax, _nsamplesPN, adc, adcG, ADCtrees, TMom::addEntry(), APDPulse, MEEBGeom::apdRefTower(), MEEEGeom::apdRefTower(), edm::DataFrameContainer::begin(), edm::SortedCollection< T, SORT >::begin(), ExpressReco_HICollisions_FallBack::cerr, channelID, EcalElectronicsId::channelId(), channelIteratorEE, channelMapEE, color, colors, corr, gather_cfg::cout, dccID, MEEEGeom::dee(), Delta01, Delta12, digiCollection_, digiPNCollection_, digiProducer_, doesABTreeExist, MEEBGeom::electronic_channel(), edm::DataFrameContainer::end(), edm::SortedCollection< T, SORT >::end(), eta, eventHeaderCollection_, eventHeaderProducer_, exception, fedID, spr::find(), edm::EventSetup::get(), TPNPulse::getAdcWithoutPedestal(), edm::Event::getByLabel(), TAPDPulse::getDelta(), EcalElectronicsMapping::getElectronicsId(), TPNPulse::getMaxSample(), TPNCor::getPNCorrectionFactor(), ecalpyutils::hashedIndex(), i, EcalPnDiodeDetId::iDCCId(), iEvent, init, EcalPnDiodeDetId::iPnId(), TMem::isMemRelevant(), TAPDPulse::isPulseOK(), TAPDPulse::isTimingQualOK(), iZ, EcalDCCHeaderBlock::LASER_DELAY_SCAN, EcalDCCHeaderBlock::LASER_GAP, EcalDCCHeaderBlock::LASER_POWER_SCAN, EcalDCCHeaderBlock::LASER_STD, laserEvents, lightside, MEEEGeom::lmmod(), MEEBGeom::lmmod(), MEEEGeom::lmr(), MEEBGeom::lmr(), MEEBGeom::localCoord(), siStripFEDMonitor_P5_cff::Max, TMem::Mem(), Mem, modules, nCrys, nevtAB, nEvtBadGain, nEvtBadTiming, nEvtTot, phi, pn, pn0, pn1, pnAmpl, pnCorrector, pnG, pnNoPed, PNPulse, edm::ESHandle< T >::product(), edm::Handle< T >::product(), runNum, runType, setGeomEB(), setGeomEE(), TAPDPulse::setPulse(), TPNPulse::setPulse(), shapana, side, findQualityFiles::size, strip(), EcalElectronicsId::stripId(), towerID, EcalElectronicsId::towerId(), EcalDCCHeaderBlock::EcalDCCEventSettings::wavelength, and EcalElectronicsId::xtalId().

                                                                            {
//========================================================================

  ++iEvent;
  
  // retrieving DCC header

  edm::Handle<EcalRawDataCollection> pDCCHeader;
  const  EcalRawDataCollection* DCCHeader=0;
  try {
    e.getByLabel(eventHeaderProducer_,eventHeaderCollection_, pDCCHeader);
    DCCHeader=pDCCHeader.product();
  }catch ( std::exception& ex ) {
    std::cerr << "Error! can't get the product retrieving DCC header" << 
      eventHeaderCollection_.c_str() <<" "<< eventHeaderProducer_.c_str() << std::endl;
  }
 
  //retrieving crystal data from Event

  edm::Handle<EBDigiCollection>  pEBDigi;
  const  EBDigiCollection* EBDigi=0;
  edm::Handle<EEDigiCollection>  pEEDigi;
  const  EEDigiCollection* EEDigi=0;


  if (_ecalPart == "EB") {
    try {
      e.getByLabel(digiProducer_,digiCollection_, pEBDigi); 
      EBDigi=pEBDigi.product(); 
    }catch ( std::exception& ex ) {
      std::cerr << "Error! can't get the product retrieving EB crystal data " << 
        digiCollection_.c_str() << std::endl;
    } 
  } else if (_ecalPart == "EE") {
    try {
      e.getByLabel(digiProducer_,digiCollection_, pEEDigi); 
      EEDigi=pEEDigi.product(); 
    }catch ( std::exception& ex ) {
      std::cerr << "Error! can't get the product retrieving EE crystal data " << 
        digiCollection_.c_str() << std::endl;
    } 
  } else {
    std::cout <<" Wrong ecalPart in cfg file " << std::endl;
    return;
  }
  
  // retrieving crystal PN diodes from Event
  
  edm::Handle<EcalPnDiodeDigiCollection>  pPNDigi;
  const  EcalPnDiodeDigiCollection* PNDigi=0;
  try {
    e.getByLabel(digiProducer_, pPNDigi);
    PNDigi=pPNDigi.product(); 
  }catch ( std::exception& ex ) {
    std::cerr << "Error! can't get the product " << digiPNCollection_.c_str() << std::endl;
  }

  // retrieving electronics mapping

  edm::ESHandle< EcalElectronicsMapping > ecalmapping;
  const EcalElectronicsMapping* TheMapping=0; 
  try{
    c.get< EcalMappingRcd >().get(ecalmapping);
    TheMapping = ecalmapping.product();
  }catch ( std::exception& ex ) {
    std::cerr << "Error! can't get the product EcalMappingRcd"<< std::endl;
  }


  // ============================
  // Decode DCCHeader Information 
  // ============================
  
  for ( EcalRawDataCollection::const_iterator headerItr= DCCHeader->begin();
        headerItr != DCCHeader->end(); ++headerItr ) {
    
    // Get run type and run number 

    int fed = headerItr->fedId();  
    if(fed!=_fedid && _fedid!=-999) continue; 
    
    runType=headerItr->getRunType();
    runNum=headerItr->getRunNumber();
    event=headerItr->getLV1();

    dccID=headerItr->getDccInTCCCommand();
    fedID=headerItr->fedId();  
    lightside=headerItr->getRtHalf();

    // Check fed corresponds to the DCC in TCC
    
    if( 600+dccID != fedID ) continue;

    // Cut on runType

    if( runType!=EcalDCCHeaderBlock::LASER_STD 
       && runType!=EcalDCCHeaderBlock::LASER_GAP 
       && runType!=EcalDCCHeaderBlock::LASER_POWER_SCAN
       && runType!=EcalDCCHeaderBlock::LASER_DELAY_SCAN ) return; 
    
    // Retrieve laser color and event number
    
    EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
    color = settings.wavelength;
    if( color<0 ) return;
    
    std::vector<int>::iterator iter = find( colors.begin(), colors.end(), color );
    if( iter==colors.end() ){
      colors.push_back( color );
      std::cout <<" new color found "<< color<<" "<< colors.size()<< std::endl;
    }
  }

  // Cut on fedID
 
  if(fedID!=_fedid && _fedid!=-999) return; 
    
  // Count laser events

  laserEvents++;


  // ======================
  // Decode PN Information
  // ======================
  
  TPNFit * pnfit = new TPNFit();
  pnfit -> init(_nsamplesPN,_firstsamplePN,  _lastsamplePN);
  
  double chi2pn=0;
  unsigned int samplemax=0;
  int  pnGain=0;
  
  std::map <int, std::vector<double> > allPNAmpl;
  std::map <int, std::vector<double> > allPNGain;

  
  // Loop on PNs digis
  
  for ( EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin();
        pnItr != PNDigi->end(); ++pnItr ) { 
    
    EcalPnDiodeDetId pnDetId = EcalPnDiodeDetId((*pnItr).id());
    
    if (_debug==1) std::cout <<"-- debug -- Inside PNDigi - pnID=" <<
                     pnDetId.iPnId()<<", dccID="<< pnDetId.iDCCId()<< std::endl;

    // Skip MEM DCC without relevant data
  
    bool isMemRelevant=Mem->isMemRelevant(pnDetId.iDCCId());
    if(!isMemRelevant) continue;

    // Loop on PN samples

    for ( int samId=0; samId < (*pnItr).size() ; samId++ ) {   
      pn[samId]=(*pnItr).sample(samId).adc();  
      pnG[samId]=(*pnItr).sample(samId).gainId(); 
      if (samId==0) pnGain=pnG[samId];
      if (samId>0) pnGain=int(TMath::Max(pnG[samId],pnGain));     
    }
    
    if( pnGain!=1 ) std::cout << "PN gain different from 1"<< std::endl;
    
    // Calculate amplitude from pulse
    
    PNPulse->setPulse(pn);
    pnNoPed=PNPulse->getAdcWithoutPedestal();
    samplemax=PNPulse->getMaxSample();
    chi2pn = pnfit -> doFit(samplemax,&pnNoPed[0]); 
    if(chi2pn == 101 || chi2pn == 102 || chi2pn == 103) pnAmpl=0.;
    else pnAmpl= pnfit -> getAmpl();
    
    // Apply linearity correction

    double corr=1.0;
    if( _docorpn ) corr=pnCorrector->getPNCorrectionFactor(pnAmpl, pnGain);
    pnAmpl*=corr;

    // Fill PN ampl vector

    allPNAmpl[pnDetId.iDCCId()].push_back(pnAmpl);
      
    if (_debug==1) std::cout <<"-- debug -- Inside PNDigi - PNampl=" << 
                     pnAmpl<<", PNgain="<< pnGain<<std::endl;  
  }
  
  // ===========================
  // Decode EBDigis Information
  // ===========================
  
  int  adcGain=0;

  if (EBDigi){
    
    // Loop on crystals
    //=================== 
    
    for ( EBDigiCollection::const_iterator digiItr= EBDigi->begin();
          digiItr != EBDigi->end(); ++digiItr ) { 
      
      // Retrieve geometry
      //=================== 
      
      EBDetId id_crystal(digiItr->id()) ;
      EBDataFrame df( *digiItr );
      EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
      
      int etaG = id_crystal.ieta() ;  // global
      int phiG = id_crystal.iphi() ;  // global
      
      std::pair<int, int> LocalCoord=MEEBGeom::localCoord( etaG , phiG );
      
      int etaL=LocalCoord.first ; // local
      int phiL=LocalCoord.second ;// local
      
      int strip=elecid_crystal.stripId();
      int xtal=elecid_crystal.xtalId();
      
      int module= MEEBGeom::lmmod(etaG, phiG);
      int tower=elecid_crystal.towerId();
           
      int apdRefTT=MEEBGeom::apdRefTower(module);      
           
      std::pair<int,int> pnpair=MEEBGeom::pn(module);
      unsigned int MyPn0=pnpair.first;
      unsigned int MyPn1=pnpair.second;
      
      int lmr=MEEBGeom::lmr( etaG,phiG ); 
      unsigned int channel=MEEBGeom::electronic_channel( etaL, phiL );
      assert( channel < nCrys );

      setGeomEB(etaG, phiG, module, tower, strip, xtal, apdRefTT, channel, lmr);      
      
      if (_debug==1) std::cout << "-- debug -- Inside EBDigi - towerID:"<< towerID<<
                       " channelID:" <<channelID<<" module:"<< module<<
                       " modules:"<<modules.size()<< std::endl;
      
      // APD Pulse
      //=========== 

      // Loop on adc samples 
     
      for (unsigned int i=0; i< (*digiItr).size() ; ++i ) {   

        EcalMGPASample samp_crystal(df.sample(i));
        adc[i]=samp_crystal.adc() ;    
        adcG[i]=samp_crystal.gainId();   
        adc[i]*=adcG[i];
        if (i==0) adcGain=adcG[i];
        if (i>0) adcGain=TMath::Max(adcG[i],adcGain);  
      }

      APDPulse->setPulse(adc);
      
      // Quality checks
      //================
      
      if(adcGain!=1) nEvtBadGain[channel]++;   
      if(!APDPulse->isTimingQualOK()) nEvtBadTiming[channel]++;
      nEvtTot[channel]++;

      // Associate PN ampl
      //===================
 
      int mem0=Mem->Mem(lmr,0);
      int mem1=Mem->Mem(lmr,1);

      if(allPNAmpl[mem0].size()>MyPn0) pn0=allPNAmpl[mem0][MyPn0];
      else pn0=0;
      if(allPNAmpl[mem1].size()>MyPn1) pn1=allPNAmpl[mem1][MyPn1];
      else pn1=0;


      // Fill if Pulse is fine
      //=======================
      
      if( APDPulse->isPulseOK() && lightside==side){
        
        ADCtrees[channel]->Fill(); 
        
        Delta01->addEntry(APDPulse->getDelta(0,1));
        Delta12->addEntry(APDPulse->getDelta(1,2));
        
        if( nevtAB[channel] < _nevtmax && _fitab ){
          if(doesABTreeExist)  shapana -> putAllVals(channel, adc, eta, phi);   
          else shapana -> putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);
          nevtAB[channel]++ ;
        }
      }
    }
    
  } else if (EEDigi) {

    // Loop on crystals
    //===================
    
    for ( EEDigiCollection::const_iterator digiItr= EEDigi->begin(); 
          digiItr != EEDigi->end(); ++digiItr ) {  
      
      // Retrieve geometry
      //=================== 
      
      EEDetId id_crystal(digiItr->id()) ;
      EEDataFrame df( *digiItr );
      EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
       
      int etaG = id_crystal.iy() ; 
      int phiG = id_crystal.ix() ;

      int iX = (phiG-1)/5+1;
      int iY = (etaG-1)/5+1;

      int tower=elecid_crystal.towerId();
      int ch=elecid_crystal.channelId()-1; 

      int module=MEEEGeom::lmmod( iX, iY );
      if( module>=18 && side==1 ) module+=2;
      int lmr=MEEEGeom::lmr( iX, iY ,iZ); 
      int dee=MEEEGeom::dee(lmr);
      int apdRefTT=MEEEGeom::apdRefTower(lmr, module);  
      
      std::pair<int,int> pnpair=MEEEGeom::pn( dee, module ) ; 
      unsigned int MyPn0=pnpair.first;
      unsigned int MyPn1=pnpair.second;

      int hashedIndex=100000*eta+phi;
      if( channelMapEE.count(hashedIndex) == 0 ){
        channelMapEE[hashedIndex]=channelIteratorEE;
        channelIteratorEE++;
      }
      unsigned int channel=channelMapEE[hashedIndex];
      assert ( channel < nCrys );

      setGeomEE(etaG, phiG, iX, iY, iZ, module, tower, ch, apdRefTT, channel, lmr);
      
      
      if (_debug==1) std::cout << "-- debug -- Inside EEDigi - towerID:"<< towerID<<
                       " channelID:" <<channelID<<" module:"<< module<<
                       " modules:"<<modules.size()<< std::endl;
      
      // APD Pulse
      //=========== 

      if( (*digiItr).size()>10) std::cout <<"SAMPLES SIZE > 10!" <<  (*digiItr).size()<< std::endl;
 
      // Loop on adc samples  

      for (unsigned int i=0; i< (*digiItr).size() ; ++i ) { 

        EcalMGPASample samp_crystal(df.sample(i));
        adc[i]=samp_crystal.adc() ;    
        adcG[i]=samp_crystal.gainId();   
        adc[i]*=adcG[i];
        
        if (i==0) adcGain=adcG[i];
        if (i>0) adcGain=TMath::Max(adcG[i],adcGain);  
      }
      
      APDPulse->setPulse(adc);
      
      // Quality checks
      //================
      
      if(adcGain!=1) nEvtBadGain[channel]++;   
      if(!APDPulse->isTimingQualOK()) nEvtBadTiming[channel]++;
      nEvtTot[channel]++;

      
      // Associate PN ampl
      //===================

      int mem0=Mem->Mem(lmr,0);
      int mem1=Mem->Mem(lmr,1);
      
      if(allPNAmpl[mem0].size()>MyPn0) pn0=allPNAmpl[mem0][MyPn0];
      else pn0=0;
      if(allPNAmpl[mem1].size()>MyPn1) pn1=allPNAmpl[mem1][MyPn1];
      else pn1=0;

      // Fill if Pulse is fine
      //=======================
      
      if( APDPulse->isPulseOK() && lightside==side){
        ADCtrees[channel]->Fill(); 
        
        Delta01->addEntry(APDPulse->getDelta(0,1));
        Delta12->addEntry(APDPulse->getDelta(1,2));
        
        if( nevtAB[channel] < _nevtmax && _fitab ){
          if(doesABTreeExist)  shapana -> putAllVals(channel, adc, eta, phi);   
          else shapana -> putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);  
          nevtAB[channel]++ ;
        }
      }
    }
  }
}

void EcalLaserAnalyzer::beginJob ( void ) [virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 196 of file EcalLaserAnalyzer.cc.

References _alpha, _beta, _chi2cut, _firstsample, _fitab, _lastsample, _nevtmax, _noise, _nsamples, _presample, adc, ADCfile, ADCFile, ADCtrees, alphafile, alphainitfile, APDfile, channelID, color, dccID, doesABTreeExist, eta, event, i, laserEvents, AlCaRecoCosmics_cfg::name, nCrys, nevtAB, NULL, phi, pn0, pn1, resdir_, resfile, shapana, side, and towerID.

                                 {
  //========================================================================

  
  // Create temporary files and trees to save adc samples
  //======================================================

  ADCfile=resdir_;
  ADCfile+="/APDSamplesLaser.root";

  APDfile=resdir_;
  APDfile+="/APDPNLaserAllEvents.root";

  ADCFile = new TFile(ADCfile.c_str(),"RECREATE");
  
  for (unsigned int i=0;i<nCrys;i++){
    
    std::stringstream name;
    name << "ADCTree" <<i+1;
    ADCtrees[i]= new TTree(name.str().c_str(),name.str().c_str());

    ADCtrees[i]->Branch( "ieta",        &eta,         "eta/I"         );
    ADCtrees[i]->Branch( "iphi",        &phi,         "phi/I"         );
    ADCtrees[i]->Branch( "side",        &side,        "side/I"        );
    ADCtrees[i]->Branch( "dccID",       &dccID,       "dccID/I"       );
    ADCtrees[i]->Branch( "towerID",     &towerID,     "towerID/I"     );
    ADCtrees[i]->Branch( "channelID",   &channelID,   "channelID/I"   );
    ADCtrees[i]->Branch( "event",       &event,       "event/I"       );
    ADCtrees[i]->Branch( "color",       &color,       "color/I"       );
    ADCtrees[i]->Branch( "adc",         &adc ,        "adc[10]/D"     );
    ADCtrees[i]->Branch( "pn0",         &pn0 ,        "pn0/D"         );
    ADCtrees[i]->Branch( "pn1",         &pn1 ,        "pn1/D"         );

     
    ADCtrees[i]->SetBranchAddress( "ieta",        &eta         );  
    ADCtrees[i]->SetBranchAddress( "iphi",        &phi         ); 
    ADCtrees[i]->SetBranchAddress( "side",        &side        ); 
    ADCtrees[i]->SetBranchAddress( "dccID",       &dccID       ); 
    ADCtrees[i]->SetBranchAddress( "towerID",     &towerID     ); 
    ADCtrees[i]->SetBranchAddress( "channelID",   &channelID   ); 
    ADCtrees[i]->SetBranchAddress( "event",       &event       );   
    ADCtrees[i]->SetBranchAddress( "color",       &color       );   
    ADCtrees[i]->SetBranchAddress( "adc",         adc          );
    ADCtrees[i]->SetBranchAddress( "pn0",         &pn0         );
    ADCtrees[i]->SetBranchAddress( "pn1",         &pn1         );    
 
    nevtAB[i]=0 ;
  } 


  // Define output results filenames and shape analyzer object (alpha,beta) 
  //=====================================================================
  
  // 1) AlphaBeta files
  
  doesABTreeExist=true;

  std::stringstream nameabinitfile;
  nameabinitfile << resdir_ <<"/ABInit.root";
  alphainitfile=nameabinitfile.str();
  
  std::stringstream nameabfile;  
  nameabfile << resdir_ <<"/AB.root";
  alphafile=nameabfile.str();
  
  FILE *test;
  if(_fitab)
    test = fopen(alphainitfile.c_str(),"r"); 
  else 
    test = fopen(alphafile.c_str(),"r"); 
  if(test == NULL) {
    doesABTreeExist=false;
    _fitab=true;
  };
  delete test;
  
  TFile *fAB=0; TTree *ABInit=0;
  if(doesABTreeExist){
    fAB=new TFile(nameabinitfile.str().c_str());
  }
  if(doesABTreeExist && fAB){
    ABInit = (TTree*) fAB->Get("ABCol0");
  }
  
  // 2) Shape analyzer
  
  if(doesABTreeExist && fAB && ABInit && ABInit->GetEntries()!=0){
    shapana= new TShapeAnalysis(ABInit, _alpha, _beta, 5.5, 1.0);
    doesABTreeExist=true;
  }else{
    shapana= new TShapeAnalysis(_alpha, _beta, 5.5, 1.0);
    doesABTreeExist=false;
    _fitab=true;  
  }  
  shapana -> set_const(_nsamples,_firstsample,_lastsample,
                       _presample, _nevtmax, _noise, _chi2cut);
  
  if(doesABTreeExist && fAB ) fAB->Close();
  

  //  2) APD file
  
  std::stringstream nameapdfile;
  nameapdfile << resdir_ <<"/APDPN_LASER.root";      
  resfile=nameapdfile.str();

  
  // Laser events counter

  laserEvents=0;

}

void EcalLaserAnalyzer::endJob ( void ) [virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 711 of file EcalLaserAnalyzer.cc.

References _debug, _ecalPart, _firstsample, _fitab, _lastsample, _niter, _nsamples, _presample, _presamplecut, _qualpercent, _saveallevents, adc, ADCfile, ADCFile, adcNoPed, ADCtrees, TPN::addEntry(), TAPD::addEntry(), alpha, alphafile, alphainitfile, APD, apdAmpl, apdAmplA, apdAmplB, APDAnal, APDfile, APDFile, APDFirstAnal, APDoAPDA, APDoAPDB, APDoPN, APDoPNA, APDoPNB, APDPulse, apdRefMap, apdTime, APDtrees, asciidump::at, beta, channelID, channelMapEE, ExpressReco_HICollisions_FallBack::chi2, color, colorref, colors, TShapeAnalysis::computeShape(), gather_cfg::cout, align_tpl::cut, dccID, Delta01, Delta12, doesABTreeExist, event, eventref, firstChanMod, flag, flagAB, TAPDPulse::getAdcWithoutPedestal(), TAPD::getAPD(), TAPD::getAPDoAPD0(), TAPD::getAPDoAPD1(), TAPD::getAPDoPN(), TAPD::getAPDoPN0(), TAPD::getAPDoPN1(), TMom::getMean(), TPN::getPN(), TPN::getPNoPN(), TPN::getPNoPN0(), TPN::getPNoPN1(), TAPD::getTime(), edm::service::getTime(), TShapeAnalysis::getVals(), i, iChannelID, idccID, iEta, iModule, init, iPhi, isGainOK, iside, TAPDPulse::isPulseOK(), IsThereDataADC, isTimingOK, iTowerID, laserEvents, moduleID, modules, AlCaRecoCosmics_cfg::name, nColor, nCrys, nEvtBadGain, nEvtBadTiming, nEvtTot, nMod, nPNPerMod, nRefChan, PN, pn0, pn1, PNAnal, PNFirstAnal, pnID, PNoPN, PNoPNA, PNoPNB, RefAPDtrees, resFile, resfile, respntrees, restrees, TAPD::set2DAPDoAPD0Cut(), TAPD::set2DAPDoAPD1Cut(), TShapeAnalysis::set_nch(), TShapeAnalysis::set_presample(), TAPD::setAPDCut(), TAPD::setAPDoPN0Cut(), TAPD::setAPDoPN1Cut(), TAPD::setAPDoPNCut(), TPN::setPNCut(), TAPDPulse::setPresamples(), TAPDPulse::setPulse(), TAPD::setTimeCut(), shapana, side, Time, towerID, wasGainOK, and wasTimingOK.

                               {
//========================================================================

  // Adjust channel numbers for EE 
  //===============================

  if( _ecalPart == "EE" ) {
    nCrys=channelMapEE.size();
    shapana->set_nch(nCrys);
  }
  
  // Set presamples number 
  //======================
  
  double delta01=Delta01->getMean();
  double delta12=Delta12->getMean();
  if(delta12>_presamplecut) {
    _presample=2;
    if(delta01>_presamplecut) _presample=1;
  }
  
  APDPulse->setPresamples(_presample);
  shapana->set_presample(_presample);

  //  Get alpha and beta
  //======================

  if(_fitab){
    std::cout <<  "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
    std::cout <<   "\t+=+     Analyzing data: getting (alpha, beta)     +=+" << std::endl;
    TFile *fAB=0; TTree *ABInit=0;
    if(doesABTreeExist){
      fAB=new TFile(alphainitfile.c_str());
    }
    if(doesABTreeExist && fAB){
      ABInit = (TTree*) fAB->Get("ABCol0");
    }
    shapana->computeShape(alphafile, ABInit);
    std::cout <<    "\t+=+    .................................... done  +=+" << std::endl;
    std::cout <<    "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
  }
  

  // Don't do anything if there is no events
  //=========================================

  if( laserEvents == 0 ){
    
    ADCFile->Close(); 
    std::stringstream del;
    del << "rm " <<ADCfile;
    system(del.str().c_str());
    std::cout << " No Laser Events "<< std::endl;
    return;
  }

  // Set quality flags for gains and timing
  //=========================================

  double BadGainEvtPercentage=0.0;
  double BadTimingEvtPercentage=0.0;
  
  int nChanBadGain=0;
  int nChanBadTiming=0;
  
  for (unsigned int i=0;i<nCrys;i++){
    if(nEvtTot[i]!=0){
      BadGainEvtPercentage=double(nEvtBadGain[i])/double(nEvtTot[i]);
      BadTimingEvtPercentage=double(nEvtBadTiming[i])/double(nEvtTot[i]);
    }
    if(BadGainEvtPercentage>_qualpercent) {
      wasGainOK[i]=false;
      nChanBadGain++;
    }
    if(BadTimingEvtPercentage>_qualpercent){
      wasTimingOK[i]=false;
      nChanBadTiming++;
    }
  }
  
  double BadGainChanPercentage=double(nChanBadGain)/double(nCrys);
  double BadTimingChanPercentage=double(nChanBadTiming)/double(nCrys);
  
  if(BadGainChanPercentage>_qualpercent) isGainOK = false;
  if(BadTimingChanPercentage>_qualpercent) isTimingOK = false;

  // Analyze adc samples to get amplitudes
  //=======================================
  
  std::cout <<  "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
  std::cout <<    "\t+=+     Analyzing laser data: getting APD, PN, APD/PN, PN/PN    +=+" << std::endl;
  
  if( !isGainOK )
  std::cout <<    "\t+=+ ............................ WARNING! APD GAIN WAS NOT 1    +=+" << std::endl;
  if( !isTimingOK )
  std::cout <<    "\t+=+ ............................ WARNING! TIMING WAS BAD        +=+" << std::endl;
    
  
  APDFile = new TFile(APDfile.c_str(),"RECREATE");
  
  int ieta, iphi, channelID, towerID, flag;

  for (unsigned int i=0;i<nCrys;i++){
    
    std::stringstream name;
    name << "APDTree" <<i+1;
    
    APDtrees[i]= new TTree(name.str().c_str(),name.str().c_str());
    
    //List of branches
    
    APDtrees[i]->Branch( "event",     &event,     "event/I"     );
    APDtrees[i]->Branch( "color",     &color,     "color/I"     );
    APDtrees[i]->Branch( "iphi",      &iphi,      "iphi/I"      );
    APDtrees[i]->Branch( "ieta",      &ieta,      "ieta/I"      );
    APDtrees[i]->Branch( "side",      &side,      "side/I"      );
    APDtrees[i]->Branch( "dccID",     &dccID,     "dccID/I"     );
    APDtrees[i]->Branch( "towerID",   &towerID,   "towerID/I"   );
    APDtrees[i]->Branch( "channelID", &channelID, "channelID/I" );
    APDtrees[i]->Branch( "apdAmpl",   &apdAmpl,   "apdAmpl/D"   );
    APDtrees[i]->Branch( "apdTime",   &apdTime,   "apdTime/D"   );
    if(_saveallevents) APDtrees[i]->Branch( "adc", &adc ,"adc[10]/D" );
    APDtrees[i]->Branch( "flag",      &flag,      "flag/I"      ); 
    APDtrees[i]->Branch( "flagAB",    &flagAB,    "flagAB/I"    ); 
    APDtrees[i]->Branch( "pn0",       &pn0,       "pn0/D"       );
    APDtrees[i]->Branch( "pn1",       &pn1,       "pn1/D"       ); 
    
    APDtrees[i]->SetBranchAddress( "event",     &event     );
    APDtrees[i]->SetBranchAddress( "color",     &color     ); 
    APDtrees[i]->SetBranchAddress( "iphi",      &iphi      );
    APDtrees[i]->SetBranchAddress( "ieta",      &ieta      );
    APDtrees[i]->SetBranchAddress( "side",      &side      );
    APDtrees[i]->SetBranchAddress( "dccID",     &dccID     );
    APDtrees[i]->SetBranchAddress( "towerID",   &towerID   );
    APDtrees[i]->SetBranchAddress( "channelID", &channelID ); 
    APDtrees[i]->SetBranchAddress( "apdAmpl",   &apdAmpl   );
    APDtrees[i]->SetBranchAddress( "apdTime",   &apdTime   );
    if(_saveallevents)APDtrees[i]->SetBranchAddress( "adc", adc );
    APDtrees[i]->SetBranchAddress( "flag",      &flag      ); 
    APDtrees[i]->SetBranchAddress( "flagAB",    &flagAB    ); 
    APDtrees[i]->SetBranchAddress( "pn0",       &pn0       );
    APDtrees[i]->SetBranchAddress( "pn1",       &pn1       ); 
  
  }


  for (unsigned int iref=0;iref<nRefChan;iref++){
    for (unsigned int imod=0;imod<nMod;imod++){
      
      int jmod=modules[imod];
      
      std::stringstream nameref;
      nameref << "refAPDTree" <<imod<<"_"<<iref;
      
      RefAPDtrees[iref][jmod]= new TTree(nameref.str().c_str(),nameref.str().c_str());
      
      RefAPDtrees[iref][jmod]->Branch( "eventref",     &eventref,     "eventref/I"     );
      RefAPDtrees[iref][jmod]->Branch( "colorref",     &colorref,     "colorref/I"     );
      if(iref==0) RefAPDtrees[iref][jmod]->Branch( "apdAmplA",   &apdAmplA,   "apdAmplA/D"   );
      if(iref==1) RefAPDtrees[iref][jmod]->Branch( "apdAmplB",   &apdAmplB,   "apdAmplB/D"   );
      
      RefAPDtrees[iref][jmod]->SetBranchAddress( "eventref",     &eventref     );
      RefAPDtrees[iref][jmod]->SetBranchAddress( "colorref",     &colorref     ); 
      if(iref==0)RefAPDtrees[iref][jmod]->SetBranchAddress( "apdAmplA",   &apdAmplA   );
      if(iref==1)RefAPDtrees[iref][jmod]->SetBranchAddress( "apdAmplB",   &apdAmplB   );
      
    } 
  }
 
  assert( colors.size()<= nColor );
  unsigned int nCol=colors.size();

  // Declare PN stuff
  //===================
  
  for (unsigned int iM=0;iM<nMod;iM++){
    unsigned int iMod=modules[iM]-1;

    for (unsigned int ich=0;ich<nPNPerMod;ich++){
      for (unsigned int icol=0;icol<nCol;icol++){
        PNFirstAnal[iMod][ich][icol]=new TPN(ich);
        PNAnal[iMod][ich][icol]=new TPN(ich);
      }

    }
  }
  

  // Declare function for APD ampl fit
  //===================================

  PulseFitWithFunction * pslsfit = new PulseFitWithFunction();
  double chi2;

  
  for (unsigned int iCry=0;iCry<nCrys;iCry++){
    for (unsigned int icol=0;icol<nCol;icol++){

      // Declare APD stuff
      //===================
      
      APDFirstAnal[iCry][icol] = new TAPD();
      IsThereDataADC[iCry][icol]=1;      
      std::stringstream cut;
      cut <<"color=="<<colors[icol];
      if(ADCtrees[iCry]->GetEntries(cut.str().c_str())<10) IsThereDataADC[iCry][icol]=0;  
      
    }

    unsigned int iMod=iModule[iCry]-1;
    double alpha, beta;
    
    // Loop on events
    //================
    
    Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry=0; jentry< ADCtrees[iCry]->GetEntriesFast();jentry++){
      nb = ADCtrees[iCry]->GetEntry(jentry);   nbytes += nb;    
      
      // Get back color
 
      unsigned int iCol=0;
      for(unsigned int i=0;i<nCol;i++){
        if(color==colors[i]) {
          iCol=i;
          i=colors.size();
        }
      }
      
      // Retreive alpha and beta
      
      std::vector<double> abvals=shapana->getVals(iCry);
      alpha=abvals[0];
      beta=abvals[1];
      flagAB=int(abvals[4]);
      iphi=iPhi[iCry];
      ieta=iEta[iCry];
      towerID=iTowerID[iCry];
      channelID=iChannelID[iCry];
      
      // Amplitude calculation

      APDPulse->setPulse(adc);
      adcNoPed=APDPulse->getAdcWithoutPedestal();

      apdAmpl=0;
      apdAmplA=0;
      apdAmplB=0;
      apdTime=0;
      
      if( APDPulse->isPulseOK()){
        
        pslsfit -> init(_nsamples,_firstsample,_lastsample,_niter,  alpha,  beta);
        chi2 = pslsfit -> doFit(&adcNoPed[0]);
        
        if( chi2 < 0. || chi2 == 102 || chi2 == 101  ) {          
          apdAmpl=0;      
          apdTime=0;
          flag=0;
        }else{
          apdAmpl = pslsfit -> getAmpl();
          apdTime = pslsfit -> getTime();
          flag=1;
        }
      }else {
        chi2=0.0;
        apdAmpl=0;        
        apdTime=0;
        flag=0;
      }
      
      if (_debug==1) std::cout <<"-- debug test -- apdAmpl="<<apdAmpl<< 
                       ", apdTime="<< apdTime<< std::endl;
      double pnmean;
      if (pn0<10 && pn1>10) {
        pnmean=pn1;
      }else if (pn1<10 && pn0>10){
        pnmean=pn0;
      }else pnmean=0.5*(pn0+pn1);
      
      if (_debug==1) std::cout <<"-- debug test -- pn0="<<pn0<<", pn1="<<pn1<< std::endl;
      
      // Fill PN stuff
      //===============

      if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){ 
        for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
          PNFirstAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
        }
      }
      
      // Fill APD stuff 
      //================

      if( APDPulse->isPulseOK() ){
        APDFirstAnal[iCry][iCol]->addEntry(apdAmpl,pnmean, pn0, pn1, apdTime);      
        APDtrees[iCry]->Fill();    
        
        // Fill reference trees
        //=====================
        
        if( apdRefMap[0][iMod+1]==iCry || apdRefMap[1][iMod+1]==iCry ) {
          
          apdAmplA=0.0;
          apdAmplB=0.0;      
          eventref=event;
          colorref=color;
          
          for(unsigned int ir=0;ir<nRefChan;ir++){
            
            if(apdRefMap[ir][iMod+1]==iCry) {
              if (ir==0) apdAmplA=apdAmpl;
              else if(ir==1) apdAmplB=apdAmpl;
              RefAPDtrees[ir][iMod+1]->Fill(); 
            }   
          }
        }
      }
    }
  }
  
  delete pslsfit;
  
  ADCFile->Close(); 
  
  // Remove temporary file
  //=======================
  std::stringstream del;
  del << "rm " <<ADCfile;
  system(del.str().c_str());
  
  
  // Create output trees
  //=====================

  resFile = new TFile(resfile.c_str(),"RECREATE");
  
  for (unsigned int iColor=0;iColor<nCol;iColor++){
    
    std::stringstream nametree;
    nametree <<"APDCol"<<colors[iColor];
    std::stringstream nametree2;
    nametree2 <<"PNCol"<<colors[iColor];
    
    restrees[iColor]= new TTree(nametree.str().c_str(),nametree.str().c_str());
    respntrees[iColor]= new TTree(nametree2.str().c_str(),nametree2.str().c_str());
    
    restrees[iColor]->Branch( "iphi",        &iphi,        "iphi/I"           );
    restrees[iColor]->Branch( "ieta",        &ieta,        "ieta/I"           );
    restrees[iColor]->Branch( "side",        &side,        "side/I"           );
    restrees[iColor]->Branch( "dccID",       &dccID,       "dccID/I"          );
    restrees[iColor]->Branch( "moduleID",    &moduleID,    "moduleID/I"       );
    restrees[iColor]->Branch( "towerID",     &towerID,     "towerID/I"        );
    restrees[iColor]->Branch( "channelID",   &channelID,   "channelID/I"      );
    restrees[iColor]->Branch( "APD",         &APD,         "APD[6]/D"         );
    restrees[iColor]->Branch( "Time",        &Time,        "Time[6]/D"        );
    restrees[iColor]->Branch( "APDoPN",      &APDoPN,      "APDoPN[6]/D"      );
    restrees[iColor]->Branch( "APDoPNA",     &APDoPNA,     "APDoPNA[6]/D"     );
    restrees[iColor]->Branch( "APDoPNB",     &APDoPNB,     "APDoPNB[6]/D"     );
    restrees[iColor]->Branch( "APDoAPDA",    &APDoAPDA,    "APDoAPDA[6]/D"    );
    restrees[iColor]->Branch( "APDoAPDB",    &APDoAPDB,    "APDoAPDB[6]/D"    );
    restrees[iColor]->Branch( "flag",        &flag,        "flag/I"           ); 

    respntrees[iColor]->Branch( "side",      &side,      "side/I"         );
    respntrees[iColor]->Branch( "moduleID",  &moduleID,  "moduleID/I"     );
    respntrees[iColor]->Branch( "pnID",      &pnID,      "pnID/I"         );
    respntrees[iColor]->Branch( "PN",        &PN,        "PN[6]/D"        ); 
    respntrees[iColor]->Branch( "PNoPN",     &PNoPN,     "PNoPN[6]/D"     );
    respntrees[iColor]->Branch( "PNoPNA",    &PNoPNA,    "PNoPNA[6]/D"    );
    respntrees[iColor]->Branch( "PNoPNB",    &PNoPNB,    "PNoPNB[6]/D"    );

    restrees[iColor]->SetBranchAddress( "iphi",        &iphi       );
    restrees[iColor]->SetBranchAddress( "ieta",        &ieta       );
    restrees[iColor]->SetBranchAddress( "side",        &side       );
    restrees[iColor]->SetBranchAddress( "dccID",       &dccID      );
    restrees[iColor]->SetBranchAddress( "moduleID",    &moduleID   );
    restrees[iColor]->SetBranchAddress( "towerID",     &towerID    );
    restrees[iColor]->SetBranchAddress( "channelID",   &channelID  );
    restrees[iColor]->SetBranchAddress( "APD",         APD         ); 
    restrees[iColor]->SetBranchAddress( "Time",        Time        );   
    restrees[iColor]->SetBranchAddress( "APDoPN",      APDoPN      ); 
    restrees[iColor]->SetBranchAddress( "APDoPNA",     APDoPNA     );
    restrees[iColor]->SetBranchAddress( "APDoPNB",     APDoPNB     );
    restrees[iColor]->SetBranchAddress( "APDoAPDA",    APDoAPDA    );
    restrees[iColor]->SetBranchAddress( "APDoAPDB",    APDoAPDB    );
    restrees[iColor]->SetBranchAddress( "flag",        &flag       ); 
   
    respntrees[iColor]->SetBranchAddress( "side",      &side     );
    respntrees[iColor]->SetBranchAddress( "moduleID",  &moduleID );
    respntrees[iColor]->SetBranchAddress( "pnID",      &pnID     );
    respntrees[iColor]->SetBranchAddress( "PN",        PN        );
    respntrees[iColor]->SetBranchAddress( "PNoPN",     PNoPN     );
    respntrees[iColor]->SetBranchAddress( "PNoPNA",    PNoPNA    );
    respntrees[iColor]->SetBranchAddress( "PNoPNB",    PNoPNB    );
  
  }

  // Set Cuts for PN stuff
  //=======================

  for (unsigned int iM=0;iM<nMod;iM++){
    unsigned int iMod=modules[iM]-1;
    
    for (unsigned int ich=0;ich<nPNPerMod;ich++){
      for (unsigned int icol=0;icol<nCol;icol++){
        PNAnal[iMod][ich][icol]->setPNCut(PNFirstAnal[iMod][ich][icol]->getPN().at(0),
                                          PNFirstAnal[iMod][ich][icol]->getPN().at(1));
      }
    }
  }
  
  // Build ref trees indexes
  //========================
  for(unsigned int imod=0;imod<nMod;imod++){
    int jmod=modules[imod];
    if( RefAPDtrees[0][jmod]->GetEntries()!=0 && RefAPDtrees[1][jmod]->GetEntries()!=0 ){
      RefAPDtrees[0][jmod]->BuildIndex("eventref");
      RefAPDtrees[1][jmod]->BuildIndex("eventref");
    }
  }
  
  
  // Final loop on crystals
  //=======================

  for (unsigned int iCry=0;iCry<nCrys;iCry++){ 
    
    unsigned int iMod=iModule[iCry]-1;
    
    // Set cuts on APD stuff
    //=======================
    
    for(unsigned int iCol=0;iCol<nCol;iCol++){   
      
      std::vector<double> lowcut;
      std::vector<double> highcut;
      double cutMin;
      double cutMax;

      cutMin=APDFirstAnal[iCry][iCol]->getAPD().at(0)-2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
      if(cutMin<0) cutMin=0;
      cutMax=APDFirstAnal[iCry][iCol]->getAPD().at(0)+2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
      
      lowcut.push_back(cutMin);
      highcut.push_back(cutMax);
      
      cutMin=APDFirstAnal[iCry][iCol]->getTime().at(0)-2.0*APDFirstAnal[iCry][iCol]->getTime().at(1);
      cutMax=APDFirstAnal[iCry][iCol]->getTime().at(0)+2.0*APDFirstAnal[iCry][iCol]->getTime().at(1); 
      lowcut.push_back(cutMin);
      highcut.push_back(cutMax);
      
      APDAnal[iCry][iCol]=new TAPD();
      APDAnal[iCry][iCol]->setAPDCut(APDFirstAnal[iCry][iCol]->getAPD().at(0),APDFirstAnal[iCry][iCol]->getAPD().at(1));
      APDAnal[iCry][iCol]->setAPDoPNCut(APDFirstAnal[iCry][iCol]->getAPDoPN().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN().at(1));
      APDAnal[iCry][iCol]->setAPDoPN0Cut(APDFirstAnal[iCry][iCol]->getAPDoPN0().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN0().at(1));
      APDAnal[iCry][iCol]->setAPDoPN1Cut(APDFirstAnal[iCry][iCol]->getAPDoPN1().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN1().at(1));
      APDAnal[iCry][iCol]->setTimeCut(APDFirstAnal[iCry][iCol]->getTime().at(0),APDFirstAnal[iCry][iCol]->getTime().at(1));
      APDAnal[iCry][iCol]->set2DAPDoAPD0Cut(lowcut,highcut);
      APDAnal[iCry][iCol]->set2DAPDoAPD1Cut(lowcut,highcut);
      
    }

    // Final loop on events
    //=======================

    Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry=0; jentry< APDtrees[iCry]->GetEntriesFast();jentry++) { 
      nb = APDtrees[iCry]->GetEntry(jentry);   nbytes += nb; 
      
      double pnmean;
      if (pn0<10 && pn1>10) {
        pnmean=pn1;
      }else if (pn1<10 && pn0>10){
        pnmean=pn0;
      }else pnmean=0.5*(pn0+pn1);
      
      // Get back color
      //================
      
      unsigned int iCol=0;
      for(unsigned int i=0;i<nCol;i++){
        if(color==colors[i]) {
          iCol=i;
          i=colors.size();
        }
      }
      
      // Fill PN stuff
      //===============
      
      if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){ 
        for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
          PNAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
        }
      }
      
      // Get ref amplitudes
      //===================

      if (_debug==1) std::cout <<"-- debug test -- Last Loop event:"<<event<<" apdAmpl:"<< apdAmpl<< std::endl;
      apdAmplA = 0.0;
      apdAmplB = 0.0;
      
      for (unsigned int iRef=0;iRef<nRefChan;iRef++){
        RefAPDtrees[iRef][iMod+1]->GetEntryWithIndex(event); 
      }
      
      if (_debug==1 ) std::cout <<"-- debug test -- Last Loop apdAmplA:"<<apdAmplA<< " apdAmplB:"<< apdAmplB<<", event:"<< event<<", eventref:"<< eventref<< std::endl;
      
      
      // Fill APD stuff
      //===============
      
      APDAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime, apdAmplA, apdAmplB); 
   
    }
    
    moduleID=iMod+1;
    
    if( moduleID>=20 ) moduleID-=2; // Trick to fix endcap specificity   
    
    // Get final results for APD
    //===========================
    
    for(unsigned int iColor=0;iColor<nCol;iColor++){
      
      std::vector<double> apdvec = APDAnal[iCry][iColor]->getAPD();
      std::vector<double> apdpnvec = APDAnal[iCry][iColor]->getAPDoPN();
      std::vector<double> apdpn0vec = APDAnal[iCry][iColor]->getAPDoPN0();
      std::vector<double> apdpn1vec = APDAnal[iCry][iColor]->getAPDoPN1();
      std::vector<double> timevec = APDAnal[iCry][iColor]->getTime();
      std::vector<double> apdapd0vec = APDAnal[iCry][iColor]->getAPDoAPD0();
      std::vector<double> apdapd1vec = APDAnal[iCry][iColor]->getAPDoAPD1();
      
      
      for(unsigned int i=0;i<apdvec.size();i++){
        
        APD[i]=apdvec.at(i);
        APDoPN[i]=apdpnvec.at(i);
        APDoPNA[i]=apdpn0vec.at(i);
        APDoPNB[i]=apdpn1vec.at(i);
        APDoAPDA[i]=apdapd0vec.at(i);
        APDoAPDB[i]=apdapd1vec.at(i);
        Time[i]=timevec.at(i);
      }
      
      
      // Fill APD results trees
      //========================

      iphi=iPhi[iCry];
      ieta=iEta[iCry];
      dccID=idccID[iCry];
      side=iside[iCry];
      towerID=iTowerID[iCry];
      channelID=iChannelID[iCry];

      
      if( !wasGainOK[iCry] || !wasTimingOK[iCry] || IsThereDataADC[iCry][iColor]==0 ){
        flag=0;
      }else flag=1;
      
      restrees[iColor]->Fill();   
    }
  }
   
  // Get final results for PN 
  //==========================
  
  for (unsigned int iM=0;iM<nMod;iM++){
    unsigned int iMod=modules[iM]-1;

    side=iside[firstChanMod[iMod]];
    
    for (unsigned int ch=0;ch<nPNPerMod;ch++){

      pnID=ch;
      moduleID=iMod+1;
      
      if( moduleID>=20 ) moduleID-=2;  // Trick to fix endcap specificity

      for(unsigned int iColor=0;iColor<nCol;iColor++){

        std::vector<double> pnvec = PNAnal[iMod][ch][iColor]->getPN();
        std::vector<double> pnopnvec = PNAnal[iMod][ch][iColor]->getPNoPN();
        std::vector<double> pnopn0vec = PNAnal[iMod][ch][iColor]->getPNoPN0();
        std::vector<double> pnopn1vec = PNAnal[iMod][ch][iColor]->getPNoPN1();
        
        for(unsigned int i=0;i<pnvec.size();i++){
          
          PN[i]=pnvec.at(i);
          PNoPN[i]=pnopnvec.at(i);
          PNoPNA[i]=pnopn0vec.at(i);
          PNoPNB[i]=pnopn1vec.at(i);
        }
        
        // Fill PN results trees
        //========================

        respntrees[iColor]->Fill();
      }
    }
  }
  
  // Remove temporary files
  //========================
  if(!_saveallevents){

    APDFile->Close();
    std::stringstream del2;
    del2 << "rm " <<APDfile;
    system(del2.str().c_str());

  }else {
    
    APDFile->cd();
    APDtrees[0]->Write();
    
    APDFile->Close();
    resFile->cd();
  }
  
  // Save results 
  //===============

  for (unsigned int i=0;i<nCol;i++){
    restrees[i]->Write();
    respntrees[i]->Write();
  }  
  
  resFile->Close(); 
  
  std::cout <<    "\t+=+    .................................................. done  +=+" << std::endl;
  std::cout <<    "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
}

void EcalLaserAnalyzer::setGeomEB	(	int	etaG,
		int	phiG,
		int	module,
		int	tower,
		int	strip,
		int	xtal,
		int	apdRefTT,
		int	channel,
		int	lmr
	)

Definition at line 1348 of file EcalLaserAnalyzer.cc.

References ME::apdRefChannels(), apdRefMap, channelID, prof2calltree::count, dccID, eta, firstChanMod, iChannelID, idccID, iEta, iModule, iPhi, isFirstChanModFilled, iside, iTowerID, module(), modules, nRefChan, phi, side, and towerID.

Referenced by analyze().

                                                                                                                                   {
  
  side=MEEBGeom::side(etaG,phiG);
  
  assert( module>=*min_element(modules.begin(),modules.end()) && module<=*max_element(modules.begin(),modules.end()) );
      
  eta = etaG;
  phi = phiG;
  channelID=5*(strip-1) + xtal-1; 
  towerID=tower;
  
  std::vector<int> apdRefChan=ME::apdRefChannels(module, lmr);      
  for (unsigned int iref=0;iref<nRefChan;iref++){
    if(channelID==apdRefChan[iref] && towerID==apdRefTT 
       && apdRefMap[iref].count(module)==0){
      apdRefMap[iref][module]=channel;
    }
  }
  
  if(isFirstChanModFilled[module-1]==0) {
    firstChanMod[module-1]=channel;
    isFirstChanModFilled[module-1]=1;
  } 
  
  iEta[channel]=eta;
  iPhi[channel]=phi;
  iModule[channel]= module ;
  iTowerID[channel]=towerID;
  iChannelID[channel]=channelID;
  idccID[channel]=dccID;
  iside[channel]=side;

}

void EcalLaserAnalyzer::setGeomEE	(	int	etaG,
		int	phiG,
		int	iX,
		int	iY,
		int	iZ,
		int	module,
		int	tower,
		int	ch,
		int	apdRefTT,
		int	channel,
		int	lmr
	)

Definition at line 1383 of file EcalLaserAnalyzer.cc.

References ME::apdRefChannels(), apdRefMap, channelID, prof2calltree::count, dccID, eta, firstChanMod, iChannelID, idccID, iEta, iModule, iPhi, isFirstChanModFilled, iside, iTowerID, module(), modules, nRefChan, phi, side, and towerID.

Referenced by analyze().

                                                                                                                                              {
  
  side=MEEEGeom::side(iX, iY, iZ);

  assert( module>=*min_element(modules.begin(),modules.end()) 
          && module<=*max_element(modules.begin(),modules.end()) );
  
  eta = etaG;
  phi = phiG;
  channelID=ch;
  towerID=tower;
  
  std::vector<int> apdRefChan=ME::apdRefChannels(module, lmr);      
  for (unsigned int iref=0;iref<nRefChan;iref++){
    if(channelID==apdRefChan[iref] && towerID==apdRefTT 
       && apdRefMap[iref].count(module)==0){
      apdRefMap[iref][module]=channel;
    }
  }
  
  if(isFirstChanModFilled[module-1]==0) {
    firstChanMod[module-1]=channel;
    isFirstChanModFilled[module-1]=1;
  } 
  
  iEta[channel]=eta;
  iPhi[channel]=phi;
  iModule[channel]= module ;
  iTowerID[channel]=towerID;
  iChannelID[channel]=channelID;
  idccID[channel]=dccID;
  iside[channel]=side;

}