Reference Manual

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