Reference Manual

Go to the documentation of this file.
00001 /* 
00002  *  \class EcalLaserAnalyzer2
00003  *  $Date: 2010/01/18 17:28:44 $
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   double chi2;
00844 
00845   for (unsigned int iCry=0;iCry<nCrys;iCry++){    
00846     for (unsigned int icol=0;icol<nCol;icol++){
00847       
00848       // Declare APD stuff
00849       //===================
00850       
00851       APDFirstAnal[iCry][icol]=new TAPD();
00852       IsThereDataADC[iCry][icol]=1;      
00853       stringstream cut;
00854       cut <<"color=="<<colors.at(icol);
00855       if(ADCtrees[iCry]->GetEntries(cut.str().c_str())<10) IsThereDataADC[iCry][icol]=0;
00856       
00857     }
00858 
00859     unsigned int iMod=iModule[iCry]-1;
00860     assert(iMod<=nMod);
00861 
00862     if(isSPRFine) psfit -> init(_nsamples,_firstsample,_lastsample,_niter, nSamplesShapes, shapesVec, _noise );
00863     
00864     // Loop on events
00865     //================
00866 
00867     Long64_t nbytes = 0, nb = 0;
00868     for (Long64_t jentry=0; jentry< ADCtrees[iCry]->GetEntriesFast();jentry++) {  // Loop on events
00869       nb = ADCtrees[iCry]->GetEntry(jentry);   nbytes += nb;    
00870       
00871       flagfit=1;
00872       apdAmpl=0.0;
00873       apdTime=0.0;
00874       ieta=eta; 
00875       iphi=phi;
00876       
00877       // Get back color
00878  
00879       unsigned int iCol=0;
00880       for(unsigned int i=0;i<nCol;i++){
00881         if(color==colors[i]) {
00882           iCol=i;
00883           i=colors.size();
00884         }
00885       }
00886       
00887       // Amplitude calculation
00888 
00889       APDPulse->setPulse(adc);
00890       adcNoPed=APDPulse->getAdcWithoutPedestal();
00891 
00892         
00893       if(isSPRFine && APDPulse->isPulseOK()) {
00894         
00895         chi2 = psfit -> doFit(&adcNoPed[0]);
00896         apdAmpl = psfit -> getAmpl();
00897         apdTime = psfit -> getTime(); 
00898         
00899       }else{
00900         
00901         chi2=0.0;
00902         apdAmpl=0;        
00903         apdTime=0;
00904         flagfit=0;
00905         
00906       }
00907       
00908       if (_debug>=1) cout <<"-- debug test -- endJob -- apdAmpl:"<<apdAmpl
00909                           <<" apdTime:"<< apdTime<< endl; 
00910 
00911       double pnmean;
00912       if (pn0<10 && pn1>10) {
00913         pnmean=pn1;
00914       }else if (pn1<10 && pn0>10){
00915         pnmean=pn0;
00916       }else pnmean=0.5*(pn0+pn1);
00917 
00918       if (_debug>=1) cout <<"-- debug test -- endJob -- pnMean:"<<pnmean << endl; 
00919 
00920       // Fill PN stuff
00921       //===============
00922 
00923       if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){ 
00924         for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
00925           PNFirstAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
00926         }
00927       }
00928       
00929       // Fill APD stuff 
00930       //================
00931  
00932       if(apdAmpl!=0.0) APDFirstAnal[iCry][iCol]->addEntry(apdAmpl,pnmean, pn0, pn1, apdTime);      
00933       if (_debug>=1) cout <<"-- debug test -- endJob -- filling APDTree"<< endl; 
00934 
00935       APDtrees[iCry]->Fill();
00936 
00937         // Fill reference trees
00938         //=====================
00939       
00940       if( apdRefMap[0][iMod+1]==iCry || apdRefMap[1][iMod+1]==iCry) {
00941         
00942         apdAmplA=0.0;
00943         apdAmplB=0.0;      
00944         eventref=event;
00945         colorref=color;
00946         
00947         for(unsigned int ir=0;ir<nRefChan;ir++){
00948           
00949           if (_debug>=1) cout <<"-- debug test -- ir:" << ir <<" tt:"<< towerID<<" refmap:"<<apdRefMap[ir][iMod+1]<< " iCry:"<<iCry<<endl;
00950           
00951           if(apdRefMap[ir][iMod+1]==iCry) {
00952             if (_debug>=1) cout <<"-- debug test -- cut passed " <<endl;
00953             if (ir==0) apdAmplA=apdAmpl;
00954             else if(ir==1) apdAmplB=apdAmpl;
00955             if (_debug>=1) cout <<"-- debug test -- apdAmplA=" <<apdAmplA<<endl;
00956             if (_debug>=1) cout <<"-- debug test -- apdAmplB=" <<apdAmplB<<endl;
00957             if (_debug>=1) cout <<"-- debug test -- color=" <<color<<", event:"<< event<<", ir:" << ir <<" tt-1:"<< towerID-1<< endl;
00958             
00959             RefAPDtrees[ir][iMod+1]->Fill(); 
00960             
00961             if (_debug>=1) cout <<"-- debug test -- tree filled"<< event<<endl;
00962           }     
00963         }
00964       }
00965     }
00966   }
00967 
00968   delete psfit;
00969   
00970   ADCFile->Close();
00971   
00972   if (_debug==1) cout <<"-- debug test -- endJob -- after apdAmpl Loop"<< endl; 
00973 
00974   // Remove temporary file
00975   //=======================
00976 
00977   stringstream del;
00978   del << "rm " <<ADCfile;
00979   system(del.str().c_str()); 
00980 
00981   
00982   // Create output trees
00983   //=====================
00984   
00985   resFile = new TFile(resfile.c_str(),"RECREATE");
00986   
00987   for (unsigned int iColor=0;iColor<nCol;iColor++){
00988 
00989     stringstream nametree;
00990     nametree <<"APDCol"<<colors.at(iColor);
00991     stringstream nametree2;
00992     nametree2 <<"PNCol"<<colors.at(iColor);
00993     
00994     restrees[iColor]= new TTree(nametree.str().c_str(),nametree.str().c_str());
00995     respntrees[iColor]= new TTree(nametree2.str().c_str(),nametree2.str().c_str());
00996 
00997     restrees[iColor]->Branch( "iphi",        &iphi,        "iphi/I"           );
00998     restrees[iColor]->Branch( "ieta",        &ieta,        "ieta/I"           );
00999     restrees[iColor]->Branch( "side",        &side,        "side/I"           );
01000     restrees[iColor]->Branch( "dccID",       &dccID,       "dccID/I"          );
01001     restrees[iColor]->Branch( "moduleID",    &moduleID,    "moduleID/I"       );
01002     restrees[iColor]->Branch( "towerID",     &towerID,     "towerID/I"        );
01003     restrees[iColor]->Branch( "channelID",   &channelID,   "channelID/I"      );
01004     restrees[iColor]->Branch( "APD",         &APD,         "APD[6]/D"         );
01005     restrees[iColor]->Branch( "Time",        &Time,        "Time[6]/D"        );
01006     restrees[iColor]->Branch( "APDoPN",      &APDoPN,      "APDoPN[6]/D"      );
01007     restrees[iColor]->Branch( "APDoPNA",     &APDoPNA,     "APDoPNA[6]/D"     );
01008     restrees[iColor]->Branch( "APDoPNB",     &APDoPNB,     "APDoPNB[6]/D"     );
01009     restrees[iColor]->Branch( "APDoAPDA",    &APDoAPDA,    "APDoAPDA[6]/D"    );
01010     restrees[iColor]->Branch( "APDoAPDB",    &APDoAPDB,    "APDoAPDB[6]/D"    );
01011     restrees[iColor]->Branch( "ShapeCor",    &ShapeCor,    "ShapeCor/D"       );
01012     restrees[iColor]->Branch( "flag",        &flag,        "flag/I"           ); 
01013 
01014 
01015     respntrees[iColor]->Branch( "moduleID",  &moduleID,  "moduleID/I"     );
01016     respntrees[iColor]->Branch( "pnID",      &pnID,      "pnID/I"         );
01017     respntrees[iColor]->Branch( "PN",        &PN,        "PN[6]/D"        ); 
01018     respntrees[iColor]->Branch( "PNoPN",     &PNoPN,     "PNoPN[6]/D"     );
01019     respntrees[iColor]->Branch( "PNoPNA",    &PNoPNA,    "PNoPNA[6]/D"    );
01020     respntrees[iColor]->Branch( "PNoPNB",    &PNoPNB,    "PNoPNB[6]/D"    );
01021 
01022     restrees[iColor]->SetBranchAddress( "iphi",        &iphi       );
01023     restrees[iColor]->SetBranchAddress( "ieta",        &ieta       );
01024     restrees[iColor]->SetBranchAddress( "dccID",       &dccID      );
01025     restrees[iColor]->SetBranchAddress( "moduleID",    &moduleID   );
01026     restrees[iColor]->SetBranchAddress( "towerID",     &towerID    );
01027     restrees[iColor]->SetBranchAddress( "channelID",   &channelID  );
01028     restrees[iColor]->SetBranchAddress( "APD",         APD         ); 
01029     restrees[iColor]->SetBranchAddress( "Time",        Time        );   
01030     restrees[iColor]->SetBranchAddress( "APDoPN",      APDoPN      ); 
01031     restrees[iColor]->SetBranchAddress( "APDoPNA",     APDoPNA     );
01032     restrees[iColor]->SetBranchAddress( "APDoPNB",     APDoPNB     );
01033     restrees[iColor]->SetBranchAddress( "APDoAPDA",    APDoAPDA    );
01034     restrees[iColor]->SetBranchAddress( "APDoAPDB",    APDoAPDB    );
01035     restrees[iColor]->SetBranchAddress( "ShapeCor",    &ShapeCor   );
01036     restrees[iColor]->SetBranchAddress( "flag",        &flag       ); 
01037    
01038     respntrees[iColor]->SetBranchAddress( "moduleID",  &moduleID );
01039     respntrees[iColor]->SetBranchAddress( "pnID",      &pnID     );
01040     respntrees[iColor]->SetBranchAddress( "PN",        PN        );
01041     respntrees[iColor]->SetBranchAddress( "PNoPN",     PNoPN     );
01042     respntrees[iColor]->SetBranchAddress( "PNoPNA",    PNoPNA    );
01043     respntrees[iColor]->SetBranchAddress( "PNoPNB",    PNoPNB    );
01044   
01045   }
01046   
01047 
01048   // Set Cuts for PN stuff
01049   //=======================
01050 
01051   for (unsigned int iM=0;iM<nMod;iM++){
01052     unsigned int iMod=modules[iM]-1;
01053     
01054     for (unsigned int ich=0;ich<nPNPerMod;ich++){
01055       for (unsigned int icol=0;icol<nCol;icol++){
01056         PNAnal[iMod][ich][icol]->setPNCut(PNFirstAnal[iMod][ich][icol]->getPN().at(0),PNFirstAnal[iMod][ich][icol]->getPN().at(1));
01057       }
01058     }
01059   }
01060 
01061   // Build ref trees indexes
01062   //========================
01063   
01064   for(unsigned int imod=0;imod<nMod;imod++){
01065     int jmod=modules[imod];
01066     if( RefAPDtrees[0][jmod]->GetEntries()!=0 && RefAPDtrees[1][jmod]->GetEntries()!=0 ){
01067       RefAPDtrees[0][jmod]->BuildIndex("eventref");
01068       RefAPDtrees[1][jmod]->BuildIndex("eventref");
01069     }
01070   }
01071   
01072   
01073   // Final loop on crystals
01074   //=======================
01075   
01076   for (unsigned int iCry=0;iCry<nCrys;iCry++){ 
01077     
01078     unsigned int iMod=iModule[iCry]-1;
01079     
01080     // Set cuts on APD stuff
01081     //=======================
01082 
01083     for(unsigned int iCol=0;iCol<nCol;iCol++){ 
01084       
01085       std::vector<double> lowcut;
01086       std::vector<double> highcut;
01087       double cutMin;
01088       double cutMax;
01089       
01090       cutMin=APDFirstAnal[iCry][iCol]->getAPD().at(0)-2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
01091       if(cutMin<0) cutMin=0;
01092       cutMax=APDFirstAnal[iCry][iCol]->getAPD().at(0)+2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
01093       
01094       lowcut.push_back(cutMin);
01095       highcut.push_back(cutMax);
01096       
01097       cutMin=APDFirstAnal[iCry][iCol]->getTime().at(0)-2.0*APDFirstAnal[iCry][iCol]->getTime().at(1);
01098       cutMax=APDFirstAnal[iCry][iCol]->getTime().at(0)+2.0*APDFirstAnal[iCry][iCol]->getTime().at(1); 
01099       lowcut.push_back(cutMin);
01100       highcut.push_back(cutMax);
01101       
01102       
01103       APDAnal[iCry][iCol]=new TAPD();
01104       APDAnal[iCry][iCol]->setAPDCut(APDFirstAnal[iCry][iCol]->getAPD().at(0),APDFirstAnal[iCry][iCol]->getAPD().at(1));
01105       APDAnal[iCry][iCol]->setAPDoPNCut(APDFirstAnal[iCry][iCol]->getAPDoPN().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN().at(1));
01106       APDAnal[iCry][iCol]->setAPDoPN0Cut(APDFirstAnal[iCry][iCol]->getAPDoPN0().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN0().at(1));
01107       APDAnal[iCry][iCol]->setAPDoPN1Cut(APDFirstAnal[iCry][iCol]->getAPDoPN1().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN1().at(1));
01108       APDAnal[iCry][iCol]->setTimeCut(APDFirstAnal[iCry][iCol]->getTime().at(0),APDFirstAnal[iCry][iCol]->getTime().at(1));
01109       
01110       
01111       APDAnal[iCry][iCol]->set2DAPDoAPD0Cut(lowcut,highcut);
01112       APDAnal[iCry][iCol]->set2DAPDoAPD1Cut(lowcut,highcut);
01113       
01114     }
01115     
01116     // Final loop on events
01117     //=======================
01118 
01119     Long64_t nbytes = 0, nb = 0;
01120     for (Long64_t jentry=0; jentry< APDtrees[iCry]->GetEntriesFast();jentry++) { 
01121       nb = APDtrees[iCry]->GetEntry(jentry);   nbytes += nb; 
01122 
01123       double pnmean;
01124       if (pn0<10 && pn1>10) {
01125         pnmean=pn1;
01126       }else if (pn1<10 && pn0>10){
01127         pnmean=pn0;
01128       }else pnmean=0.5*(pn0+pn1);
01129 
01130       // Get back color
01131       //===============
01132  
01133       unsigned int iCol=0;
01134       for(unsigned int i=0;i<nCol;i++){
01135         if(color==colors[i]) {
01136           iCol=i;
01137           i=colors.size();
01138         }
01139       }
01140 
01141       // Fill PN stuff
01142       //===============
01143       
01144       if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){ 
01145         for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
01146           PNAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
01147         }
01148       }
01149       
01150       // Get ref amplitudes
01151       //===================
01152       
01153       if (_debug>=1) cout <<"-- debug test -- LastLoop event:"<<event<<" apdAmpl:"<< apdAmpl<< endl;
01154       apdAmplA = 0.0;
01155       apdAmplB = 0.0;
01156       
01157       for (unsigned int iRef=0;iRef<nRefChan;iRef++){
01158         RefAPDtrees[iRef][iMod+1]->GetEntryWithIndex(event); 
01159       }
01160       
01161       if (_debug==1 ) cout <<"-- debug test -- LastLoop apdAmplA:"<<apdAmplA<< " apdAmplB:"<< apdAmplB<<", event:"<< event<<", eventref:"<< eventref<< endl;
01162 
01163       // Fill APD stuff
01164       //===============
01165 
01166       APDAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime, apdAmplA, apdAmplB);
01167     }
01168 
01169     moduleID=iMod+1;
01170     
01171     if( moduleID>=20 ) moduleID-=2;  // Trick to fix endcap specificity
01172 
01173     // Get final results for APD 
01174     //===========================
01175     
01176     for(unsigned int iColor=0;iColor<nCol;iColor++){
01177       
01178 
01179       std::vector<double> apdvec = APDAnal[iCry][iColor]->getAPD();
01180       std::vector<double> apdpnvec = APDAnal[iCry][iColor]->getAPDoPN();
01181       std::vector<double> apdpn0vec = APDAnal[iCry][iColor]->getAPDoPN0();
01182       std::vector<double> apdpn1vec = APDAnal[iCry][iColor]->getAPDoPN1();
01183       std::vector<double> timevec = APDAnal[iCry][iColor]->getTime();
01184       std::vector<double> apdapd0vec = APDAnal[iCry][iColor]->getAPDoAPD0();
01185       std::vector<double> apdapd1vec = APDAnal[iCry][iColor]->getAPDoAPD1();
01186       
01187       
01188       for(unsigned int i=0;i<apdvec.size();i++){
01189         
01190         APD[i]=apdvec.at(i);
01191         APDoPN[i]=apdpnvec.at(i);
01192         APDoPNA[i]=apdpn0vec.at(i);
01193         APDoPNB[i]=apdpn1vec.at(i);
01194         APDoAPDA[i]=apdapd0vec.at(i);
01195         APDoAPDB[i]=apdapd1vec.at(i);
01196         Time[i]=timevec.at(i);
01197         ShapeCor=shapeCorrection;
01198 
01199       }
01200       
01201       // Fill APD results trees
01202       //========================
01203 
01204       iphi=iPhi[iCry];
01205       ieta=iEta[iCry];
01206       dccID=idccID[iCry];
01207       towerID=iTowerID[iCry];
01208       side=iside[iCry];
01209       channelID=iChannelID[iCry];
01210       
01211       if( !wasGainOK[iCry] || !wasTimingOK[iCry] || IsThereDataADC[iCry][iColor]==0 ) flag=1;
01212       else flag=0;
01213       
01214       if (_debug>=1) cout <<"-- debug test -- endJob -- APD[0]"<< APD[0]<<" APDoPN[0] "<<APDoPN[0]<<" APDoAPDA[0] "<<APDoAPDA[0]<< " flag "<< flag<< endl; 
01215       restrees[iColor]->Fill(); 
01216       
01217     }
01218    
01219   }
01220   
01221   // Get final results for PN 
01222   //==========================
01223 
01224   for (unsigned int iM=0;iM<nMod;iM++){
01225     unsigned int iMod=modules[iM]-1;
01226     
01227     side=iside[firstChanMod[iMod]];
01228 
01229     for (unsigned int ch=0;ch<nPNPerMod;ch++){
01230       
01231       pnID=ch;
01232       moduleID=iMod+1;
01233 
01234       if( moduleID>=20 ) moduleID-=2;  // Trick to fix endcap specificity
01235 
01236       for(unsigned int iColor=0;iColor<nCol;iColor++){
01237 
01238         std::vector<double> pnvec = PNAnal[iMod][ch][iColor]->getPN();
01239         std::vector<double> pnopnvec = PNAnal[iMod][ch][iColor]->getPNoPN();
01240         std::vector<double> pnopn0vec = PNAnal[iMod][ch][iColor]->getPNoPN0();
01241         std::vector<double> pnopn1vec = PNAnal[iMod][ch][iColor]->getPNoPN1();
01242         
01243         for(unsigned int i=0;i<pnvec.size();i++){
01244           
01245           PN[i]=pnvec.at(i);
01246           PNoPN[i]=pnopnvec.at(i);
01247           PNoPNA[i]=pnopn0vec.at(i);
01248           PNoPNB[i]=pnopn1vec.at(i);
01249         }
01250         
01251         if (_debug>=1) cout <<"-- debug test -- endJob -- filling pn results'tree: PN[0]:"<<PN[0]<<" iModule:" << iMod<<" iColor:"<<iColor<<" ch:"<< ch<< endl; 
01252         
01253         // Fill PN results trees
01254         //========================
01255         
01256         respntrees[iColor]->Fill();
01257       }
01258     }
01259   }
01260 
01261   // Remove temporary files
01262   //========================
01263  
01264   if(!_saveallevents){
01265     
01266     APDFile->Close(); 
01267     stringstream del2;
01268     del2 << "rm " <<APDfile;
01269     system(del2.str().c_str());
01270 
01271   }else {
01272     
01273     APDFile->cd();
01274     APDtrees[0]->Write();
01275     APDFile->Close();  
01276     resFile->cd();  
01277   }
01278 
01279 
01280   // Save results 
01281   //===============
01282 
01283   for (unsigned int i=0;i<nCol;i++){
01284     restrees[i]->Write();
01285     respntrees[i]->Write();
01286   }
01287   
01288   resFile->Close(); 
01289     
01290   cout <<    "\t+=+    .................................................. done  +=+" << endl;
01291   cout <<    "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
01292 }
01293 
01294 
01295 
01296 //========================================================================
01297 bool EcalLaserAnalyzer2::getShapes() {
01298 //========================================================================
01299   
01300 
01301   // Get Pulse From Matacq Analysis:
01302   //================================
01303 
01304   bool IsMatacqOK=false;
01305 
01306   int doesMatFileExist=0;
01307   int doesMatShapeExist=0;
01308   FILE *test2;   
01309   TProfile *laserShape=0;
01310   test2 = fopen(matfile.c_str(),"r");
01311   if (test2) doesMatFileExist=1; 
01312   
01313   TFile *MatShapeFile;
01314   if (doesMatFileExist==1){
01315     MatShapeFile = new TFile(matfile.c_str());
01316     laserShape= (TProfile*) MatShapeFile->Get("shapeLaser");
01317     if(laserShape){
01318       doesMatShapeExist=1;
01319       double y=laserShape->Integral("w");
01320       if(y!=0)laserShape->Scale(1.0/y);
01321     }
01322   }else{
01323 
01324     cout <<" ERROR! Matacq shape file not found !"<< endl;
01325      
01326   }
01327   if (doesMatShapeExist) IsMatacqOK=true;
01328 
01329   // Get SPR from the average elec shape in SM6:
01330   //============================================
01331 
01332   int doesElecFileExist=0;
01333   FILE *test; 
01334   test = fopen(elecfile_.c_str(),"r");
01335   if (test) doesElecFileExist=1; 
01336   
01337   TFile *ElecShapesFile;
01338   TH1D* elecShape=0 ;
01339 
01340   if (doesElecFileExist==1){
01341     ElecShapesFile = new TFile(elecfile_.c_str());
01342     stringstream name;
01343     name << "MeanElecShape";
01344     elecShape=(TH1D*) ElecShapesFile->Get(name.str().c_str());
01345     if(elecShape && doesMatShapeExist==1){
01346       double x=elecShape->GetMaximum();
01347       if (x!=0) elecShape->Scale(1.0/x);
01348       isSPRFine=true;
01349     }else{
01350       isSPRFine=false;
01351     }
01352     
01353   }else{
01354     
01355     cout <<" ERROR! Elec shape file not found !"<< endl;
01356     
01357   }
01358   
01359 
01360   if(IsMatacqOK){
01361     
01362     ShapeFile = new TFile(shapefile.c_str(),"RECREATE");
01363     
01364     unsigned int nBins=int(laserShape->GetEntries());
01365     assert( nSamplesShapes == nBins);
01366     double elec_jj, laser_iiMinusjj;
01367     double sum_jj;
01368     
01369     if(isSPRFine==true){
01370       
01371       unsigned int nBins2=int(elecShape->GetNbinsX());
01372       
01373       if(nBins2<nBins){  
01374         cout<< "EcalLaserAnalyzer2::getShapes: wrong configuration of the shapes' number of bins"<< std::endl;
01375         isSPRFine=false;
01376       }
01377       assert( nSamplesShapes == nBins2 );
01378       
01379       stringstream name;
01380       name << "PulseShape";
01381       
01382       PulseShape=new TProfile(name.str().c_str(),name.str().c_str(),nBins,-0.5,double(nBins)-0.5);
01383       
01384       // shift shapes to have max close to the real APD max
01385       
01386       for(int ii=0;ii<50;ii++){
01387         shapes[ii]=0.0;
01388         PulseShape->Fill(ii,0.0);
01389       }
01390       
01391       for(unsigned int ii=0;ii<nBins-50;ii++){
01392         sum_jj=0.0;
01393         for(unsigned int jj=0;jj<ii;jj++){
01394           elec_jj=elecShape->GetBinContent(jj+1);
01395           laser_iiMinusjj=laserShape->GetBinContent(ii-jj+1);
01396           sum_jj+=elec_jj*laser_iiMinusjj;
01397         }
01398         PulseShape->Fill(ii+50,sum_jj);
01399         shapes[ii+50]=sum_jj;
01400       }
01401       
01402       double scale= PulseShape->GetMaximum();
01403       shapeCorrection=scale;
01404       
01405       if(scale!=0) {
01406         PulseShape->Scale(1.0/scale);
01407         for(unsigned int ii=0;ii<nBins;ii++){
01408           shapesVec.push_back(shapes[ii]/scale);
01409         }
01410       }
01411       
01412       if(_saveshapes) PulseShape->Write();
01413     }
01414   }
01415   ShapeFile->Close();
01416 
01417   if(!_saveshapes) {
01418 
01419     stringstream del;
01420     del << "rm " <<shapefile;
01421     system(del.str().c_str()); 
01422     
01423   }
01424   
01425   return IsMatacqOK;
01426 }
01427 
01428 
01429 void EcalLaserAnalyzer2::setGeomEB(int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr){
01430   
01431   side=MEEBGeom::side(etaG,phiG);
01432   
01433   assert( module>=*min_element(modules.begin(),modules.end()) && module<=*max_element(modules.begin(),modules.end()) );
01434       
01435   eta = etaG;
01436   phi = phiG;
01437   channelID=5*(strip-1) + xtal-1; 
01438   towerID=tower;
01439   
01440   vector<int> apdRefChan=ME::apdRefChannels(module, lmr);      
01441   for (unsigned int iref=0;iref<nRefChan;iref++){
01442     if(channelID==apdRefChan[iref] && towerID==apdRefTT 
01443        && apdRefMap[iref].count(module)==0){
01444       apdRefMap[iref][module]=channel;
01445     }
01446   }
01447   
01448   if(isFirstChanModFilled[module-1]==0) {
01449     firstChanMod[module-1]=channel;
01450     isFirstChanModFilled[module-1]=1;
01451   } 
01452   
01453   iEta[channel]=eta;
01454   iPhi[channel]=phi;
01455   iModule[channel]= module ;
01456   iTowerID[channel]=towerID;
01457   iChannelID[channel]=channelID;
01458   idccID[channel]=dccID;
01459   iside[channel]=side;
01460 
01461 }
01462 
01463 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){
01464   
01465   side=MEEEGeom::side(iX, iY, iZ);
01466 
01467   assert( module>=*min_element(modules.begin(),modules.end()) 
01468           && module<=*max_element(modules.begin(),modules.end()) );
01469   
01470   eta = etaG;
01471   phi = phiG;
01472   channelID=ch;
01473   towerID=tower;
01474   
01475   vector<int> apdRefChan=ME::apdRefChannels(module, lmr);      
01476   for (unsigned int iref=0;iref<nRefChan;iref++){
01477     if(channelID==apdRefChan[iref] && towerID==apdRefTT 
01478        && apdRefMap[iref].count(module)==0){
01479       apdRefMap[iref][module]=channel;
01480     }
01481   }
01482   
01483   if(isFirstChanModFilled[module-1]==0) {
01484     firstChanMod[module-1]=channel;
01485     isFirstChanModFilled[module-1]=1;
01486   } 
01487   
01488   iEta[channel]=eta;
01489   iPhi[channel]=phi;
01490   iModule[channel]= module ;
01491   iTowerID[channel]=towerID;
01492   iChannelID[channel]=channelID;
01493   idccID[channel]=dccID;
01494   iside[channel]=side;
01495 
01496 }
01497 DEFINE_FWK_MODULE(EcalLaserAnalyzer2);
01498
CMSSW_4_4_3_patch1/src/CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalLaserAnalyzer2.cc
