Reference Manual

Go to the documentation of this file.
00001 /* 
00002  *  \class EcalLaserAnalyzer
00003  *
00004  *  $Date: 2012/02/09 10:07:36 $
00005  *  primary author: Julie Malcles - CEA/Saclay
00006  *  author: Gautier Hamel De Monchenault - CEA/Saclay
00007  */
00008 
00009 #include <TAxis.h>
00010 #include <TH1.h>
00011 #include <TProfile.h>
00012 #include <TTree.h>
00013 #include <TChain.h>
00014 #include <TFile.h>
00015 #include <TMath.h>
00016 
00017 #include "CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalLaserAnalyzer.h"
00018 
00019 #include <sstream>
00020 #include <iostream>
00021 #include <fstream>
00022 #include <iomanip>
00023 
00024 #include <FWCore/MessageLogger/interface/MessageLogger.h>
00025 #include <FWCore/Utilities/interface/Exception.h>
00026 
00027 #include <FWCore/Framework/interface/ESHandle.h>
00028 #include <FWCore/Framework/interface/EventSetup.h>
00029 
00030 #include <Geometry/EcalMapping/interface/EcalElectronicsMapping.h>
00031 #include <Geometry/EcalMapping/interface/EcalMappingRcd.h>
00032 
00033 #include <FWCore/Framework/interface/Event.h>
00034 #include <FWCore/Framework/interface/MakerMacros.h>
00035 #include <FWCore/ParameterSet/interface/ParameterSet.h>
00036 #include <DataFormats/EcalDigi/interface/EcalDigiCollections.h>
00037 #include <DataFormats/EcalDetId/interface/EcalElectronicsId.h>
00038 #include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
00039 #include <DataFormats/EcalRawData/interface/EcalRawDataCollections.h>
00040 
00041 #include <vector>
00042 
00043 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TShapeAnalysis.h>
00044 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h>
00045 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
00046 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TAPD.h>
00047 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TAPDPulse.h>
00048 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPN.h>
00049 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNPulse.h>
00050 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNCor.h>
00051 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMem.h>
00052 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h>
00053 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h>
00054 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/MEGeom.h>
00055 
00056 
00057 //========================================================================
00058 EcalLaserAnalyzer::EcalLaserAnalyzer(const edm::ParameterSet& iConfig)
00059 //========================================================================
00060   :
00061 iEvent(0),
00062 
00063 // Framework parameters with default values
00064 
00065 _nsamples(        iConfig.getUntrackedParameter< unsigned int >( "nSamples",           10 ) ),
00066 _presample(       iConfig.getUntrackedParameter< unsigned int >( "nPresamples",         2 ) ),
00067 _firstsample(     iConfig.getUntrackedParameter< unsigned int >( "firstSample",         1 ) ),
00068 _lastsample(      iConfig.getUntrackedParameter< unsigned int >( "lastSample",          2 ) ),
00069 _nsamplesPN(      iConfig.getUntrackedParameter< unsigned int >( "nSamplesPN",         50 ) ), 
00070 _presamplePN(     iConfig.getUntrackedParameter< unsigned int >( "nPresamplesPN",       6 ) ),
00071 _firstsamplePN(   iConfig.getUntrackedParameter< unsigned int >( "firstSamplePN",       7 ) ),
00072 _lastsamplePN(    iConfig.getUntrackedParameter< unsigned int >( "lastSamplePN",        8 ) ),
00073 _timingcutlow(    iConfig.getUntrackedParameter< unsigned int >( "timingCutLow",        2 ) ),
00074 _timingcuthigh(   iConfig.getUntrackedParameter< unsigned int >( "timingCutHigh",       9 ) ),
00075 _timingquallow(   iConfig.getUntrackedParameter< unsigned int >( "timingQualLow",       3 ) ),
00076 _timingqualhigh(  iConfig.getUntrackedParameter< unsigned int >( "timingQualHigh",      8 ) ),
00077 _ratiomincutlow(  iConfig.getUntrackedParameter< double       >( "ratioMinCutLow",    0.4 ) ),
00078 _ratiomincuthigh( iConfig.getUntrackedParameter< double       >( "ratioMinCutHigh",  0.95 ) ),
00079 _ratiomaxcutlow(  iConfig.getUntrackedParameter< double       >( "ratioMaxCutLow",    0.8 ) ),
00080 _presamplecut(    iConfig.getUntrackedParameter< double       >( "presampleCut",      5.0 ) ),
00081 _niter(           iConfig.getUntrackedParameter< unsigned int >( "nIter",               3 ) ),
00082 _fitab(           iConfig.getUntrackedParameter< bool         >( "fitAB",           false ) ),
00083 _alpha(           iConfig.getUntrackedParameter< double       >( "alpha",       1.5076494 ) ),
00084 _beta(            iConfig.getUntrackedParameter< double       >( "beta",        1.5136036 ) ),
00085 _nevtmax(         iConfig.getUntrackedParameter< unsigned int >( "nEvtMax",           200 ) ),
00086 _noise(           iConfig.getUntrackedParameter< double       >( "noise",             2.0 ) ),
00087 _chi2cut(         iConfig.getUntrackedParameter< double       >( "chi2cut",          10.0 ) ), 
00088 _ecalPart(        iConfig.getUntrackedParameter< std::string  >( "ecalPart",         "EB" ) ),
00089 _docorpn(         iConfig.getUntrackedParameter< bool         >( "doCorPN",         false ) ),
00090 _fedid(           iConfig.getUntrackedParameter< int          >( "fedID",            -999 ) ),
00091 _saveallevents(   iConfig.getUntrackedParameter< bool         >( "saveAllEvents",   false ) ),
00092 _qualpercent(     iConfig.getUntrackedParameter< double       >( "qualPercent",       0.2 ) ),
00093 _debug(           iConfig.getUntrackedParameter< int          >( "debug",              0  ) ),
00094 nCrys(                                                                               NCRYSEB),
00095 nPNPerMod(                                                                         NPNPERMOD),
00096 nMod(                                                                                 NMODEE),
00097 nSides(                                                                               NSIDES),
00098 runType(-1), runNum(0), fedID(-1), dccID(-1), side(2), lightside(2), iZ(1),
00099 phi(-1), eta(-1), event(0), color(-1),pn0(0), pn1(0), apdAmpl(0), apdAmplA(0), 
00100 apdAmplB(0),apdTime(0),pnAmpl(0),
00101 pnID(-1), moduleID(-1), channelIteratorEE(0)
00102 
00103 
00104   //========================================================================
00105 
00106 
00107 {
00108 
00109   // Initialization from cfg file
00110 
00111   resdir_                 = iConfig.getUntrackedParameter<std::string>("resDir");
00112   pncorfile_              = iConfig.getUntrackedParameter<std::string>("pnCorFile");
00113 
00114   digiCollection_         = iConfig.getParameter<std::string>("digiCollection");
00115   digiPNCollection_       = iConfig.getParameter<std::string>("digiPNCollection");
00116   digiProducer_           = iConfig.getParameter<std::string>("digiProducer");
00117 
00118   eventHeaderCollection_  = iConfig.getParameter<std::string>("eventHeaderCollection");
00119   eventHeaderProducer_    = iConfig.getParameter<std::string>("eventHeaderProducer");
00120 
00121 
00122   // Geometrical constants initialization 
00123 
00124   if (_ecalPart == "EB") {
00125     nCrys    = NCRYSEB;
00126   } else {
00127     nCrys    = NCRYSEE;
00128   }
00129   iZ         =  1;
00130   if(_fedid <= 609 ) 
00131     iZ       = -1;
00132   modules    = ME::lmmodFromDcc(_fedid);
00133   nMod       = modules.size(); 
00134   nRefChan   = NREFCHAN;
00135   
00136   for(unsigned int j=0;j<nCrys;j++){
00137     iEta[j]=-1;
00138     iPhi[j]=-1;
00139     iModule[j]=10;
00140     iTowerID[j]=-1;
00141     iChannelID[j]=-1;
00142     idccID[j]=-1;
00143     iside[j]=-1;
00144     wasTimingOK[j]=true;
00145     wasGainOK[j]=true;
00146   }
00147   
00148   for(unsigned int j=0;j<nMod;j++){
00149     int ii= modules[j];
00150     firstChanMod[ii-1]=0;
00151     isFirstChanModFilled[ii-1]=0;
00152   }
00153   
00154 
00155   // Quality check flags
00156   
00157   isGainOK=true;
00158   isTimingOK=true;
00159 
00160   // PN linearity corrector
00161 
00162   pnCorrector = new TPNCor(pncorfile_.c_str());
00163 
00164 
00165   // Objects dealing with pulses
00166 
00167   APDPulse = new TAPDPulse(_nsamples, _presample, _firstsample, _lastsample, 
00168                            _timingcutlow, _timingcuthigh, _timingquallow, _timingqualhigh,
00169                            _ratiomincutlow,_ratiomincuthigh, _ratiomaxcutlow);
00170   PNPulse = new TPNPulse(_nsamplesPN, _presamplePN);
00171 
00172   // Object dealing with MEM numbering
00173 
00174   Mem = new TMem(_fedid);
00175 
00176   // Objects needed for npresample calculation
00177 
00178   Delta01=new TMom();
00179   Delta12=new TMom();
00180 
00181 }
00182 
00183 //========================================================================
00184 EcalLaserAnalyzer::~EcalLaserAnalyzer(){
00185   //========================================================================
00186 
00187 
00188   // do anything here that needs to be done at destruction time
00189   // (e.g. close files, deallocate resources etc.)
00190 
00191 }
00192 
00193 
00194 
00195 //========================================================================
00196 void EcalLaserAnalyzer::beginJob() {
00197   //========================================================================
00198 
00199   
00200   // Create temporary files and trees to save adc samples
00201   //======================================================
00202 
00203   ADCfile=resdir_;
00204   ADCfile+="/APDSamplesLaser.root";
00205 
00206   APDfile=resdir_;
00207   APDfile+="/APDPNLaserAllEvents.root";
00208 
00209   ADCFile = new TFile(ADCfile.c_str(),"RECREATE");
00210   
00211   for (unsigned int i=0;i<nCrys;i++){
00212     
00213     std::stringstream name;
00214     name << "ADCTree" <<i+1;
00215     ADCtrees[i]= new TTree(name.str().c_str(),name.str().c_str());
00216 
00217     ADCtrees[i]->Branch( "ieta",        &eta,         "eta/I"         );
00218     ADCtrees[i]->Branch( "iphi",        &phi,         "phi/I"         );
00219     ADCtrees[i]->Branch( "side",        &side,        "side/I"        );
00220     ADCtrees[i]->Branch( "dccID",       &dccID,       "dccID/I"       );
00221     ADCtrees[i]->Branch( "towerID",     &towerID,     "towerID/I"     );
00222     ADCtrees[i]->Branch( "channelID",   &channelID,   "channelID/I"   );
00223     ADCtrees[i]->Branch( "event",       &event,       "event/I"       );
00224     ADCtrees[i]->Branch( "color",       &color,       "color/I"       );
00225     ADCtrees[i]->Branch( "adc",         &adc ,        "adc[10]/D"     );
00226     ADCtrees[i]->Branch( "pn0",         &pn0 ,        "pn0/D"         );
00227     ADCtrees[i]->Branch( "pn1",         &pn1 ,        "pn1/D"         );
00228 
00229      
00230     ADCtrees[i]->SetBranchAddress( "ieta",        &eta         );  
00231     ADCtrees[i]->SetBranchAddress( "iphi",        &phi         ); 
00232     ADCtrees[i]->SetBranchAddress( "side",        &side        ); 
00233     ADCtrees[i]->SetBranchAddress( "dccID",       &dccID       ); 
00234     ADCtrees[i]->SetBranchAddress( "towerID",     &towerID     ); 
00235     ADCtrees[i]->SetBranchAddress( "channelID",   &channelID   ); 
00236     ADCtrees[i]->SetBranchAddress( "event",       &event       );   
00237     ADCtrees[i]->SetBranchAddress( "color",       &color       );   
00238     ADCtrees[i]->SetBranchAddress( "adc",         adc          );
00239     ADCtrees[i]->SetBranchAddress( "pn0",         &pn0         );
00240     ADCtrees[i]->SetBranchAddress( "pn1",         &pn1         );    
00241  
00242     nevtAB[i]=0 ;
00243   } 
00244 
00245 
00246   // Define output results filenames and shape analyzer object (alpha,beta) 
00247   //=====================================================================
00248   
00249   // 1) AlphaBeta files
00250   
00251   doesABTreeExist=true;
00252 
00253   std::stringstream nameabinitfile;
00254   nameabinitfile << resdir_ <<"/ABInit.root";
00255   alphainitfile=nameabinitfile.str();
00256   
00257   std::stringstream nameabfile;  
00258   nameabfile << resdir_ <<"/AB.root";
00259   alphafile=nameabfile.str();
00260   
00261   FILE *test;
00262   if(_fitab)
00263     test = fopen(alphainitfile.c_str(),"r"); 
00264   else 
00265     test = fopen(alphafile.c_str(),"r"); 
00266   if(test == NULL) {
00267     doesABTreeExist=false;
00268     _fitab=true;
00269   };
00270   delete test;
00271   
00272   TFile *fAB=0; TTree *ABInit=0;
00273   if(doesABTreeExist){
00274     fAB=new TFile(nameabinitfile.str().c_str());
00275   }
00276   if(doesABTreeExist && fAB){
00277     ABInit = (TTree*) fAB->Get("ABCol0");
00278   }
00279   
00280   // 2) Shape analyzer
00281   
00282   if(doesABTreeExist && fAB && ABInit && ABInit->GetEntries()!=0){
00283     shapana= new TShapeAnalysis(ABInit, _alpha, _beta, 5.5, 1.0);
00284     doesABTreeExist=true;
00285   }else{
00286     shapana= new TShapeAnalysis(_alpha, _beta, 5.5, 1.0);
00287     doesABTreeExist=false;
00288     _fitab=true;  
00289   }  
00290   shapana -> set_const(_nsamples,_firstsample,_lastsample,
00291                        _presample, _nevtmax, _noise, _chi2cut);
00292   
00293   if(doesABTreeExist && fAB ) fAB->Close();
00294   
00295 
00296   //  2) APD file
00297   
00298   std::stringstream nameapdfile;
00299   nameapdfile << resdir_ <<"/APDPN_LASER.root";      
00300   resfile=nameapdfile.str();
00301 
00302   
00303   // Laser events counter
00304 
00305   laserEvents=0;
00306 
00307 }
00308 
00309 //========================================================================
00310 void EcalLaserAnalyzer:: analyze( const edm::Event & e, const  edm::EventSetup& c){
00311 //========================================================================
00312 
00313   ++iEvent;
00314   
00315   // retrieving DCC header
00316 
00317   edm::Handle<EcalRawDataCollection> pDCCHeader;
00318   const  EcalRawDataCollection* DCCHeader=0;
00319   try {
00320     e.getByLabel(eventHeaderProducer_,eventHeaderCollection_, pDCCHeader);
00321     DCCHeader=pDCCHeader.product();
00322   }catch ( std::exception& ex ) {
00323     std::cerr << "Error! can't get the product retrieving DCC header" << 
00324       eventHeaderCollection_.c_str() <<" "<< eventHeaderProducer_.c_str() << std::endl;
00325   }
00326  
00327   //retrieving crystal data from Event
00328 
00329   edm::Handle<EBDigiCollection>  pEBDigi;
00330   const  EBDigiCollection* EBDigi=0;
00331   edm::Handle<EEDigiCollection>  pEEDigi;
00332   const  EEDigiCollection* EEDigi=0;
00333 
00334 
00335   if (_ecalPart == "EB") {
00336     try {
00337       e.getByLabel(digiProducer_,digiCollection_, pEBDigi); 
00338       EBDigi=pEBDigi.product(); 
00339     }catch ( std::exception& ex ) {
00340       std::cerr << "Error! can't get the product retrieving EB crystal data " << 
00341         digiCollection_.c_str() << std::endl;
00342     } 
00343   } else if (_ecalPart == "EE") {
00344     try {
00345       e.getByLabel(digiProducer_,digiCollection_, pEEDigi); 
00346       EEDigi=pEEDigi.product(); 
00347     }catch ( std::exception& ex ) {
00348       std::cerr << "Error! can't get the product retrieving EE crystal data " << 
00349         digiCollection_.c_str() << std::endl;
00350     } 
00351   } else {
00352     std::cout <<" Wrong ecalPart in cfg file " << std::endl;
00353     return;
00354   }
00355   
00356   // retrieving crystal PN diodes from Event
00357   
00358   edm::Handle<EcalPnDiodeDigiCollection>  pPNDigi;
00359   const  EcalPnDiodeDigiCollection* PNDigi=0;
00360   try {
00361     e.getByLabel(digiProducer_, pPNDigi);
00362     PNDigi=pPNDigi.product(); 
00363   }catch ( std::exception& ex ) {
00364     std::cerr << "Error! can't get the product " << digiPNCollection_.c_str() << std::endl;
00365   }
00366 
00367   // retrieving electronics mapping
00368 
00369   edm::ESHandle< EcalElectronicsMapping > ecalmapping;
00370   const EcalElectronicsMapping* TheMapping=0; 
00371   try{
00372     c.get< EcalMappingRcd >().get(ecalmapping);
00373     TheMapping = ecalmapping.product();
00374   }catch ( std::exception& ex ) {
00375     std::cerr << "Error! can't get the product EcalMappingRcd"<< std::endl;
00376   }
00377 
00378 
00379   // ============================
00380   // Decode DCCHeader Information 
00381   // ============================
00382   
00383   for ( EcalRawDataCollection::const_iterator headerItr= DCCHeader->begin();
00384         headerItr != DCCHeader->end(); ++headerItr ) {
00385     
00386     // Get run type and run number 
00387 
00388     int fed = headerItr->fedId();  
00389     if(fed!=_fedid && _fedid!=-999) continue; 
00390     
00391     runType=headerItr->getRunType();
00392     runNum=headerItr->getRunNumber();
00393     event=headerItr->getLV1();
00394 
00395     dccID=headerItr->getDccInTCCCommand();
00396     fedID=headerItr->fedId();  
00397     lightside=headerItr->getRtHalf();
00398 
00399     // Check fed corresponds to the DCC in TCC
00400     
00401     if( 600+dccID != fedID ) continue;
00402 
00403     // Cut on runType
00404 
00405     if( runType!=EcalDCCHeaderBlock::LASER_STD 
00406        && runType!=EcalDCCHeaderBlock::LASER_GAP 
00407        && runType!=EcalDCCHeaderBlock::LASER_POWER_SCAN
00408        && runType!=EcalDCCHeaderBlock::LASER_DELAY_SCAN ) return; 
00409     
00410     // Retrieve laser color and event number
00411     
00412     EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
00413     color = settings.wavelength;
00414     if( color<0 ) return;
00415     
00416     std::vector<int>::iterator iter = find( colors.begin(), colors.end(), color );
00417     if( iter==colors.end() ){
00418       colors.push_back( color );
00419       std::cout <<" new color found "<< color<<" "<< colors.size()<< std::endl;
00420     }
00421   }
00422 
00423   // Cut on fedID
00424  
00425   if(fedID!=_fedid && _fedid!=-999) return; 
00426     
00427   // Count laser events
00428 
00429   laserEvents++;
00430 
00431 
00432   // ======================
00433   // Decode PN Information
00434   // ======================
00435   
00436   TPNFit * pnfit = new TPNFit();
00437   pnfit -> init(_nsamplesPN,_firstsamplePN,  _lastsamplePN);
00438   
00439   double chi2pn=0;
00440   unsigned int samplemax=0;
00441   int  pnGain=0;
00442   
00443   std::map <int, std::vector<double> > allPNAmpl;
00444   std::map <int, std::vector<double> > allPNGain;
00445 
00446   
00447   // Loop on PNs digis
00448   
00449   for ( EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin();
00450         pnItr != PNDigi->end(); ++pnItr ) { 
00451     
00452     EcalPnDiodeDetId pnDetId = EcalPnDiodeDetId((*pnItr).id());
00453     
00454     if (_debug==1) std::cout <<"-- debug -- Inside PNDigi - pnID=" <<
00455                      pnDetId.iPnId()<<", dccID="<< pnDetId.iDCCId()<< std::endl;
00456 
00457     // Skip MEM DCC without relevant data
00458   
00459     bool isMemRelevant=Mem->isMemRelevant(pnDetId.iDCCId());
00460     if(!isMemRelevant) continue;
00461 
00462     // Loop on PN samples
00463 
00464     for ( int samId=0; samId < (*pnItr).size() ; samId++ ) {   
00465       pn[samId]=(*pnItr).sample(samId).adc();  
00466       pnG[samId]=(*pnItr).sample(samId).gainId(); 
00467       if (samId==0) pnGain=pnG[samId];
00468       if (samId>0) pnGain=int(TMath::Max(pnG[samId],pnGain));     
00469     }
00470     
00471     if( pnGain!=1 ) std::cout << "PN gain different from 1"<< std::endl;
00472     
00473     // Calculate amplitude from pulse
00474     
00475     PNPulse->setPulse(pn);
00476     pnNoPed=PNPulse->getAdcWithoutPedestal();
00477     samplemax=PNPulse->getMaxSample();
00478     chi2pn = pnfit -> doFit(samplemax,&pnNoPed[0]); 
00479     if(chi2pn == 101 || chi2pn == 102 || chi2pn == 103) pnAmpl=0.;
00480     else pnAmpl= pnfit -> getAmpl();
00481     
00482     // Apply linearity correction
00483 
00484     double corr=1.0;
00485     if( _docorpn ) corr=pnCorrector->getPNCorrectionFactor(pnAmpl, pnGain);
00486     pnAmpl*=corr;
00487 
00488     // Fill PN ampl vector
00489 
00490     allPNAmpl[pnDetId.iDCCId()].push_back(pnAmpl);
00491       
00492     if (_debug==1) std::cout <<"-- debug -- Inside PNDigi - PNampl=" << 
00493                      pnAmpl<<", PNgain="<< pnGain<<std::endl;  
00494   }
00495   
00496   // ===========================
00497   // Decode EBDigis Information
00498   // ===========================
00499   
00500   int  adcGain=0;
00501 
00502   if (EBDigi){
00503     
00504     // Loop on crystals
00505     //=================== 
00506     
00507     for ( EBDigiCollection::const_iterator digiItr= EBDigi->begin();
00508           digiItr != EBDigi->end(); ++digiItr ) { 
00509       
00510       // Retrieve geometry
00511       //=================== 
00512       
00513       EBDetId id_crystal(digiItr->id()) ;
00514       EBDataFrame df( *digiItr );
00515       EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
00516       
00517       int etaG = id_crystal.ieta() ;  // global
00518       int phiG = id_crystal.iphi() ;  // global
00519       
00520       std::pair<int, int> LocalCoord=MEEBGeom::localCoord( etaG , phiG );
00521       
00522       int etaL=LocalCoord.first ; // local
00523       int phiL=LocalCoord.second ;// local
00524       
00525       int strip=elecid_crystal.stripId();
00526       int xtal=elecid_crystal.xtalId();
00527       
00528       int module= MEEBGeom::lmmod(etaG, phiG);
00529       int tower=elecid_crystal.towerId();
00530            
00531       int apdRefTT=MEEBGeom::apdRefTower(module);      
00532            
00533       std::pair<int,int> pnpair=MEEBGeom::pn(module);
00534       unsigned int MyPn0=pnpair.first;
00535       unsigned int MyPn1=pnpair.second;
00536       
00537       int lmr=MEEBGeom::lmr( etaG,phiG ); 
00538       unsigned int channel=MEEBGeom::electronic_channel( etaL, phiL );
00539       assert( channel < nCrys );
00540 
00541       setGeomEB(etaG, phiG, module, tower, strip, xtal, apdRefTT, channel, lmr);      
00542       
00543       if (_debug==1) std::cout << "-- debug -- Inside EBDigi - towerID:"<< towerID<<
00544                        " channelID:" <<channelID<<" module:"<< module<<
00545                        " modules:"<<modules.size()<< std::endl;
00546       
00547       // APD Pulse
00548       //=========== 
00549 
00550       // Loop on adc samples 
00551      
00552       for (unsigned int i=0; i< (*digiItr).size() ; ++i ) {   
00553 
00554         EcalMGPASample samp_crystal(df.sample(i));
00555         adc[i]=samp_crystal.adc() ;    
00556         adcG[i]=samp_crystal.gainId();   
00557         adc[i]*=adcG[i];
00558         if (i==0) adcGain=adcG[i];
00559         if (i>0) adcGain=TMath::Max(adcG[i],adcGain);  
00560       }
00561 
00562       APDPulse->setPulse(adc);
00563       
00564       // Quality checks
00565       //================
00566       
00567       if(adcGain!=1) nEvtBadGain[channel]++;   
00568       if(!APDPulse->isTimingQualOK()) nEvtBadTiming[channel]++;
00569       nEvtTot[channel]++;
00570 
00571       // Associate PN ampl
00572       //===================
00573  
00574       int mem0=Mem->Mem(lmr,0);
00575       int mem1=Mem->Mem(lmr,1);
00576 
00577       if(allPNAmpl[mem0].size()>MyPn0) pn0=allPNAmpl[mem0][MyPn0];
00578       else pn0=0;
00579       if(allPNAmpl[mem1].size()>MyPn1) pn1=allPNAmpl[mem1][MyPn1];
00580       else pn1=0;
00581 
00582 
00583       // Fill if Pulse is fine
00584       //=======================
00585       
00586       if( APDPulse->isPulseOK() && lightside==side){
00587         
00588         ADCtrees[channel]->Fill(); 
00589         
00590         Delta01->addEntry(APDPulse->getDelta(0,1));
00591         Delta12->addEntry(APDPulse->getDelta(1,2));
00592         
00593         if( nevtAB[channel] < _nevtmax && _fitab ){
00594           if(doesABTreeExist)  shapana -> putAllVals(channel, adc, eta, phi);   
00595           else shapana -> putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);
00596           nevtAB[channel]++ ;
00597         }
00598       }
00599     }
00600     
00601   } else if (EEDigi) {
00602 
00603     // Loop on crystals
00604     //===================
00605     
00606     for ( EEDigiCollection::const_iterator digiItr= EEDigi->begin(); 
00607           digiItr != EEDigi->end(); ++digiItr ) {  
00608       
00609       // Retrieve geometry
00610       //=================== 
00611       
00612       EEDetId id_crystal(digiItr->id()) ;
00613       EEDataFrame df( *digiItr );
00614       EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
00615        
00616       int etaG = id_crystal.iy() ; 
00617       int phiG = id_crystal.ix() ;
00618 
00619       int iX = (phiG-1)/5+1;
00620       int iY = (etaG-1)/5+1;
00621 
00622       int tower=elecid_crystal.towerId();
00623       int ch=elecid_crystal.channelId()-1; 
00624 
00625       int module=MEEEGeom::lmmod( iX, iY );
00626       if( module>=18 && side==1 ) module+=2;
00627       int lmr=MEEEGeom::lmr( iX, iY ,iZ); 
00628       int dee=MEEEGeom::dee(lmr);
00629       int apdRefTT=MEEEGeom::apdRefTower(lmr, module);  
00630       
00631       std::pair<int,int> pnpair=MEEEGeom::pn( dee, module ) ; 
00632       unsigned int MyPn0=pnpair.first;
00633       unsigned int MyPn1=pnpair.second;
00634 
00635       int hashedIndex=100000*eta+phi;
00636       if( channelMapEE.count(hashedIndex) == 0 ){
00637         channelMapEE[hashedIndex]=channelIteratorEE;
00638         channelIteratorEE++;
00639       }
00640       unsigned int channel=channelMapEE[hashedIndex];
00641       assert ( channel < nCrys );
00642 
00643       setGeomEE(etaG, phiG, iX, iY, iZ, module, tower, ch, apdRefTT, channel, lmr);
00644       
00645       
00646       if (_debug==1) std::cout << "-- debug -- Inside EEDigi - towerID:"<< towerID<<
00647                        " channelID:" <<channelID<<" module:"<< module<<
00648                        " modules:"<<modules.size()<< std::endl;
00649       
00650       // APD Pulse
00651       //=========== 
00652 
00653       if( (*digiItr).size()>10) std::cout <<"SAMPLES SIZE > 10!" <<  (*digiItr).size()<< std::endl;
00654  
00655       // Loop on adc samples  
00656 
00657       for (unsigned int i=0; i< (*digiItr).size() ; ++i ) { 
00658 
00659         EcalMGPASample samp_crystal(df.sample(i));
00660         adc[i]=samp_crystal.adc() ;    
00661         adcG[i]=samp_crystal.gainId();   
00662         adc[i]*=adcG[i];
00663         
00664         if (i==0) adcGain=adcG[i];
00665         if (i>0) adcGain=TMath::Max(adcG[i],adcGain);  
00666       }
00667       
00668       APDPulse->setPulse(adc);
00669       
00670       // Quality checks
00671       //================
00672       
00673       if(adcGain!=1) nEvtBadGain[channel]++;   
00674       if(!APDPulse->isTimingQualOK()) nEvtBadTiming[channel]++;
00675       nEvtTot[channel]++;
00676 
00677       
00678       // Associate PN ampl
00679       //===================
00680 
00681       int mem0=Mem->Mem(lmr,0);
00682       int mem1=Mem->Mem(lmr,1);
00683       
00684       if(allPNAmpl[mem0].size()>MyPn0) pn0=allPNAmpl[mem0][MyPn0];
00685       else pn0=0;
00686       if(allPNAmpl[mem1].size()>MyPn1) pn1=allPNAmpl[mem1][MyPn1];
00687       else pn1=0;
00688 
00689       // Fill if Pulse is fine
00690       //=======================
00691       
00692       if( APDPulse->isPulseOK() && lightside==side){
00693         ADCtrees[channel]->Fill(); 
00694         
00695         Delta01->addEntry(APDPulse->getDelta(0,1));
00696         Delta12->addEntry(APDPulse->getDelta(1,2));
00697         
00698         if( nevtAB[channel] < _nevtmax && _fitab ){
00699           if(doesABTreeExist)  shapana -> putAllVals(channel, adc, eta, phi);   
00700           else shapana -> putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);  
00701           nevtAB[channel]++ ;
00702         }
00703       }
00704     }
00705   }
00706 }
00707 // analyze
00708 
00709 
00710 //========================================================================
00711 void EcalLaserAnalyzer::endJob() {
00712 //========================================================================
00713 
00714   // Adjust channel numbers for EE 
00715   //===============================
00716 
00717   if( _ecalPart == "EE" ) {
00718     nCrys=channelMapEE.size();
00719     shapana->set_nch(nCrys);
00720   }
00721   
00722   // Set presamples number 
00723   //======================
00724   
00725   double delta01=Delta01->getMean();
00726   double delta12=Delta12->getMean();
00727   if(delta12>_presamplecut) {
00728     _presample=2;
00729     if(delta01>_presamplecut) _presample=1;
00730   }
00731   
00732   APDPulse->setPresamples(_presample);
00733   shapana->set_presample(_presample);
00734 
00735   //  Get alpha and beta
00736   //======================
00737 
00738   if(_fitab){
00739     std::cout <<  "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
00740     std::cout <<   "\t+=+     Analyzing data: getting (alpha, beta)     +=+" << std::endl;
00741     TFile *fAB=0; TTree *ABInit=0;
00742     if(doesABTreeExist){
00743       fAB=new TFile(alphainitfile.c_str());
00744     }
00745     if(doesABTreeExist && fAB){
00746       ABInit = (TTree*) fAB->Get("ABCol0");
00747     }
00748     shapana->computeShape(alphafile, ABInit);
00749     std::cout <<    "\t+=+    .................................... done  +=+" << std::endl;
00750     std::cout <<    "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
00751   }
00752   
00753 
00754   // Don't do anything if there is no events
00755   //=========================================
00756 
00757   if( laserEvents == 0 ){
00758     
00759     ADCFile->Close(); 
00760     std::stringstream del;
00761     del << "rm " <<ADCfile;
00762     system(del.str().c_str());
00763     std::cout << " No Laser Events "<< std::endl;
00764     return;
00765   }
00766 
00767   // Set quality flags for gains and timing
00768   //=========================================
00769 
00770   double BadGainEvtPercentage=0.0;
00771   double BadTimingEvtPercentage=0.0;
00772   
00773   int nChanBadGain=0;
00774   int nChanBadTiming=0;
00775   
00776   for (unsigned int i=0;i<nCrys;i++){
00777     if(nEvtTot[i]!=0){
00778       BadGainEvtPercentage=double(nEvtBadGain[i])/double(nEvtTot[i]);
00779       BadTimingEvtPercentage=double(nEvtBadTiming[i])/double(nEvtTot[i]);
00780     }
00781     if(BadGainEvtPercentage>_qualpercent) {
00782       wasGainOK[i]=false;
00783       nChanBadGain++;
00784     }
00785     if(BadTimingEvtPercentage>_qualpercent){
00786       wasTimingOK[i]=false;
00787       nChanBadTiming++;
00788     }
00789   }
00790   
00791   double BadGainChanPercentage=double(nChanBadGain)/double(nCrys);
00792   double BadTimingChanPercentage=double(nChanBadTiming)/double(nCrys);
00793   
00794   if(BadGainChanPercentage>_qualpercent) isGainOK = false;
00795   if(BadTimingChanPercentage>_qualpercent) isTimingOK = false;
00796 
00797   // Analyze adc samples to get amplitudes
00798   //=======================================
00799   
00800   std::cout <<  "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
00801   std::cout <<    "\t+=+     Analyzing laser data: getting APD, PN, APD/PN, PN/PN    +=+" << std::endl;
00802   
00803   if( !isGainOK )
00804   std::cout <<    "\t+=+ ............................ WARNING! APD GAIN WAS NOT 1    +=+" << std::endl;
00805   if( !isTimingOK )
00806   std::cout <<    "\t+=+ ............................ WARNING! TIMING WAS BAD        +=+" << std::endl;
00807     
00808   
00809   APDFile = new TFile(APDfile.c_str(),"RECREATE");
00810   
00811   int ieta, iphi, channelID, towerID, flag;
00812 
00813   for (unsigned int i=0;i<nCrys;i++){
00814     
00815     std::stringstream name;
00816     name << "APDTree" <<i+1;
00817     
00818     APDtrees[i]= new TTree(name.str().c_str(),name.str().c_str());
00819     
00820     //List of branches
00821     
00822     APDtrees[i]->Branch( "event",     &event,     "event/I"     );
00823     APDtrees[i]->Branch( "color",     &color,     "color/I"     );
00824     APDtrees[i]->Branch( "iphi",      &iphi,      "iphi/I"      );
00825     APDtrees[i]->Branch( "ieta",      &ieta,      "ieta/I"      );
00826     APDtrees[i]->Branch( "side",      &side,      "side/I"      );
00827     APDtrees[i]->Branch( "dccID",     &dccID,     "dccID/I"     );
00828     APDtrees[i]->Branch( "towerID",   &towerID,   "towerID/I"   );
00829     APDtrees[i]->Branch( "channelID", &channelID, "channelID/I" );
00830     APDtrees[i]->Branch( "apdAmpl",   &apdAmpl,   "apdAmpl/D"   );
00831     APDtrees[i]->Branch( "apdTime",   &apdTime,   "apdTime/D"   );
00832     if(_saveallevents) APDtrees[i]->Branch( "adc", &adc ,"adc[10]/D" );
00833     APDtrees[i]->Branch( "flag",      &flag,      "flag/I"      ); 
00834     APDtrees[i]->Branch( "flagAB",    &flagAB,    "flagAB/I"    ); 
00835     APDtrees[i]->Branch( "pn0",       &pn0,       "pn0/D"       );
00836     APDtrees[i]->Branch( "pn1",       &pn1,       "pn1/D"       ); 
00837     
00838     APDtrees[i]->SetBranchAddress( "event",     &event     );
00839     APDtrees[i]->SetBranchAddress( "color",     &color     ); 
00840     APDtrees[i]->SetBranchAddress( "iphi",      &iphi      );
00841     APDtrees[i]->SetBranchAddress( "ieta",      &ieta      );
00842     APDtrees[i]->SetBranchAddress( "side",      &side      );
00843     APDtrees[i]->SetBranchAddress( "dccID",     &dccID     );
00844     APDtrees[i]->SetBranchAddress( "towerID",   &towerID   );
00845     APDtrees[i]->SetBranchAddress( "channelID", &channelID ); 
00846     APDtrees[i]->SetBranchAddress( "apdAmpl",   &apdAmpl   );
00847     APDtrees[i]->SetBranchAddress( "apdTime",   &apdTime   );
00848     if(_saveallevents)APDtrees[i]->SetBranchAddress( "adc", adc );
00849     APDtrees[i]->SetBranchAddress( "flag",      &flag      ); 
00850     APDtrees[i]->SetBranchAddress( "flagAB",    &flagAB    ); 
00851     APDtrees[i]->SetBranchAddress( "pn0",       &pn0       );
00852     APDtrees[i]->SetBranchAddress( "pn1",       &pn1       ); 
00853   
00854   }
00855 
00856 
00857   for (unsigned int iref=0;iref<nRefChan;iref++){
00858     for (unsigned int imod=0;imod<nMod;imod++){
00859       
00860       int jmod=modules[imod];
00861       
00862       std::stringstream nameref;
00863       nameref << "refAPDTree" <<imod<<"_"<<iref;
00864       
00865       RefAPDtrees[iref][jmod]= new TTree(nameref.str().c_str(),nameref.str().c_str());
00866       
00867       RefAPDtrees[iref][jmod]->Branch( "eventref",     &eventref,     "eventref/I"     );
00868       RefAPDtrees[iref][jmod]->Branch( "colorref",     &colorref,     "colorref/I"     );
00869       if(iref==0) RefAPDtrees[iref][jmod]->Branch( "apdAmplA",   &apdAmplA,   "apdAmplA/D"   );
00870       if(iref==1) RefAPDtrees[iref][jmod]->Branch( "apdAmplB",   &apdAmplB,   "apdAmplB/D"   );
00871       
00872       RefAPDtrees[iref][jmod]->SetBranchAddress( "eventref",     &eventref     );
00873       RefAPDtrees[iref][jmod]->SetBranchAddress( "colorref",     &colorref     ); 
00874       if(iref==0)RefAPDtrees[iref][jmod]->SetBranchAddress( "apdAmplA",   &apdAmplA   );
00875       if(iref==1)RefAPDtrees[iref][jmod]->SetBranchAddress( "apdAmplB",   &apdAmplB   );
00876       
00877     } 
00878   }
00879  
00880   assert( colors.size()<= nColor );
00881   unsigned int nCol=colors.size();
00882 
00883   // Declare PN stuff
00884   //===================
00885   
00886   for (unsigned int iM=0;iM<nMod;iM++){
00887     unsigned int iMod=modules[iM]-1;
00888 
00889     for (unsigned int ich=0;ich<nPNPerMod;ich++){
00890       for (unsigned int icol=0;icol<nCol;icol++){
00891         PNFirstAnal[iMod][ich][icol]=new TPN(ich);
00892         PNAnal[iMod][ich][icol]=new TPN(ich);
00893       }
00894 
00895     }
00896   }
00897   
00898 
00899   // Declare function for APD ampl fit
00900   //===================================
00901 
00902   PulseFitWithFunction * pslsfit = new PulseFitWithFunction();
00903   double chi2;
00904 
00905   
00906   for (unsigned int iCry=0;iCry<nCrys;iCry++){
00907     for (unsigned int icol=0;icol<nCol;icol++){
00908 
00909       // Declare APD stuff
00910       //===================
00911       
00912       APDFirstAnal[iCry][icol] = new TAPD();
00913       IsThereDataADC[iCry][icol]=1;      
00914       std::stringstream cut;
00915       cut <<"color=="<<colors[icol];
00916       if(ADCtrees[iCry]->GetEntries(cut.str().c_str())<10) IsThereDataADC[iCry][icol]=0;  
00917       
00918     }
00919 
00920     unsigned int iMod=iModule[iCry]-1;
00921     double alpha, beta;
00922     
00923     // Loop on events
00924     //================
00925     
00926     Long64_t nbytes = 0, nb = 0;
00927     for (Long64_t jentry=0; jentry< ADCtrees[iCry]->GetEntriesFast();jentry++){
00928       nb = ADCtrees[iCry]->GetEntry(jentry);   nbytes += nb;    
00929       
00930       // Get back color
00931  
00932       unsigned int iCol=0;
00933       for(unsigned int i=0;i<nCol;i++){
00934         if(color==colors[i]) {
00935           iCol=i;
00936           i=colors.size();
00937         }
00938       }
00939       
00940       // Retreive alpha and beta
00941       
00942       std::vector<double> abvals=shapana->getVals(iCry);
00943       alpha=abvals[0];
00944       beta=abvals[1];
00945       flagAB=int(abvals[4]);
00946       iphi=iPhi[iCry];
00947       ieta=iEta[iCry];
00948       towerID=iTowerID[iCry];
00949       channelID=iChannelID[iCry];
00950       
00951       // Amplitude calculation
00952 
00953       APDPulse->setPulse(adc);
00954       adcNoPed=APDPulse->getAdcWithoutPedestal();
00955 
00956       apdAmpl=0;
00957       apdAmplA=0;
00958       apdAmplB=0;
00959       apdTime=0;
00960       
00961       if( APDPulse->isPulseOK()){
00962         
00963         pslsfit -> init(_nsamples,_firstsample,_lastsample,_niter,  alpha,  beta);
00964         chi2 = pslsfit -> doFit(&adcNoPed[0]);
00965         
00966         if( chi2 < 0. || chi2 == 102 || chi2 == 101  ) {          
00967           apdAmpl=0;      
00968           apdTime=0;
00969           flag=0;
00970         }else{
00971           apdAmpl = pslsfit -> getAmpl();
00972           apdTime = pslsfit -> getTime();
00973           flag=1;
00974         }
00975       }else {
00976         chi2=0.0;
00977         apdAmpl=0;        
00978         apdTime=0;
00979         flag=0;
00980       }
00981       
00982       if (_debug==1) std::cout <<"-- debug test -- apdAmpl="<<apdAmpl<< 
00983                        ", apdTime="<< apdTime<< std::endl;
00984       double pnmean;
00985       if (pn0<10 && pn1>10) {
00986         pnmean=pn1;
00987       }else if (pn1<10 && pn0>10){
00988         pnmean=pn0;
00989       }else pnmean=0.5*(pn0+pn1);
00990       
00991       if (_debug==1) std::cout <<"-- debug test -- pn0="<<pn0<<", pn1="<<pn1<< std::endl;
00992       
00993       // Fill PN stuff
00994       //===============
00995 
00996       if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){ 
00997         for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
00998           PNFirstAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
00999         }
01000       }
01001       
01002       // Fill APD stuff 
01003       //================
01004 
01005       if( APDPulse->isPulseOK() ){
01006         APDFirstAnal[iCry][iCol]->addEntry(apdAmpl,pnmean, pn0, pn1, apdTime);      
01007         APDtrees[iCry]->Fill();    
01008         
01009         // Fill reference trees
01010         //=====================
01011         
01012         if( apdRefMap[0][iMod+1]==iCry || apdRefMap[1][iMod+1]==iCry ) {
01013           
01014           apdAmplA=0.0;
01015           apdAmplB=0.0;      
01016           eventref=event;
01017           colorref=color;
01018           
01019           for(unsigned int ir=0;ir<nRefChan;ir++){
01020             
01021             if(apdRefMap[ir][iMod+1]==iCry) {
01022               if (ir==0) apdAmplA=apdAmpl;
01023               else if(ir==1) apdAmplB=apdAmpl;
01024               RefAPDtrees[ir][iMod+1]->Fill(); 
01025             }   
01026           }
01027         }
01028       }
01029     }
01030   }
01031   
01032   delete pslsfit;
01033   
01034   ADCFile->Close(); 
01035   
01036   // Remove temporary file
01037   //=======================
01038   std::stringstream del;
01039   del << "rm " <<ADCfile;
01040   system(del.str().c_str());
01041   
01042   
01043   // Create output trees
01044   //=====================
01045 
01046   resFile = new TFile(resfile.c_str(),"RECREATE");
01047   
01048   for (unsigned int iColor=0;iColor<nCol;iColor++){
01049     
01050     std::stringstream nametree;
01051     nametree <<"APDCol"<<colors[iColor];
01052     std::stringstream nametree2;
01053     nametree2 <<"PNCol"<<colors[iColor];
01054     
01055     restrees[iColor]= new TTree(nametree.str().c_str(),nametree.str().c_str());
01056     respntrees[iColor]= new TTree(nametree2.str().c_str(),nametree2.str().c_str());
01057     
01058     restrees[iColor]->Branch( "iphi",        &iphi,        "iphi/I"           );
01059     restrees[iColor]->Branch( "ieta",        &ieta,        "ieta/I"           );
01060     restrees[iColor]->Branch( "side",        &side,        "side/I"           );
01061     restrees[iColor]->Branch( "dccID",       &dccID,       "dccID/I"          );
01062     restrees[iColor]->Branch( "moduleID",    &moduleID,    "moduleID/I"       );
01063     restrees[iColor]->Branch( "towerID",     &towerID,     "towerID/I"        );
01064     restrees[iColor]->Branch( "channelID",   &channelID,   "channelID/I"      );
01065     restrees[iColor]->Branch( "APD",         &APD,         "APD[6]/D"         );
01066     restrees[iColor]->Branch( "Time",        &Time,        "Time[6]/D"        );
01067     restrees[iColor]->Branch( "APDoPN",      &APDoPN,      "APDoPN[6]/D"      );
01068     restrees[iColor]->Branch( "APDoPNA",     &APDoPNA,     "APDoPNA[6]/D"     );
01069     restrees[iColor]->Branch( "APDoPNB",     &APDoPNB,     "APDoPNB[6]/D"     );
01070     restrees[iColor]->Branch( "APDoAPDA",    &APDoAPDA,    "APDoAPDA[6]/D"    );
01071     restrees[iColor]->Branch( "APDoAPDB",    &APDoAPDB,    "APDoAPDB[6]/D"    );
01072     restrees[iColor]->Branch( "flag",        &flag,        "flag/I"           ); 
01073 
01074     respntrees[iColor]->Branch( "side",      &side,      "side/I"         );
01075     respntrees[iColor]->Branch( "moduleID",  &moduleID,  "moduleID/I"     );
01076     respntrees[iColor]->Branch( "pnID",      &pnID,      "pnID/I"         );
01077     respntrees[iColor]->Branch( "PN",        &PN,        "PN[6]/D"        ); 
01078     respntrees[iColor]->Branch( "PNoPN",     &PNoPN,     "PNoPN[6]/D"     );
01079     respntrees[iColor]->Branch( "PNoPNA",    &PNoPNA,    "PNoPNA[6]/D"    );
01080     respntrees[iColor]->Branch( "PNoPNB",    &PNoPNB,    "PNoPNB[6]/D"    );
01081 
01082     restrees[iColor]->SetBranchAddress( "iphi",        &iphi       );
01083     restrees[iColor]->SetBranchAddress( "ieta",        &ieta       );
01084     restrees[iColor]->SetBranchAddress( "side",        &side       );
01085     restrees[iColor]->SetBranchAddress( "dccID",       &dccID      );
01086     restrees[iColor]->SetBranchAddress( "moduleID",    &moduleID   );
01087     restrees[iColor]->SetBranchAddress( "towerID",     &towerID    );
01088     restrees[iColor]->SetBranchAddress( "channelID",   &channelID  );
01089     restrees[iColor]->SetBranchAddress( "APD",         APD         ); 
01090     restrees[iColor]->SetBranchAddress( "Time",        Time        );   
01091     restrees[iColor]->SetBranchAddress( "APDoPN",      APDoPN      ); 
01092     restrees[iColor]->SetBranchAddress( "APDoPNA",     APDoPNA     );
01093     restrees[iColor]->SetBranchAddress( "APDoPNB",     APDoPNB     );
01094     restrees[iColor]->SetBranchAddress( "APDoAPDA",    APDoAPDA    );
01095     restrees[iColor]->SetBranchAddress( "APDoAPDB",    APDoAPDB    );
01096     restrees[iColor]->SetBranchAddress( "flag",        &flag       ); 
01097    
01098     respntrees[iColor]->SetBranchAddress( "side",      &side     );
01099     respntrees[iColor]->SetBranchAddress( "moduleID",  &moduleID );
01100     respntrees[iColor]->SetBranchAddress( "pnID",      &pnID     );
01101     respntrees[iColor]->SetBranchAddress( "PN",        PN        );
01102     respntrees[iColor]->SetBranchAddress( "PNoPN",     PNoPN     );
01103     respntrees[iColor]->SetBranchAddress( "PNoPNA",    PNoPNA    );
01104     respntrees[iColor]->SetBranchAddress( "PNoPNB",    PNoPNB    );
01105   
01106   }
01107 
01108   // Set Cuts for PN stuff
01109   //=======================
01110 
01111   for (unsigned int iM=0;iM<nMod;iM++){
01112     unsigned int iMod=modules[iM]-1;
01113     
01114     for (unsigned int ich=0;ich<nPNPerMod;ich++){
01115       for (unsigned int icol=0;icol<nCol;icol++){
01116         PNAnal[iMod][ich][icol]->setPNCut(PNFirstAnal[iMod][ich][icol]->getPN().at(0),
01117                                           PNFirstAnal[iMod][ich][icol]->getPN().at(1));
01118       }
01119     }
01120   }
01121   
01122   // Build ref trees indexes
01123   //========================
01124   for(unsigned int imod=0;imod<nMod;imod++){
01125     int jmod=modules[imod];
01126     if( RefAPDtrees[0][jmod]->GetEntries()!=0 && RefAPDtrees[1][jmod]->GetEntries()!=0 ){
01127       RefAPDtrees[0][jmod]->BuildIndex("eventref");
01128       RefAPDtrees[1][jmod]->BuildIndex("eventref");
01129     }
01130   }
01131   
01132   
01133   // Final loop on crystals
01134   //=======================
01135 
01136   for (unsigned int iCry=0;iCry<nCrys;iCry++){ 
01137     
01138     unsigned int iMod=iModule[iCry]-1;
01139     
01140     // Set cuts on APD stuff
01141     //=======================
01142     
01143     for(unsigned int iCol=0;iCol<nCol;iCol++){   
01144       
01145       std::vector<double> lowcut;
01146       std::vector<double> highcut;
01147       double cutMin;
01148       double cutMax;
01149 
01150       cutMin=APDFirstAnal[iCry][iCol]->getAPD().at(0)-2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
01151       if(cutMin<0) cutMin=0;
01152       cutMax=APDFirstAnal[iCry][iCol]->getAPD().at(0)+2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
01153       
01154       lowcut.push_back(cutMin);
01155       highcut.push_back(cutMax);
01156       
01157       cutMin=APDFirstAnal[iCry][iCol]->getTime().at(0)-2.0*APDFirstAnal[iCry][iCol]->getTime().at(1);
01158       cutMax=APDFirstAnal[iCry][iCol]->getTime().at(0)+2.0*APDFirstAnal[iCry][iCol]->getTime().at(1); 
01159       lowcut.push_back(cutMin);
01160       highcut.push_back(cutMax);
01161       
01162       APDAnal[iCry][iCol]=new TAPD();
01163       APDAnal[iCry][iCol]->setAPDCut(APDFirstAnal[iCry][iCol]->getAPD().at(0),APDFirstAnal[iCry][iCol]->getAPD().at(1));
01164       APDAnal[iCry][iCol]->setAPDoPNCut(APDFirstAnal[iCry][iCol]->getAPDoPN().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN().at(1));
01165       APDAnal[iCry][iCol]->setAPDoPN0Cut(APDFirstAnal[iCry][iCol]->getAPDoPN0().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN0().at(1));
01166       APDAnal[iCry][iCol]->setAPDoPN1Cut(APDFirstAnal[iCry][iCol]->getAPDoPN1().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN1().at(1));
01167       APDAnal[iCry][iCol]->setTimeCut(APDFirstAnal[iCry][iCol]->getTime().at(0),APDFirstAnal[iCry][iCol]->getTime().at(1));
01168       APDAnal[iCry][iCol]->set2DAPDoAPD0Cut(lowcut,highcut);
01169       APDAnal[iCry][iCol]->set2DAPDoAPD1Cut(lowcut,highcut);
01170       
01171     }
01172 
01173     // Final loop on events
01174     //=======================
01175 
01176     Long64_t nbytes = 0, nb = 0;
01177     for (Long64_t jentry=0; jentry< APDtrees[iCry]->GetEntriesFast();jentry++) { 
01178       nb = APDtrees[iCry]->GetEntry(jentry);   nbytes += nb; 
01179       
01180       double pnmean;
01181       if (pn0<10 && pn1>10) {
01182         pnmean=pn1;
01183       }else if (pn1<10 && pn0>10){
01184         pnmean=pn0;
01185       }else pnmean=0.5*(pn0+pn1);
01186       
01187       // Get back color
01188       //================
01189       
01190       unsigned int iCol=0;
01191       for(unsigned int i=0;i<nCol;i++){
01192         if(color==colors[i]) {
01193           iCol=i;
01194           i=colors.size();
01195         }
01196       }
01197       
01198       // Fill PN stuff
01199       //===============
01200       
01201       if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){ 
01202         for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
01203           PNAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
01204         }
01205       }
01206       
01207       // Get ref amplitudes
01208       //===================
01209 
01210       if (_debug==1) std::cout <<"-- debug test -- Last Loop event:"<<event<<" apdAmpl:"<< apdAmpl<< std::endl;
01211       apdAmplA = 0.0;
01212       apdAmplB = 0.0;
01213       
01214       for (unsigned int iRef=0;iRef<nRefChan;iRef++){
01215         RefAPDtrees[iRef][iMod+1]->GetEntryWithIndex(event); 
01216       }
01217       
01218       if (_debug==1 ) std::cout <<"-- debug test -- Last Loop apdAmplA:"<<apdAmplA<< " apdAmplB:"<< apdAmplB<<", event:"<< event<<", eventref:"<< eventref<< std::endl;
01219       
01220       
01221       // Fill APD stuff
01222       //===============
01223       
01224       APDAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime, apdAmplA, apdAmplB); 
01225    
01226     }
01227     
01228     moduleID=iMod+1;
01229     
01230     if( moduleID>=20 ) moduleID-=2; // Trick to fix endcap specificity   
01231     
01232     // Get final results for APD
01233     //===========================
01234     
01235     for(unsigned int iColor=0;iColor<nCol;iColor++){
01236       
01237       std::vector<double> apdvec = APDAnal[iCry][iColor]->getAPD();
01238       std::vector<double> apdpnvec = APDAnal[iCry][iColor]->getAPDoPN();
01239       std::vector<double> apdpn0vec = APDAnal[iCry][iColor]->getAPDoPN0();
01240       std::vector<double> apdpn1vec = APDAnal[iCry][iColor]->getAPDoPN1();
01241       std::vector<double> timevec = APDAnal[iCry][iColor]->getTime();
01242       std::vector<double> apdapd0vec = APDAnal[iCry][iColor]->getAPDoAPD0();
01243       std::vector<double> apdapd1vec = APDAnal[iCry][iColor]->getAPDoAPD1();
01244       
01245       
01246       for(unsigned int i=0;i<apdvec.size();i++){
01247         
01248         APD[i]=apdvec.at(i);
01249         APDoPN[i]=apdpnvec.at(i);
01250         APDoPNA[i]=apdpn0vec.at(i);
01251         APDoPNB[i]=apdpn1vec.at(i);
01252         APDoAPDA[i]=apdapd0vec.at(i);
01253         APDoAPDB[i]=apdapd1vec.at(i);
01254         Time[i]=timevec.at(i);
01255       }
01256       
01257       
01258       // Fill APD results trees
01259       //========================
01260 
01261       iphi=iPhi[iCry];
01262       ieta=iEta[iCry];
01263       dccID=idccID[iCry];
01264       side=iside[iCry];
01265       towerID=iTowerID[iCry];
01266       channelID=iChannelID[iCry];
01267 
01268       
01269       if( !wasGainOK[iCry] || !wasTimingOK[iCry] || IsThereDataADC[iCry][iColor]==0 ){
01270         flag=0;
01271       }else flag=1;
01272       
01273       restrees[iColor]->Fill();   
01274     }
01275   }
01276    
01277   // Get final results for PN 
01278   //==========================
01279   
01280   for (unsigned int iM=0;iM<nMod;iM++){
01281     unsigned int iMod=modules[iM]-1;
01282 
01283     side=iside[firstChanMod[iMod]];
01284     
01285     for (unsigned int ch=0;ch<nPNPerMod;ch++){
01286 
01287       pnID=ch;
01288       moduleID=iMod+1;
01289       
01290       if( moduleID>=20 ) moduleID-=2;  // Trick to fix endcap specificity
01291 
01292       for(unsigned int iColor=0;iColor<nCol;iColor++){
01293 
01294         std::vector<double> pnvec = PNAnal[iMod][ch][iColor]->getPN();
01295         std::vector<double> pnopnvec = PNAnal[iMod][ch][iColor]->getPNoPN();
01296         std::vector<double> pnopn0vec = PNAnal[iMod][ch][iColor]->getPNoPN0();
01297         std::vector<double> pnopn1vec = PNAnal[iMod][ch][iColor]->getPNoPN1();
01298         
01299         for(unsigned int i=0;i<pnvec.size();i++){
01300           
01301           PN[i]=pnvec.at(i);
01302           PNoPN[i]=pnopnvec.at(i);
01303           PNoPNA[i]=pnopn0vec.at(i);
01304           PNoPNB[i]=pnopn1vec.at(i);
01305         }
01306         
01307         // Fill PN results trees
01308         //========================
01309 
01310         respntrees[iColor]->Fill();
01311       }
01312     }
01313   }
01314   
01315   // Remove temporary files
01316   //========================
01317   if(!_saveallevents){
01318 
01319     APDFile->Close();
01320     std::stringstream del2;
01321     del2 << "rm " <<APDfile;
01322     system(del2.str().c_str());
01323 
01324   }else {
01325     
01326     APDFile->cd();
01327     APDtrees[0]->Write();
01328     
01329     APDFile->Close();
01330     resFile->cd();
01331   }
01332   
01333   // Save results 
01334   //===============
01335 
01336   for (unsigned int i=0;i<nCol;i++){
01337     restrees[i]->Write();
01338     respntrees[i]->Write();
01339   }  
01340   
01341   resFile->Close(); 
01342   
01343   std::cout <<    "\t+=+    .................................................. done  +=+" << std::endl;
01344   std::cout <<    "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
01345 }
01346 
01347 
01348 void EcalLaserAnalyzer::setGeomEB(int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr){
01349   
01350   side=MEEBGeom::side(etaG,phiG);
01351   
01352   assert( module>=*min_element(modules.begin(),modules.end()) && module<=*max_element(modules.begin(),modules.end()) );
01353       
01354   eta = etaG;
01355   phi = phiG;
01356   channelID=5*(strip-1) + xtal-1; 
01357   towerID=tower;
01358   
01359   std::vector<int> apdRefChan=ME::apdRefChannels(module, lmr);      
01360   for (unsigned int iref=0;iref<nRefChan;iref++){
01361     if(channelID==apdRefChan[iref] && towerID==apdRefTT 
01362        && apdRefMap[iref].count(module)==0){
01363       apdRefMap[iref][module]=channel;
01364     }
01365   }
01366   
01367   if(isFirstChanModFilled[module-1]==0) {
01368     firstChanMod[module-1]=channel;
01369     isFirstChanModFilled[module-1]=1;
01370   } 
01371   
01372   iEta[channel]=eta;
01373   iPhi[channel]=phi;
01374   iModule[channel]= module ;
01375   iTowerID[channel]=towerID;
01376   iChannelID[channel]=channelID;
01377   idccID[channel]=dccID;
01378   iside[channel]=side;
01379 
01380 }
01381 
01382 
01383 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){
01384   
01385   side=MEEEGeom::side(iX, iY, iZ);
01386 
01387   assert( module>=*min_element(modules.begin(),modules.end()) 
01388           && module<=*max_element(modules.begin(),modules.end()) );
01389   
01390   eta = etaG;
01391   phi = phiG;
01392   channelID=ch;
01393   towerID=tower;
01394   
01395   std::vector<int> apdRefChan=ME::apdRefChannels(module, lmr);      
01396   for (unsigned int iref=0;iref<nRefChan;iref++){
01397     if(channelID==apdRefChan[iref] && towerID==apdRefTT 
01398        && apdRefMap[iref].count(module)==0){
01399       apdRefMap[iref][module]=channel;
01400     }
01401   }
01402   
01403   if(isFirstChanModFilled[module-1]==0) {
01404     firstChanMod[module-1]=channel;
01405     isFirstChanModFilled[module-1]=1;
01406   } 
01407   
01408   iEta[channel]=eta;
01409   iPhi[channel]=phi;
01410   iModule[channel]= module ;
01411   iTowerID[channel]=towerID;
01412   iChannelID[channel]=channelID;
01413   idccID[channel]=dccID;
01414   iside[channel]=side;
01415 
01416 }
01417 
01418 
01419 DEFINE_FWK_MODULE(EcalLaserAnalyzer);
01420
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalLaserAnalyzer.cc
