HIP Killing and PreMixing for 81X

DataFormats/StdDictionaries/src/classes_def_others.xml

@@ -5,6 +5,7 @@
  <class name="std::allocator<short>"/>
  <class name="std::basic_string<char>"/>
  <class name="std::bidirectional_iterator_tag"/>
+ <class name="std::bitset<6>"/>
  <class name="std::bitset<7>"/>
  <class name="std::bitset<15>"/>
  <class name="std::deque<int>"/>

DataFormats/StdDictionaries/src/classes_def_vector.xml

@@ -43,4 +43,6 @@
  <class name="std::vector<unsigned long long>"/>
  <class name="std::vector<unsigned short>"/>
  <class name="std::vector<const void*>" />
+ <class name="std::vector<std::pair<int,std::bitset<6> > >" />
+ <class name="std::pair<int,std::bitset<6>>"/>
 </lcgdict>

DataFormats/StdDictionaries/src/classes_others.h

@@ -15,6 +15,7 @@ namespace DataFormats_StdDictionaries {
   std::allocator<short> ashort;
   std::basic_string<char> bschar;
   std::bidirectional_iterator_tag bidirectiter;
+  std::bitset<6> dummybitset6;
   std::bitset<7> dummybitset7;
   std::bitset<15> dummybitset15;
   std::deque<int> dummy18;

DataFormats/StdDictionaries/src/classes_vector.h

@@ -65,5 +65,6 @@ namespace DataFormats_StdDictionaries {
   std::vector<unsigned short>::iterator itus;
   std::vector<void *>::iterator itvp;
   std::vector<const void *> dummyVCPtr;
+  std::vector<std::pair<int,std::bitset<6> > > v_p_i_b;
   };
 }

DataFormats/WrappedStdDictionaries/src/classes.h

@@ -6,6 +6,7 @@
 #include <deque>
 #include <set>
 #include <string>
+#include <bitset>
 
 namespace DataFormats_WrappedStdDictionaries {
   struct dictionary {
@@ -39,6 +40,8 @@ namespace DataFormats_WrappedStdDictionaries {
 
   edm::Wrapper<std::set<int> > dummy19;
 
+  edm::Wrapper<std::bitset<6>> dummybitset;
+
   edm::Wrapper<std::pair<unsigned long, unsigned long> > dymmywp1;
   edm::Wrapper<std::pair<unsigned int, unsigned int> > dymmywp2;
   edm::Wrapper<std::pair<unsigned int, int> > dymmywp2_1;
@@ -86,5 +89,6 @@ namespace DataFormats_WrappedStdDictionaries {
   edm::Wrapper<bool> dummyw13;
   edm::Wrapper<unsigned long long> dummyw14;
   edm::Wrapper<long long> dummyw15;
+  edm::Wrapper<std::vector<std::pair<int,std::bitset<6> > > > v_p_i_b;
 };
 }

DataFormats/WrappedStdDictionaries/src/classes_def.xml

@@ -52,4 +52,5 @@
  <class name="edm::Wrapper<unsigned long>"/>
  <class name="edm::Wrapper<unsigned long long>"/>
  <class name="edm::Wrapper<unsigned short>"/>
+ <class name="edm::Wrapper<std::vector<std::pair<int,std::bitset<6> > > >" />
 </lcgdict>

SimGeneral/Configuration/python/SimGeneral_EventContent_cff.py

@@ -11,6 +11,7 @@
 SimGeneralRAW = cms.PSet(
     outputCommands = cms.untracked.vstring('keep CrossingFramePlaybackInfoNew_*_*_*',
                                            'keep PileupSummaryInfos_*_*_*',
+                                           'keep int6stdbitsetstdpairs_*_AffectedAPVList_*',
                                            'keep int_*_bunchSpacing_*')
 )
 #RECO content

SimGeneral/DataMixingModule/plugins/DataMixingModule.h

@@ -59,7 +59,7 @@ namespace edm {
       virtual ~DataMixingModule();
 
       // copies, with EventSetup
-      virtual void checkSignal(const edm::Event &e) {}
+      virtual void checkSignal(const edm::Event &e) override {}; 
       virtual void createnewEDProduct() {}
       virtual void addSignals(const edm::Event &e, const edm::EventSetup& ES); 
       virtual void doPileUp(edm::Event &e,const edm::EventSetup& ES) override;

SimGeneral/DataMixingModule/plugins/DataMixingSiStripMCDigiWorker.cc

@@ -17,6 +17,7 @@
 #include "CondFormats/SiStripObjects/interface/SiStripBadStrip.h"
 #include "CalibTracker/Records/interface/SiStripDependentRecords.h"
 //
+#include "CLHEP/Random/RandFlat.h"
 //
 #include "DataMixingSiStripMCDigiWorker.h"
 
@@ -34,9 +35,11 @@ namespace edm
 							       edm::ConsumesCollector && iC) : 
     label_(ps.getParameter<std::string>("Label")),
     gainLabel(ps.getParameter<std::string>("Gain")),
+    SingleStripNoise(ps.getParameter<bool>("SingleStripNoise")),
     peakMode(ps.getParameter<bool>("APVpeakmode")),
     theThreshold(ps.getParameter<double>("NoiseSigmaThreshold")),
     theElectronPerADC(ps.getParameter<double>( peakMode ? "electronPerAdcPeak" : "electronPerAdcDec" )),
+    APVSaturationFromHIP_(ps.getParameter<bool>("APVSaturationFromHIP")),
     theFedAlgo(ps.getParameter<int>("FedAlgorithm_PM")),
     geometryType(ps.getParameter<std::string>("GeometryType")),
     theSiZeroSuppress(new SiStripFedZeroSuppression(theFedAlgo)),
@@ -53,7 +56,14 @@ namespace edm
     SiStripPileInputTag_ = ps.getParameter<edm::InputTag>("SiStripPileInputTag");
 
     SiStripDigiCollectionDM_  = ps.getParameter<std::string>("SiStripDigiCollectionDM");
+    SistripAPVListDM_= ps.getParameter<std::string>("SiStripAPVListDM");
 
+
+    if(APVSaturationFromHIP_) { 
+      SistripAPVLabelSig_ = ps.getParameter<edm::InputTag>("SistripAPVLabelSig");
+      SiStripAPVPileInputTag_ = ps.getParameter<edm::InputTag>("SistripAPVPileInputTag");
+      iC.consumes< std::vector<std::pair<int,std::bitset<6>> > >(SistripAPVLabelSig_);
+    }
     iC.consumes<edm::DetSetVector<SiStripDigi>>(SistripLabelSig_);
     // clear local storage for this event                                                                     
     SiHitStorage_.clear();
@@ -97,8 +107,6 @@ namespace edm
 	DMinitializeDetUnit(stripdet, iSetup);
       }
     }
-
-
   }
 
 
@@ -113,8 +121,12 @@ namespace edm
     SiStripBadStrip::Range detBadStripRange = deadChannelHandle->getRange(detId);
     //storing the bad strip of the the module. the module is not removed but just signal put to 0
     std::vector<bool>& badChannels = allBadChannels[detId];
+    std::vector<bool>& hipChannels = allHIPChannels[detId];
     badChannels.clear();
     badChannels.insert(badChannels.begin(), numStrips, false);
+    hipChannels.clear();
+    hipChannels.insert(hipChannels.begin(), numStrips, false);
+
     for(SiStripBadStrip::ContainerIterator it = detBadStripRange.first; it != detBadStripRange.second; ++it) {
       SiStripBadStrip::data fs = deadChannelHandle->decode(*it);
       for(int strip = fs.firstStrip; strip < fs.firstStrip + fs.range; ++strip) badChannels[strip] = true;
@@ -148,6 +160,20 @@ namespace edm
       }
 
     }
+
+    // keep here for future reference.  In current implementation, HIP killing is done once in PU file
+    /*  if(APVSaturationFromHIP_) {
+      Handle<std::vector<std::pair<int,std::bitset<6>> > >  APVinput;
+
+      if( e.getByLabel(SistripAPVLabelSig_,APVinput) ) {
+
+	std::vector<std::pair<int,std::bitset<6>> >::const_iterator entry = APVinput->begin();
+	for( ; entry != APVinput->end(); entry++) {
+	  theAffectedAPVmap_.insert(APVMap::value_type(entry->first, entry->second));
+	}
+      }
+      } */
+
   } // end of addSiStripSignals
 
 
@@ -204,6 +230,21 @@ namespace edm
 	  SiHitStorage_.insert( SiGlobalIndex::value_type( DSViter->id, LocalMap ) );
 	}
       }
+
+      if(APVSaturationFromHIP_) {
+	std::shared_ptr<Wrapper<std::vector<std::pair<int,std::bitset<6>> >  >  const> inputAPVPTR =
+	  getProductByTag< std::vector<std::pair<int,std::bitset<6>> > >(*ep, SiStripAPVPileInputTag_, mcc);
+
+	if(inputAPVPTR) {
+
+	  const std::vector<std::pair<int,std::bitset<6>> >   *APVinput = const_cast< std::vector<std::pair<int,std::bitset<6>> >  * >(inputAPVPTR->product());
+
+	  std::vector<std::pair<int,std::bitset<6>> >::const_iterator entry = APVinput->begin();
+	  for( ; entry != APVinput->end(); entry++) {
+	    theAffectedAPVmap_.insert( APVMap::value_type(entry->first, entry->second ));
+	  }
+	}
+      }
     }
   }
 
@@ -222,6 +263,13 @@ namespace edm
     iSetup.get<SiStripThresholdRcd>().get(thresholdHandle);
     iSetup.get<SiStripPedestalsRcd>().get(pedestalHandle);
 
+    edm::Service<edm::RandomNumberGenerator> rng;
+    CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
+
+    std::map< int,std::bitset<6>>  DeadAPVList;
+    DeadAPVList.clear();
+
+
     // First, have to convert all ADC counts to raw pulse heights so that values can be added properly
     // In PreMixing, pulse heights are saved with ADC = sqrt(9.0*PulseHeight) - have to undo. 
 
@@ -275,7 +323,57 @@ namespace edm
       SiRawDigis_.insert( SiGlobalRawIndex::value_type( detID, LocalRawMap ) );
     }
 
-    //  Ok, done with merging raw signals - now add signals on duplicate strips 
+    // If we are killing APVs, merge list of dead ones before we digitize
+
+    int NumberOfBxBetweenHIPandEvent=1e3;
+
+    if(APVSaturationFromHIP_) { 
+
+      // calculate affected BX parameter
+
+      bool HasAtleastOneAffectedAPV=false;
+      while(!HasAtleastOneAffectedAPV){
+	for(int bx=floor(300.0/25.0);bx>0;bx--){ //Reminder: make these numbers not hard coded!!
+	  float temp=CLHEP::RandFlat::shoot(engine)<0.5?1:0;
+	  if(temp==1 && bx<NumberOfBxBetweenHIPandEvent){
+	    NumberOfBxBetweenHIPandEvent=bx;
+	    HasAtleastOneAffectedAPV=true;
+	  }
+	}
+      }
+
+      APVMap::const_iterator iAPVchk;
+      uint32_t formerID = 0;
+      uint32_t currentID;
+      std::bitset<6> NewAPVBits;
+
+      for(APVMap::const_iterator iAPV  = theAffectedAPVmap_.begin();
+	  iAPV != theAffectedAPVmap_.end(); ++iAPV) {
+
+	currentID = iAPV->first;
+
+	if (currentID == formerID) { // we have to OR these
+	  for( int ibit=0; ibit<6; ++ibit){
+	    NewAPVBits[ibit] = NewAPVBits[ibit]||(iAPV->second)[ibit];
+	  }
+	}
+	else {
+	  DeadAPVList[currentID]=NewAPVBits;          
+	  //save pointers for next iteration                                 
+	  formerID = currentID;
+	  NewAPVBits = iAPV->second;
+	}
+
+	iAPVchk = iAPV;
+	if((++iAPVchk) == theAffectedAPVmap_.end()) {  //make sure not to lose the last one 
+	  DeadAPVList[currentID]=NewAPVBits;          
+	}
+      }
+
+    }
+    // 
+
+    //  Ok, done with merging raw signals and APVs - now add signals on duplicate strips 
 
     // collection of Digis to put in the event
     std::vector< edm::DetSet<SiStripDigi> > vSiStripDigi;
@@ -324,8 +422,7 @@ namespace edm
         }
 
         iLocalchk = iLocal;
-        if((++iLocalchk) == LocalMap.end()) {  //make sure not to lose the last one  
-
+	if((++iLocalchk) == LocalMap.end()) {  //make sure not to lose the last one  
           Signals.insert( std::make_pair(formerStrip, ADCSum));
         }
       }
@@ -360,7 +457,28 @@ namespace edm
 
 	//removing signal from the dead (and HIP effected) strips
 	std::vector<bool>& badChannels = allBadChannels[detID];
-	for(int strip =0; strip < numStrips; ++strip) if(badChannels[strip]) detAmpl[strip] = 0.;
+
+        for(int strip =0; strip < numStrips; ++strip) {
+          if(badChannels[strip]) detAmpl[strip] = 0.;
+	}
+
+	if(APVSaturationFromHIP_) {
+	  std::bitset<6> & bs=DeadAPVList[detID];
+
+	  if(bs.any()){
+	    // Here below is the scaling function which describes the evolution of the baseline (i.e. how the charge is suppressed).
+	    // This must be replaced as soon as we have a proper modeling of the baseline evolution from VR runs
+	    float Shift=1-NumberOfBxBetweenHIPandEvent/floor(300.0/25.0); //Reminder: make these numbers not hardcoded!! 
+	    float randomX=CLHEP::RandFlat::shoot(engine);
+	    float scalingValue=(randomX-Shift)*10.0/7.0-3.0/7.0;
+
+	    for(int strip =0; strip < numStrips; ++strip) {
+	      if(!badChannels[strip] &&  bs[strip/128]==1){
+		detAmpl[strip] *=scalingValue>0?scalingValue:0.0;
+	      }
+	    }
+	  }
+	}
 
 	SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detID);
 	SiStripApvGain::Range detGainRange = gainHandle->getRange(detID);
@@ -388,17 +506,32 @@ namespace edm
         size_t firstChannelWithSignal = 0;
         size_t lastChannelWithSignal = numStrips;
 
-	int RefStrip = int(numStrips/2.);
-	while(RefStrip<numStrips&&badChannels[RefStrip]){ //if the refstrip is bad, I move up to when I don't find it
-	  RefStrip++;
-	}
-	if(RefStrip<numStrips){
-	  float noiseRMS = noiseHandle->getNoise(RefStrip,detNoiseRange);
-	  float gainValue = gainHandle->getStripGain(RefStrip, detGainRange);
-	  edm::Service<edm::RandomNumberGenerator> rng;
-	  CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
-	  theSiNoiseAdder->addNoise(detAmpl,firstChannelWithSignal,lastChannelWithSignal,numStrips,noiseRMS*theElectronPerADC/gainValue,engine);
-	}
+        if(SingleStripNoise){
+          //      std::cout<<"In SSN, detId="<<detID<<std::endl;                                                                                                   
+	  std::vector<float> noiseRMSv;
+          noiseRMSv.clear();
+          noiseRMSv.insert(noiseRMSv.begin(),numStrips,0.);
+          for(int strip=0; strip< numStrips; ++strip){
+            if(!badChannels[strip]){
+              float gainValue = gainHandle->getStripGain(strip, detGainRange);
+              noiseRMSv[strip] = (noiseHandle->getNoise(strip,detNoiseRange))* theElectronPerADC/gainValue;
+              //std::cout<<"<SiStripDigitizerAlgorithm::digitize>: gainValue: "<<gainValue<<"\tnoiseRMSv["<<strip<<"]: "<<noiseRMSv[strip]<<std::endl;             
+            }
+          }
+          theSiNoiseAdder->addNoiseVR(detAmpl, noiseRMSv, engine);
+        } else {
+          int RefStrip = int(numStrips/2.);
+          while(RefStrip<numStrips&&badChannels[RefStrip]){ //if the refstrip is bad, I move up to when I don't find it                                            
+            RefStrip++;
+          }
+          if(RefStrip<numStrips){
+            float RefgainValue = gainHandle->getStripGain(RefStrip, detGainRange);
+            float RefnoiseRMS = noiseHandle->getNoise(RefStrip,detNoiseRange) *theElectronPerADC/RefgainValue;
+
+            theSiNoiseAdder->addNoise(detAmpl,firstChannelWithSignal,lastChannelWithSignal,numStrips,RefnoiseRMS, engine);
+            //std::cout<<"<SiStripDigitizerAlgorithm::digitize>: RefgainValue: "<<RefgainValue<<"\tRefnoiseRMS: "<<RefnoiseRMS<<std::endl;                         
+          }
+        }
 
 	DigitalVecType digis;
 	theSiZeroSuppress->suppress(theSiDigitalConverter->convert(detAmpl, gainHandle, detID), digis, detID,noiseHandle,thresholdHandle);

SimGeneral/DataMixingModule/plugins/DataMixingSiStripMCDigiWorker.h

@@ -78,6 +78,11 @@ namespace edm
       edm::InputTag SiStripPileInputTag_ ;    // InputTag for pileup strips
       std::string SiStripDigiCollectionDM_  ; // secondary name to be given to new SiStrip digis
 
+      edm::InputTag SistripAPVLabelSig_;      // where to find vector of dead APVs
+      edm::InputTag SiStripAPVPileInputTag_;
+      std::string SistripAPVListDM_;          // output tag
+
+
       // 
 
       typedef float Amplitude;
@@ -110,14 +115,21 @@ namespace edm
 
       signalMaps signals_;
 
+      // to keep track of dead APVs from HIP interactions
+      typedef std::multimap< uint32_t, std::bitset<6> > APVMap;
+
+      APVMap theAffectedAPVmap_;
+
       // for noise adding:
 
       std::string label_;
 
       std::string gainLabel;
+      bool SingleStripNoise;
       bool peakMode;
       double theThreshold;
       double theElectronPerADC;
+      bool APVSaturationFromHIP_;
       int theFedAlgo;
       std::string geometryType;
 
@@ -129,6 +141,8 @@ namespace edm
 
       // bad channels for each detector ID
       std::map<unsigned int, std::vector<bool> > allBadChannels;
+      // channels killed by HIP interactions for each detector ID
+      std::map<unsigned int, std::vector<bool> > allHIPChannels;
       // first and last channel wit signal for each detector ID
       std::map<unsigned int, size_t> firstChannelsWithSignal;
       std::map<unsigned int, size_t> lastChannelsWithSignal;

SimGeneral/DataMixingModule/python/mixOne_simraw_on_sim_cfi.py

@@ -177,7 +177,12 @@
     CSCStripDigiSimLinkPileInputTag = cms.InputTag("simMuonCSCDigis","MuonCSCStripDigiSimLinks"),
     CSCWireDigiSimLinkPileInputTag = cms.InputTag("simMuonCSCDigis","MuonCSCWireDigiSimLinks"),
 
-    #
+    # Dead APV Vector
+    SistripAPVPileInputTag = cms.InputTag("mix","AffectedAPVList"),
+    SistripAPVLabelSig = cms.InputTag("mix","AffectedAPVList"),
+
+    # Note: elements with "@MIXING" in the input tag are generated by
+    # running Raw2Digi in the input step on the Secondary input stream
     EBPileInputTag = cms.InputTag("ecalDigis","ebDigis","@MIXING"),
     EEPileInputTag = cms.InputTag("ecalDigis","eeDigis","@MIXING"),
     ESPileInputTag = cms.InputTag("ecalPreshowerDigis","","@MIXING"),
@@ -231,7 +236,7 @@
     RPCDigiSimLinkDM = cms.string('RPCDigiSimLink'),
     CSCStripDigiSimLinkDM = cms.string('MuonCSCStripDigiSimLinks'),
     CSCWireDigiSimLinkDM = cms.string('MuonCSCWireDigiSimLinks'),
-
+    SiStripAPVListDM = cms.string('SiStripAPVList'),
 
     #
     #  Calorimeter Digis