00001 #include "RecoLocalTracker/SiStripZeroSuppression/plugins/SiStripZeroSuppression.h"
00002
00003 #include "FWCore/Framework/interface/Event.h"
00004 #include "FWCore/Framework/interface/EventSetup.h"
00005 #include "DataFormats/Common/interface/Handle.h"
00006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00007 #include "DataFormats/Common/interface/DetSetVector.h"
00008 #include "DataFormats/SiStripDigi/interface/SiStripDigi.h"
00009 #include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h"
00010 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripRawProcessingFactory.h"
00011 #include <memory>
00012
00013 SiStripZeroSuppression::
00014 SiStripZeroSuppression(edm::ParameterSet const& conf)
00015 : inputTags(conf.getParameter<std::vector<edm::InputTag> >("RawDigiProducersList")),
00016 algorithms(SiStripRawProcessingFactory::create(conf.getParameter<edm::ParameterSet>("Algorithms"))),
00017 storeCM(conf.getParameter<bool>("storeCM")),
00018 mergeCollections(conf.getParameter<bool>("mergeCollections"))
00019
00020 {
00021
00022 produceRawDigis = conf.getParameter<bool>("produceRawDigis");
00023 produceCalculatedBaseline = conf.getParameter<bool>("produceCalculatedBaseline");
00024 produceBaselinePoints = conf.getParameter<bool>("produceBaselinePoints");
00025 storeInZScollBadAPV = conf.getParameter<bool>("storeInZScollBadAPV");
00026 fixCM = conf.getParameter<bool>("fixCM");
00027
00028 if(mergeCollections){
00029 storeCM = false;
00030 produceRawDigis = false;
00031 DigisToMergeZS = conf.getParameter<edm::InputTag>("DigisToMergeZS");
00032 DigisToMergeVR = conf.getParameter<edm::InputTag>("DigisToMergeVR");
00033 produces< edm::DetSetVector<SiStripDigi> > ("ZeroSuppressed");
00034 }
00035
00036 for(tag_iterator_t inputTag = inputTags.begin(); inputTag != inputTags.end(); ++inputTag ){
00037 produces< edm::DetSetVector<SiStripDigi> > (inputTag->instance());
00038 if(produceRawDigis) produces< edm::DetSetVector<SiStripRawDigi> > (inputTag->instance());
00039 if(produceCalculatedBaseline) produces< edm::DetSetVector<SiStripProcessedRawDigi> > ("BADAPVBASELINE"+inputTag->instance());
00040 if(produceBaselinePoints) produces< edm::DetSetVector<SiStripDigi> > ("BADAPVBASELINEPOINTS"+inputTag->instance());
00041 if(storeCM) produces< edm::DetSetVector<SiStripProcessedRawDigi> > ("APVCM"+inputTag->instance());
00042 }
00043
00044
00045
00046
00047
00048 }
00049
00050 void SiStripZeroSuppression::
00051 produce(edm::Event& e, const edm::EventSetup& es) {
00052
00053 algorithms->initialize(es, e);
00054
00055 if(mergeCollections){
00056 this->MergeCollectionsZeroSuppression(e);
00057 }else{
00058 this->StandardZeroSuppression(e);
00059 }
00060 }
00061
00062 inline void SiStripZeroSuppression::StandardZeroSuppression(edm::Event& e){
00063
00064 for(tag_iterator_t inputTag = inputTags.begin(); inputTag != inputTags.end(); ++inputTag ) {
00065
00066 edm::Handle< edm::DetSetVector<SiStripRawDigi> > input;
00067 e.getByLabel(*inputTag,input);
00068
00069 if (input->size())
00070 processRaw(*inputTag, *input);
00071
00072 std::auto_ptr< edm::DetSetVector<SiStripDigi> > output(new edm::DetSetVector<SiStripDigi>(output_base) );
00073 e.put( output, inputTag->instance() );
00074
00075 if(produceRawDigis){
00076 std::auto_ptr< edm::DetSetVector<SiStripRawDigi> > outputraw(new edm::DetSetVector<SiStripRawDigi>(output_base_raw) );
00077 e.put(outputraw, inputTag->instance() );
00078 }
00079
00080 if(produceCalculatedBaseline){
00081 std::auto_ptr< edm::DetSetVector<SiStripProcessedRawDigi> > outputbaseline(new edm::DetSetVector<SiStripProcessedRawDigi>(output_baseline) );
00082 e.put(outputbaseline, "BADAPVBASELINE"+inputTag->instance() );
00083 }
00084
00085 if(produceBaselinePoints){
00086 std::auto_ptr< edm::DetSetVector<SiStripDigi> > outputbaselinepoints(new edm::DetSetVector<SiStripDigi>(output_baseline_points) );
00087 e.put(outputbaselinepoints, "BADAPVBASELINEPOINTS"+inputTag->instance() );
00088 }
00089
00090 if(storeCM){
00091 std::auto_ptr< edm::DetSetVector<SiStripProcessedRawDigi> > outputAPVCM(new edm::DetSetVector<SiStripProcessedRawDigi>(output_apvcm) );
00092 e.put( outputAPVCM,"APVCM"+inputTag->instance());
00093 }
00094
00095 }
00096 }
00097
00098
00099
00100
00101 inline
00102 void SiStripZeroSuppression::
00103 processRaw(const edm::InputTag& inputTag, const edm::DetSetVector<SiStripRawDigi>& input ) {
00104
00105 output_apvcm.clear();
00106 output_baseline.clear();
00107 output_baseline_points.clear();
00108 output_base.clear();
00109 output_base_raw.clear();
00110
00111 if(storeCM) output_apvcm.reserve(16000);
00112 if(produceCalculatedBaseline) output_baseline.reserve(16000);
00113 if(produceBaselinePoints) output_baseline_points.reserve(16000);
00114 if(produceRawDigis) output_base_raw.reserve(16000);
00115 output_base.reserve(16000);
00116
00117
00118 for ( edm::DetSetVector<SiStripRawDigi>::const_iterator
00119 rawDigis = input.begin(); rawDigis != input.end(); ++rawDigis) {
00120
00121 edm::DetSet<SiStripDigi> suppressedDigis(rawDigis->id);
00122 int16_t nAPVflagged = 0;
00123
00124 if ( "ProcessedRaw" == inputTag.instance()) nAPVflagged = algorithms->SuppressProcessedRawData(*rawDigis, suppressedDigis);
00125 else if ( "VirginRaw" == inputTag.instance()) nAPVflagged = algorithms->SuppressVirginRawData(*rawDigis, suppressedDigis);
00126 else
00127 throw cms::Exception("Unknown input type")
00128 << inputTag.instance() << " unknown. SiStripZeroZuppression can only process types \"VirginRaw\" and \"ProcessedRaw\" ";
00129
00130
00131 this->storeExtraOutput(rawDigis->id, nAPVflagged);
00132 if (suppressedDigis.size() && (storeInZScollBadAPV || nAPVflagged ==0))
00133 output_base.push_back(suppressedDigis);
00134
00135 if (produceRawDigis && nAPVflagged > 0){
00136 edm::DetSet<SiStripRawDigi> outRawDigis(rawDigis->id);
00137 this->formatRawDigis(rawDigis, outRawDigis);
00138 output_base_raw.push_back(outRawDigis);
00139 }
00140
00141 }
00142
00143 }
00144
00145 inline
00146 void SiStripZeroSuppression::formatRawDigis(edm::DetSetVector<SiStripRawDigi>::const_iterator
00147 rawDigis, edm::DetSet<SiStripRawDigi>& outRawDigis){
00148
00149 const std::vector<bool>& apvf = algorithms->GetAPVFlags();
00150 edm::DetSet<SiStripRawDigi>::const_iterator itRawDigis = rawDigis->begin();
00151
00152 uint32_t strip=0;
00153 for (; itRawDigis != rawDigis->end(); ++itRawDigis){
00154 int16_t APVn = strip/128;
00155 if(apvf[APVn]) outRawDigis.push_back(*itRawDigis);
00156 else outRawDigis.push_back(SiStripRawDigi(0));
00157 ++strip;
00158 }
00159
00160 }
00161
00162
00163 inline
00164 void SiStripZeroSuppression::storeExtraOutput(uint32_t id, int16_t nAPVflagged){
00165
00166 const std::vector< std::pair<short,float> >& vmedians = algorithms->getAPVsCM();
00167 if(storeCM) this->storeCMN(id, vmedians);
00168 if(nAPVflagged > 0){
00169 if(produceCalculatedBaseline) this->storeBaseline(id, vmedians);
00170 if(produceBaselinePoints) this->storeBaselinePoints(id);
00171 }
00172 }
00173
00174
00175 inline
00176 void SiStripZeroSuppression::storeBaseline(uint32_t id, const std::vector< std::pair<short,float> >& vmedians){
00177
00178 std::map< uint16_t, std::vector < int16_t> >& baselinemap = algorithms->GetBaselineMap();
00179
00180 edm::DetSet<SiStripProcessedRawDigi> baselineDetSet(id);
00181 std::map< uint16_t, std::vector < int16_t> >::iterator itBaselineMap;
00182
00183 for(size_t i=0;i<vmedians.size();++i){
00184 uint16_t APVn = vmedians[i].first;
00185 float median = vmedians[i].second;
00186 itBaselineMap = baselinemap.find(APVn);
00187 if(itBaselineMap==baselinemap.end()){
00188 for(size_t strip=0; strip < 128; ++strip) baselineDetSet.push_back(SiStripProcessedRawDigi(median));
00189 } else {
00190 for(size_t strip=0; strip < 128; ++strip) baselineDetSet.push_back(SiStripProcessedRawDigi((itBaselineMap->second)[strip]));
00191 }
00192
00193 }
00194
00195 if(baselineDetSet.size())
00196 output_baseline.push_back(baselineDetSet);
00197
00198 }
00199
00200 inline
00201 void SiStripZeroSuppression::storeBaselinePoints(uint32_t id){
00202
00203 std::map< uint16_t, std::map< uint16_t, int16_t> >& BasPointVec = algorithms->GetSmoothedPoints();
00204 std::map< uint16_t, std::map< uint16_t, int16_t> >::iterator itBasPointVect = BasPointVec.begin() ;
00205 std::map< uint16_t, int16_t>::iterator itBaselinePointMap;
00206
00207 edm::DetSet<SiStripDigi> baspointDetSet(id);
00208
00209 for(; itBasPointVect != BasPointVec.end(); ++itBasPointVect){
00210 uint16_t APVn= itBasPointVect->first;
00211 itBaselinePointMap =itBasPointVect->second.begin();
00212 for(;itBaselinePointMap != itBasPointVect->second.end(); ++itBaselinePointMap){
00213 uint16_t bpstrip = (itBaselinePointMap->first) + APVn*128;
00214 int16_t bp = itBaselinePointMap->second;
00215 baspointDetSet.push_back(SiStripDigi(bpstrip,bp+128));
00216
00217 }
00218 }
00219
00220
00221 if(baspointDetSet.size())
00222 output_baseline_points.push_back(baspointDetSet);
00223
00224 }
00225
00226 inline
00227 void SiStripZeroSuppression::storeCMN(uint32_t id, const std::vector< std::pair<short,float> >& vmedians){
00228
00229 edm::DetSet<SiStripProcessedRawDigi> apvDetSet(id);
00230 short apvNb=0;
00231
00232 std::vector<bool> apvf;
00233 apvf.clear();
00234 apvf.insert(apvf.begin(), 6, false);
00235
00236 if(fixCM){
00237 std::vector<bool>& apvFlagged = algorithms->GetAPVFlags();
00238 for(uint16_t it=0; it< apvFlagged.size(); ++it) apvf[it] = apvFlagged[it];
00239 }
00240
00241 for(size_t i=0;i<vmedians.size();++i){
00242 if(vmedians[i].first>apvNb){
00243 for(int i=0;i<vmedians[i].first-apvNb;++i){
00244 apvDetSet.push_back(SiStripProcessedRawDigi(-999.));
00245 apvNb++;
00246 }
00247 }
00248
00249 if(fixCM&&apvf[vmedians[i].first]){
00250 apvDetSet.push_back(SiStripProcessedRawDigi(-999.));
00251 }else{
00252 apvDetSet.push_back(SiStripProcessedRawDigi(vmedians[i].second));
00253 }
00254 apvNb++;
00255 }
00256
00257 if(apvDetSet.size())
00258 output_apvcm.push_back(apvDetSet);
00259
00260 }
00261
00262
00263 inline void SiStripZeroSuppression::MergeCollectionsZeroSuppression(edm::Event& e){
00264
00265 std::cout<< "starting Merging" << std::endl;
00266 edm::Handle< edm::DetSetVector<SiStripDigi> > inputdigi;
00267 edm::Handle< edm::DetSetVector<SiStripRawDigi> > inputraw;
00268 e.getByLabel(DigisToMergeZS,inputdigi);
00269 e.getByLabel(DigisToMergeVR,inputraw);
00270
00271 std::cout << inputdigi->size() << " " << inputraw->size() << std::endl;
00272
00273 if (inputraw->size()){
00274
00275 std::vector<edm::DetSet<SiStripDigi> > outputdigi;
00276 outputdigi.clear();
00277
00278
00279 for ( edm::DetSetVector<SiStripDigi>::const_iterator Digis = inputdigi->begin(); Digis != inputdigi->end(); ++Digis) outputdigi.push_back(*Digis);
00280
00281
00282 std::cout << "looping over the raw data collection" << std::endl;
00283 for ( edm::DetSetVector<SiStripRawDigi>::const_iterator rawDigis = inputraw->begin(); rawDigis != inputraw->end(); ++rawDigis) {
00284
00285 edm::DetSet<SiStripRawDigi>::const_iterator itRawDigis = rawDigis->begin();
00286 uint16_t nAPV = rawDigis->size()/128;
00287 uint32_t rawDetId = rawDigis->id;
00288
00289 std::vector<bool> restoredAPV;
00290 restoredAPV.clear();
00291 restoredAPV.insert(restoredAPV.begin(), nAPV, false);
00292
00293
00294 bool isModuleRestored = false;
00295 for( uint16_t strip =0; strip < rawDigis->size();++strip){
00296 if(itRawDigis[strip].adc()!=0){
00297 restoredAPV[strip/128] = true;
00298 isModuleRestored = true;
00299 }
00300 }
00301
00302
00303 if(isModuleRestored){
00304 std::cout << "apply the ZS to the raw data collection" << std::endl;
00305 edm::DetSet<SiStripDigi> suppressedDigis(rawDetId);
00306 std::vector<int16_t> processedRawDigis(rawDigis->size());
00307 algorithms->SuppressVirginRawData(*rawDigis, suppressedDigis);
00308
00309 if(suppressedDigis.size()){
00310 std::cout << "looking for the detId with the new ZS in the collection of the zero suppressed data" << std::endl;
00311 std::vector<edm::DetSet<SiStripDigi> >::iterator zsModule = outputdigi.begin();
00312
00313
00314 uint32_t zsDetId = zsModule->id;
00315 bool isModuleInZscollection = false;
00316 while((zsDetId <= rawDetId)&&(zsModule != outputdigi.end())&&(!isModuleInZscollection)){
00317
00318 if( zsDetId == rawDetId){
00319 isModuleInZscollection = true;
00320 }else{
00321 ++zsModule;
00322 zsDetId = zsModule->id;
00323 }
00324 }
00325 std::cout << "after the look " << rawDetId << " ==== " << zsDetId << std::endl;
00326 std::cout << "exiting looking for the detId with the new ZS in the collection of the zero suppressed data" << std::endl;
00327
00328
00329 std::vector<uint16_t> MergedRawDigis;
00330 MergedRawDigis.clear();
00331 MergedRawDigis.insert(MergedRawDigis.begin(), nAPV*128, 0);
00332
00333 uint32_t count=0;
00334
00335 edm::DetSet<SiStripDigi> newDigiToIndert(rawDetId);
00336 if(!isModuleInZscollection){
00337 std::cout << "WE HAVE A PROBLEM, THE MODULE IS NTOT FOUND" << std::endl;
00338 outputdigi.insert(zsModule, newDigiToIndert);
00339 --zsModule;
00340 std::cout << "new module id -1 " << zsModule->id << std::endl;
00341 ++zsModule;
00342 std::cout << "new module id " << zsModule->id << std::endl;
00343 ++zsModule;
00344 std::cout << "new module id +1 " << zsModule->id << std::endl;
00345 --zsModule;
00346
00347 } else {
00348 std::cout << "inserting only the digis for not restored APVs" << std::endl;
00349 std::cout << "size : " << zsModule->size() << std::endl;
00350 edm::DetSet<SiStripDigi>::iterator itZsModule = zsModule->begin();
00351 for(; itZsModule != zsModule->end(); ++itZsModule){
00352 uint16_t adc = itZsModule->adc();
00353 uint16_t strip = itZsModule->strip();
00354 if(!restoredAPV[strip/128]){
00355 MergedRawDigis[strip] = adc;
00356 ++count;
00357 std::cout << "original count: "<< count << " strip: " << strip << " adc: " << adc << std::endl;
00358 }
00359 }
00360
00361 }
00362
00363 std::cout << "size of digis to keep: " << count << std::endl;
00364 std::cout << "inserting only the digis for the restored APVs" << std::endl;
00365 std::cout << "size : " << suppressedDigis.size() << std::endl;
00366 edm::DetSet<SiStripDigi>::iterator itSuppDigi = suppressedDigis.begin();
00367 for(; itSuppDigi != suppressedDigis.end(); ++itSuppDigi){
00368 uint16_t adc = itSuppDigi->adc();
00369 uint16_t strip = itSuppDigi->strip();
00370 if(restoredAPV[strip/128]){
00371 MergedRawDigis[strip] = adc;
00372 std::cout << "new suppressed strip: " << strip << " adc: " << adc << std::endl;
00373 }
00374 }
00375
00376
00377
00378 std::cout << "suppressing the raw digis" << std::endl;
00379 zsModule->clear();
00380 for(uint16_t strip=0; strip < MergedRawDigis.size(); ++strip){
00381 uint16_t adc = MergedRawDigis[strip];
00382 if(adc) zsModule->push_back(SiStripDigi(strip, adc));
00383 }
00384 std::cout << "size zsModule after the merging: " << zsModule->size() << std::endl;
00385 if((count + suppressedDigis.size()) != zsModule->size()) std::cout << "WE HAVE A PROBLEM!!!! THE NUMBER OF DIGIS IS NOT RIGHT==============" << std::endl;
00386 std::cout << "exiting suppressing the raw digis" << std::endl;
00387 }
00388 }
00389 }
00390
00391 uint32_t oldid =0;
00392 for(edm::DetSetVector<SiStripDigi>::const_iterator dg = outputdigi.begin(); dg != outputdigi.end(); ++dg){
00393 uint32_t iddg = dg->id;
00394 if(iddg < oldid){
00395 std::cout<< "NOT IN THE RIGHT ORGER" << std:: endl;
00396 std::cout<< "======================="<< std:: endl;
00397 }
00398 oldid = iddg;
00399 }
00400
00401
00402 std::cout << "write the output vector" << std::endl;
00403 std::auto_ptr< edm::DetSetVector<SiStripDigi> > output(new edm::DetSetVector<SiStripDigi>(outputdigi) );
00404 e.put( output, "ZeroSuppressed" );
00405
00406
00407 }
00408
00409
00410
00411 }
00412
00413 inline void SiStripZeroSuppression::CollectionMergedZeroSuppression(edm::Event& e){
00414
00415 for(tag_iterator_t inputTag = inputTags.begin(); inputTag != inputTags.end(); ++inputTag ) {
00416
00417 edm::Handle< edm::DetSetVector<SiStripDigi> > inputdigi;
00418 edm::Handle< edm::DetSetVector<SiStripRawDigi> > inputraw;
00419 e.getByLabel(*inputTag,inputdigi);
00420 e.getByLabel(*inputTag,inputraw);
00421
00422 std::vector<edm::DetSet<SiStripDigi> > outputdigi;
00423 std::vector<edm::DetSet<SiStripRawDigi> > outputraw;
00424
00425
00426 if (inputraw->size())
00427 processRaw(*inputTag, *inputraw);
00428
00429
00430 for ( std::vector<edm::DetSet<SiStripDigi> >::const_iterator itinputdigi = inputdigi->begin(); itinputdigi !=inputdigi->end(); ++itinputdigi) {
00431 output_base.push_back(*itinputdigi);
00432 }
00433
00434 std::auto_ptr< edm::DetSetVector<SiStripDigi> > output(new edm::DetSetVector<SiStripDigi>(output_base) );
00435 e.put( output, inputTag->instance() );
00436
00437 }
00438
00439 }
1.7.3