Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:11:18

0001 // user include files
0002 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0003 
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 
0009 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0010 
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 
0013 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0014 #include "FastSimulation/Event/interface/FSimEvent.h"
0015 #include "FastSimulation/Event/interface/FSimTrack.h"
0016 #include "FastSimulation/Event/interface/FSimVertex.h"
0017 #include "FastSimDataFormats/NuclearInteractions/interface/NUEvent.h"
0018 
0019 #include <memory>
0020 #include <string>
0021 #include <vector>
0022 #include "TFile.h"
0023 #include "TTree.h"
0024 #include "TProcessID.h"
0025 
0026 class testNuclearInteractions : public DQMEDAnalyzer {
0027 public:
0028   explicit testNuclearInteractions(const edm::ParameterSet&);
0029   ~testNuclearInteractions();
0030 
0031   virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
0032   virtual void dqmBeginRun(edm::Run const&, const edm::EventSetup&) override;
0033   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0034 
0035 private:
0036   // See RecoParticleFlow/PFProducer/interface/PFProducer.h
0037   const edm::ParameterSet particleFilter_;
0038   const bool saveNU;
0039   const int maxNU;
0040   const edm::ESGetToken<HepPDT::ParticleDataTable, PDTRecord> tok_pdt_;
0041   const edm::EDGetTokenT<std::vector<SimTrack>> tok_fullSimTk_;
0042   const edm::EDGetTokenT<std::vector<SimVertex>> tok_fullSimVx_;
0043   const edm::EDGetTokenT<std::vector<SimTrack>> tok_fastSimTk_;
0044   const edm::EDGetTokenT<std::vector<SimVertex>> tok_fastSimVx_;
0045   std::vector<std::unique_ptr<FSimEvent>> mySimEvent;
0046   NUEvent* nuEvent;
0047   TTree* nuTree;
0048   TFile* outFile;
0049   int ObjectNumber;
0050   std::string simModuleLabel_;
0051   // Histograms
0052 
0053   std::vector<MonitorElement*> h0;
0054   std::vector<MonitorElement*> h1;
0055   std::vector<MonitorElement*> h2;
0056   std::vector<MonitorElement*> h3;
0057   std::vector<MonitorElement*> h4;
0058   std::vector<MonitorElement*> h5;
0059   std::vector<MonitorElement*> h6;
0060   std::vector<MonitorElement*> h7;
0061   std::vector<MonitorElement*> h8;
0062   std::vector<MonitorElement*> h9;
0063   std::vector<MonitorElement*> h10;
0064   std::vector<MonitorElement*> htmp;
0065   std::vector<MonitorElement*> totalCharge;
0066 
0067   std::vector<std::vector<MonitorElement*>> h100;
0068   std::vector<std::vector<MonitorElement*>> h200;
0069   std::vector<std::vector<MonitorElement*>> h300;
0070   std::vector<std::vector<double>> trackerRadius;
0071   std::vector<std::vector<double>> trackerLength;
0072   std::vector<std::vector<double>> blockTrackerRadius;
0073   std::vector<std::vector<double>> blockTrackerLength;
0074   std::vector<std::vector<double>> subTrackerRadius;
0075   std::vector<std::vector<double>> subTrackerLength;
0076   std::vector<double> tmpRadius;
0077   std::vector<double> tmpLength;
0078 
0079   std::vector<unsigned> stoppedPions;
0080   std::vector<unsigned> interactingPions;
0081   /*
0082   std::vector<MonitorElement*> h6;
0083   std::vector<MonitorElement*> h7;
0084   std::vector<MonitorElement*> h8;
0085   std::vector<MonitorElement*> h9;
0086   std::vector<MonitorElement*> h10;
0087   std::vector<MonitorElement*> h11;
0088   std::vector<MonitorElement*> h12;
0089   */
0090   int intfull;
0091   int intfast;
0092 
0093   std::string NUEventFileName;
0094 
0095   int totalNEvt;
0096   int totalNU;
0097 };
0098 
0099 testNuclearInteractions::testNuclearInteractions(const edm::ParameterSet& p)
0100     : particleFilter_(p.getParameter<edm::ParameterSet>("TestParticleFilter")),
0101       saveNU(p.getParameter<bool>("SaveNuclearInteractions")),
0102       maxNU(p.getParameter<unsigned>("MaxNumberOfNuclearInteractions")),
0103       tok_pdt_(esConsumes<HepPDT::ParticleDataTable, PDTRecord>()),
0104       tok_fullSimTk_(consumes<std::vector<SimTrack>>(edm::InputTag("g4SimHits"))),
0105       tok_fullSimVx_(consumes<std::vector<SimVertex>>(edm::InputTag("g4SimHits"))),
0106       tok_fastSimTk_(consumes<std::vector<SimTrack>>(edm::InputTag("famosSimHits"))),
0107       tok_fastSimVx_(consumes<std::vector<SimVertex>>(edm::InputTag("famosSimHits"))),
0108       h0(2, static_cast<MonitorElement*>(0)),
0109       h1(2, static_cast<MonitorElement*>(0)),
0110       h2(2, static_cast<MonitorElement*>(0)),
0111       h3(2, static_cast<MonitorElement*>(0)),
0112       h4(2, static_cast<MonitorElement*>(0)),
0113       h5(2, static_cast<MonitorElement*>(0)),
0114       h6(2, static_cast<MonitorElement*>(0)),
0115       h7(2, static_cast<MonitorElement*>(0)),
0116       h8(2, static_cast<MonitorElement*>(0)),
0117       h9(2, static_cast<MonitorElement*>(0)),
0118       h10(2, static_cast<MonitorElement*>(0)),
0119       htmp(2, static_cast<MonitorElement*>(0)),
0120       totalCharge(2, static_cast<MonitorElement*>(0)),
0121       tmpRadius(2, static_cast<double>(0.)),
0122       tmpLength(2, static_cast<double>(0.)),
0123       stoppedPions(2, static_cast<unsigned>(0)),
0124       interactingPions(2, static_cast<unsigned>(0)),
0125       /*
0126   h6(2,static_cast<MonitorElement*>(0)),
0127   h7(2,static_cast<MonitorElement*>(0)),
0128   h8(2,static_cast<MonitorElement*>(0)),
0129   h9(2,static_cast<MonitorElement*>(0)),
0130   h10(2,static_cast<MonitorElement*>(0)),
0131   h11(2,static_cast<MonitorElement*>(0)),
0132   h12(2,static_cast<MonitorElement*>(0)),
0133  */
0134       intfull(0),
0135       intfast(0),
0136       totalNEvt(0),
0137       totalNU(0) {
0138   // Let's just initialize the SimEvent's
0139   // For the full sim
0140   mySimEvent.emplace_back(std::make_unique<FSimEvent>(particleFilter_));
0141   // For the fast sim
0142   mySimEvent.emplace_back(std::make_unique<FSimEvent>(particleFilter_));
0143 
0144   // Do we save the nuclear interactions?
0145   if (saveNU)
0146     std::cout << "Nuclear Interactions will be saved ! " << std::endl;
0147 
0148   // Where the nuclear interactions are saved;
0149   NUEventFileName = "none";
0150   if (saveNU) {
0151     nuEvent = new NUEvent();
0152 
0153     NUEventFileName = p.getUntrackedParameter<std::string>("NUEventFile", "NuclearInteractionsTest.root");
0154     //    std::string outFileName = "NuclearInteractionsTest.root";
0155     outFile = new TFile(NUEventFileName.c_str(), "recreate");
0156 
0157     // Open the tree
0158     nuTree = new TTree("NuclearInteractions", "");
0159     nuTree->Branch("nuEvent", "NUEvent", &nuEvent, 32000, 99);
0160   }
0161 
0162   // ObjectNumber
0163   ObjectNumber = -1;
0164 
0165   // Beam Pipe
0166   std::vector<double> tmpRadius = p.getUntrackedParameter<std::vector<double>>("BPCylinderRadius");
0167   std::vector<double> tmpLength = p.getUntrackedParameter<std::vector<double>>("BPCylinderLength");
0168   trackerRadius.push_back(tmpRadius);
0169   trackerLength.push_back(tmpLength);
0170 
0171   // Beam Pipe (cont'd)
0172   subTrackerRadius.push_back(tmpRadius);
0173   subTrackerLength.push_back(tmpLength);
0174 
0175   // PIXB1
0176   tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXB1CylinderRadius");
0177   tmpLength = p.getUntrackedParameter<std::vector<double>>("PXB1CylinderLength");
0178   trackerRadius.push_back(tmpRadius);
0179   trackerLength.push_back(tmpLength);
0180 
0181   // PIXB2
0182   tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXB2CylinderRadius");
0183   tmpLength = p.getUntrackedParameter<std::vector<double>>("PXB2CylinderLength");
0184   trackerRadius.push_back(tmpRadius);
0185   trackerLength.push_back(tmpLength);
0186 
0187   // PIXB3
0188   tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXB3CylinderRadius");
0189   tmpLength = p.getUntrackedParameter<std::vector<double>>("PXB3CylinderLength");
0190   trackerRadius.push_back(tmpRadius);
0191   trackerLength.push_back(tmpLength);
0192 
0193   // PIXB Cables
0194   tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXBCablesCylinderRadius");
0195   tmpLength = p.getUntrackedParameter<std::vector<double>>("PXBCablesCylinderLength");
0196   trackerRadius.push_back(tmpRadius);
0197   trackerLength.push_back(tmpLength);
0198 
0199   // All Pixel Barrel
0200   blockTrackerRadius.push_back(tmpRadius);
0201   blockTrackerLength.push_back(tmpLength);
0202 
0203   // PIXD1
0204   tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXD1CylinderRadius");
0205   tmpLength = p.getUntrackedParameter<std::vector<double>>("PXD1CylinderLength");
0206   trackerRadius.push_back(tmpRadius);
0207   trackerLength.push_back(tmpLength);
0208 
0209   // PIXD2
0210   tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXD2CylinderRadius");
0211   tmpLength = p.getUntrackedParameter<std::vector<double>>("PXD2CylinderLength");
0212   trackerRadius.push_back(tmpRadius);
0213   trackerLength.push_back(tmpLength);
0214 
0215   // PIXD Cables
0216   tmpRadius = p.getUntrackedParameter<std::vector<double>>("PXDCablesCylinderRadius");
0217   tmpLength = p.getUntrackedParameter<std::vector<double>>("PXDCablesCylinderLength");
0218   trackerRadius.push_back(tmpRadius);
0219   trackerLength.push_back(tmpLength);
0220 
0221   // All Pixel Disks
0222   blockTrackerRadius.push_back(tmpRadius);
0223   blockTrackerLength.push_back(tmpLength);
0224 
0225   // All Pixel
0226   subTrackerRadius.push_back(tmpRadius);
0227   subTrackerLength.push_back(tmpLength);
0228 
0229   // TIB1
0230   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIB1CylinderRadius");
0231   tmpLength = p.getUntrackedParameter<std::vector<double>>("TIB1CylinderLength");
0232   trackerRadius.push_back(tmpRadius);
0233   trackerLength.push_back(tmpLength);
0234 
0235   // TIB2
0236   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIB2CylinderRadius");
0237   tmpLength = p.getUntrackedParameter<std::vector<double>>("TIB2CylinderLength");
0238   trackerRadius.push_back(tmpRadius);
0239   trackerLength.push_back(tmpLength);
0240 
0241   // TIB3
0242   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIB3CylinderRadius");
0243   tmpLength = p.getUntrackedParameter<std::vector<double>>("TIB3CylinderLength");
0244   trackerRadius.push_back(tmpRadius);
0245   trackerLength.push_back(tmpLength);
0246 
0247   // TIB4
0248   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIB4CylinderRadius");
0249   tmpLength = p.getUntrackedParameter<std::vector<double>>("TIB4CylinderLength");
0250   trackerRadius.push_back(tmpRadius);
0251   trackerLength.push_back(tmpLength);
0252 
0253   // TIB Cables
0254   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIBCablesCylinderRadius");
0255   tmpLength = p.getUntrackedParameter<std::vector<double>>("TIBCablesCylinderLength");
0256   trackerRadius.push_back(tmpRadius);
0257   trackerLength.push_back(tmpLength);
0258 
0259   // All TIB
0260   blockTrackerRadius.push_back(tmpRadius);
0261   blockTrackerLength.push_back(tmpLength);
0262 
0263   // TID1
0264   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TID1CylinderRadius");
0265   tmpLength = p.getUntrackedParameter<std::vector<double>>("TID1CylinderLength");
0266   trackerRadius.push_back(tmpRadius);
0267   trackerLength.push_back(tmpLength);
0268 
0269   // TID2
0270   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TID2CylinderRadius");
0271   tmpLength = p.getUntrackedParameter<std::vector<double>>("TID2CylinderLength");
0272   trackerRadius.push_back(tmpRadius);
0273   trackerLength.push_back(tmpLength);
0274 
0275   // TID3
0276   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TID3CylinderRadius");
0277   tmpLength = p.getUntrackedParameter<std::vector<double>>("TID3CylinderLength");
0278   trackerRadius.push_back(tmpRadius);
0279   trackerLength.push_back(tmpLength);
0280 
0281   // TID Cables
0282   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TIDCablesCylinderRadius");
0283   tmpLength = p.getUntrackedParameter<std::vector<double>>("TIDCablesCylinderLength");
0284   trackerRadius.push_back(tmpRadius);
0285   trackerLength.push_back(tmpLength);
0286 
0287   // All TID
0288   blockTrackerRadius.push_back(tmpRadius);
0289   blockTrackerLength.push_back(tmpLength);
0290 
0291   // All Inner Tracker
0292   subTrackerRadius.push_back(tmpRadius);
0293   subTrackerLength.push_back(tmpLength);
0294 
0295   // TOB1
0296   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB1CylinderRadius");
0297   tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB1CylinderLength");
0298   trackerRadius.push_back(tmpRadius);
0299   trackerLength.push_back(tmpLength);
0300 
0301   // TOB2
0302   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB2CylinderRadius");
0303   tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB2CylinderLength");
0304   trackerRadius.push_back(tmpRadius);
0305   trackerLength.push_back(tmpLength);
0306 
0307   // TOB3
0308   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB3CylinderRadius");
0309   tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB3CylinderLength");
0310   trackerRadius.push_back(tmpRadius);
0311   trackerLength.push_back(tmpLength);
0312 
0313   // TOB4
0314   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB4CylinderRadius");
0315   tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB4CylinderLength");
0316   trackerRadius.push_back(tmpRadius);
0317   trackerLength.push_back(tmpLength);
0318 
0319   // TOB5
0320   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB5CylinderRadius");
0321   tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB5CylinderLength");
0322   trackerRadius.push_back(tmpRadius);
0323   trackerLength.push_back(tmpLength);
0324 
0325   // TOB6
0326   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOB6CylinderRadius");
0327   tmpLength = p.getUntrackedParameter<std::vector<double>>("TOB6CylinderLength");
0328   trackerRadius.push_back(tmpRadius);
0329   trackerLength.push_back(tmpLength);
0330 
0331   // TOB Cables
0332   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TOBCablesCylinderRadius");
0333   tmpLength = p.getUntrackedParameter<std::vector<double>>("TOBCablesCylinderLength");
0334   trackerRadius.push_back(tmpRadius);
0335   trackerLength.push_back(tmpLength);
0336 
0337   // All TOB
0338   blockTrackerRadius.push_back(tmpRadius);
0339   blockTrackerLength.push_back(tmpLength);
0340 
0341   // TEC1
0342   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC1CylinderRadius");
0343   tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC1CylinderLength");
0344   trackerRadius.push_back(tmpRadius);
0345   trackerLength.push_back(tmpLength);
0346 
0347   // TEC2
0348   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC2CylinderRadius");
0349   tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC2CylinderLength");
0350   trackerRadius.push_back(tmpRadius);
0351   trackerLength.push_back(tmpLength);
0352 
0353   // TEC3
0354   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC3CylinderRadius");
0355   tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC3CylinderLength");
0356   trackerRadius.push_back(tmpRadius);
0357   trackerLength.push_back(tmpLength);
0358 
0359   // TEC4
0360   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC4CylinderRadius");
0361   tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC4CylinderLength");
0362   trackerRadius.push_back(tmpRadius);
0363   trackerLength.push_back(tmpLength);
0364 
0365   // TEC5
0366   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC5CylinderRadius");
0367   tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC5CylinderLength");
0368   trackerRadius.push_back(tmpRadius);
0369   trackerLength.push_back(tmpLength);
0370 
0371   // TEC6
0372   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC6CylinderRadius");
0373   tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC6CylinderLength");
0374   trackerRadius.push_back(tmpRadius);
0375   trackerLength.push_back(tmpLength);
0376 
0377   // TEC7
0378   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC7CylinderRadius");
0379   tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC7CylinderLength");
0380   trackerRadius.push_back(tmpRadius);
0381   trackerLength.push_back(tmpLength);
0382 
0383   // TEC8
0384   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC8CylinderRadius");
0385   tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC8CylinderLength");
0386   trackerRadius.push_back(tmpRadius);
0387   trackerLength.push_back(tmpLength);
0388 
0389   // TEC9
0390   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TEC9CylinderRadius");
0391   tmpLength = p.getUntrackedParameter<std::vector<double>>("TEC9CylinderLength");
0392   trackerRadius.push_back(tmpRadius);
0393   trackerLength.push_back(tmpLength);
0394 
0395   // All TEC
0396   blockTrackerRadius.push_back(tmpRadius);
0397   blockTrackerLength.push_back(tmpLength);
0398 
0399   // All Outer
0400   subTrackerRadius.push_back(tmpRadius);
0401   subTrackerLength.push_back(tmpLength);
0402 
0403   // Outer Cables
0404   tmpRadius = p.getUntrackedParameter<std::vector<double>>("TrackerCablesCylinderRadius");
0405   tmpLength = p.getUntrackedParameter<std::vector<double>>("TrackerCablesCylinderLength");
0406   trackerRadius.push_back(tmpRadius);
0407   trackerLength.push_back(tmpLength);
0408 
0409   // All TEC
0410   h200.push_back(htmp);
0411   blockTrackerRadius.push_back(tmpRadius);
0412   blockTrackerLength.push_back(tmpLength);
0413 
0414   // All
0415   h300.push_back(htmp);
0416   subTrackerRadius.push_back(tmpRadius);
0417   subTrackerLength.push_back(tmpLength);
0418 }
0419 
0420 void testNuclearInteractions::bookHistograms(DQMStore::IBooker& ibooker,
0421                                              edm::Run const& iRun,
0422                                              edm::EventSetup const& iSetup) {
0423   ibooker.setCurrentFolder("testNuclearInteractions");
0424 
0425   h0[0] = ibooker.book2D("radioFull", "Full Tracker radiography", 1000, 0., 320., 1000, 0., 150.);
0426   h0[1] = ibooker.book2D("radioFast", "Fast Tracker radiography", 1000, 0., 320., 1000, 0., 150.);
0427   h1[0] = ibooker.book1D("vertexFull", "Full Nb of Vertices", 20, -0.5, 19.5);
0428   h1[1] = ibooker.book1D("vertexFast", "Fast Nb of Vertices", 20, -0.5, 19.5);
0429   h2[0] = ibooker.book1D("daughterFull", "Full Nb of daughters", 20, -0.5, 19.5);
0430   h2[1] = ibooker.book1D("daughterFast", "Fast Nb of daughters", 20, -0.5, 19.5);
0431   h3[0] = ibooker.book1D("ecmFull", "Full centre-of-mass energy", 100, 0., 10.);
0432   h3[1] = ibooker.book1D("ecmFast", "Fast centre-of-mass energy", 100, 0., 10.);
0433   h4[0] = ibooker.book1D("FecmFull", "Full c.m. energy fraction", 100, 0., 2.);
0434   h4[1] = ibooker.book1D("FecmFast", "Fast c.m. energy fraction", 100, 0., 2.);
0435   h5[0] = ibooker.book1D("FmomFull", "Full momemtum", 100, 0., 10.);
0436   h5[1] = ibooker.book1D("FmomFast", "Fast momemtum", 100, 0., 10.);
0437   h6[0] = ibooker.book1D("DeltaEFull4", "Full DeltaE", 2000, -1., 4.);
0438   h6[1] = ibooker.book1D("DeltaEFast4", "Fast DetlaE", 2000, -1., 4.);
0439   h7[0] = ibooker.book1D("DeltaEFull3", "Full DeltaE 3 daugh", 2000, -1., 4.);
0440   h7[1] = ibooker.book1D("DeltaEFast3", "Fast DetlaE 3 daugh", 2000, -1., 4.);
0441   h8[0] = ibooker.book1D("DeltaMFull4", "Full DeltaE", 2000, -10., 40.);
0442   h8[1] = ibooker.book1D("DeltaMFast4", "Fast DetlaE", 2000, -10., 40.);
0443   h9[0] = ibooker.book1D("DeltaMFull3", "Full DeltaE 3 daugh", 2000, -10., 40.);
0444   h9[1] = ibooker.book1D("DeltaMFast3", "Fast DetlaE 3 daugh", 2000, -10., 40.);
0445   h10[0] = ibooker.book1D("EafterFull", "E(after)/E(before) full", 200, 0., 4.);
0446   h10[1] = ibooker.book1D("EafterFast", "E(after)/E(before) fast", 200, 0., 4.);
0447   /*
0448   h6[0] = ibooker.book2D("radioFullRem1", "Full Tracker radiography", 1000, 0.,320.,1000,0., 150. );
0449   h6[1] = ibooker.book2D("radioFastRem1", "Fast Tracker radiography", 1000, 0.,320.,1000,0., 150. );
0450   h7[0] = ibooker.book2D("radioFullRem2", "Full Tracker radiography", 1000, 0.,320.,1000,0., 150. );
0451   h7[1] = ibooker.book2D("radioFullRem2", "Fast Tracker radiography", 1000, 0.,320.,1000,0., 150. );
0452   h8[0] = ibooker.book2D("radioFullBP", "Full BP radiography", 1000, 0.,320.,1000,0., 150. );
0453   h8[1] = ibooker.book2D("radioFastBP", "Fast BP radiography", 1000, 0.,320.,1000,0., 150. );
0454   h9[0] = ibooker.book2D("radioFullPX", "Full PX radiography", 1000, 0.,320.,1000,0., 150. );
0455   h9[1] = ibooker.book2D("radioFastPX", "Fast PX radiography", 1000, 0.,320.,1000,0., 150. );
0456   h10[0] = ibooker.book2D("radioFullTI", "Full TI radiography", 1000, 0.,320.,1000,0., 150. );
0457   h10[1] = ibooker.book2D("radioFastTI", "Fast TI radiography", 1000, 0.,320.,1000,0., 150. );
0458   h11[0] = ibooker.book2D("radioFullTO", "Full TO radiography", 1000, 0.,320.,1000,0., 150. );
0459   h11[1] = ibooker.book2D("radioFastTO", "Fast TO radiography", 1000, 0.,320.,1000,0., 150. );
0460   h12[0] = ibooker.book2D("radioFullCA", "Full CA radiography", 1000, 0.,320.,1000,0., 150. );
0461   h12[1] = ibooker.book2D("radioFastCA", "Fast CA radiography", 1000, 0.,320.,1000,0., 150. );
0462   */
0463 
0464   totalCharge[0] = ibooker.book1D("ChargeFull", "Total Charge (full)", 19, -9.5, 9.5);
0465   totalCharge[1] = ibooker.book1D("ChargeFast", "Total Charge (fast)", 19, -9.5, 9.5);
0466 
0467   // Beam Pipe
0468   htmp[0] = ibooker.book1D("BeamPipeFull", "Full Beam Pipe", 120, 0., 3.);
0469   htmp[1] = ibooker.book1D("BeamPipeFast", "Fast Beam Pipe", 120, 0., 3.);
0470   h100.push_back(htmp);
0471 
0472   // Beam Pipe (cont'd)
0473   htmp[0] = ibooker.book1D("BPFull", "Full Beam Pipe", 120, 0., 3.);
0474   htmp[1] = ibooker.book1D("BPFast", "Fast Beam Pipe", 120, 0., 3.);
0475   h300.push_back(htmp);
0476 
0477   // PIXB1
0478   htmp[0] = ibooker.book1D("PXB1Full", "Full Pixel Barrel 1", 120, 0., 3.);
0479   htmp[1] = ibooker.book1D("PXB1Fast", "Fast Pixel Barrel 1", 120, 0., 3.);
0480   h100.push_back(htmp);
0481 
0482   // PIXB2
0483   htmp[0] = ibooker.book1D("PXB2Full", "Full Pixel Barrel 2", 120, 0., 3.);
0484   htmp[1] = ibooker.book1D("PXB2Fast", "Fast Pixel Barrel 2", 120, 0., 3.);
0485   h100.push_back(htmp);
0486 
0487   // PIXB3
0488   htmp[0] = ibooker.book1D("PXB3Full", "Full Pixel Barrel 3", 120, 0., 3.);
0489   htmp[1] = ibooker.book1D("PXB3Fast", "Fast Pixel Barrel 3", 120, 0., 3.);
0490   h100.push_back(htmp);
0491 
0492   // PIXB Cables
0493   htmp[0] = ibooker.book1D("PXBCFull", "Full Pixel Barrel Cables", 120, 0., 3.);
0494   htmp[1] = ibooker.book1D("PXBCFast", "Fast Pixel Barrel Cables", 120, 0., 3.);
0495   h100.push_back(htmp);
0496 
0497   // All Pixel Barrel
0498   htmp[0] = ibooker.book1D("PXBFull", "Full Pixel Barrel", 120, 0., 3.);
0499   htmp[1] = ibooker.book1D("PXBFast", "Fast Pixel Barrel", 120, 0., 3.);
0500   h200.push_back(htmp);
0501 
0502   // PIXD1
0503   htmp[0] = ibooker.book1D("PXD1Full", "Full Pixel Disk 1", 120, 0., 3.);
0504   htmp[1] = ibooker.book1D("PXD1Fast", "Fast Pixel Disk 1", 120, 0., 3.);
0505   h100.push_back(htmp);
0506 
0507   // PIXD2
0508   htmp[0] = ibooker.book1D("PXD2Full", "Full Pixel Disk 2", 120, 0., 3.);
0509   htmp[1] = ibooker.book1D("PXD2Fast", "Fast Pixel Disk 2", 120, 0., 3.);
0510   h100.push_back(htmp);
0511 
0512   // PIXD Cables
0513   htmp[0] = ibooker.book1D("PXDCFull", "Full Pixel Disk Cables", 120, 0., 3.);
0514   htmp[1] = ibooker.book1D("PXDCFast", "Fast Pixel Disk Cables", 120, 0., 3.);
0515   h100.push_back(htmp);
0516 
0517   // All Pixel Disks
0518   htmp[0] = ibooker.book1D("PXDFull", "Full Pixel Disk", 120, 0., 3.);
0519   htmp[1] = ibooker.book1D("PXDFast", "Fast Pixel Disk", 120, 0., 3.);
0520   h200.push_back(htmp);
0521 
0522   // All Pixel
0523   htmp[0] = ibooker.book1D("PixelFull", "Full Pixel", 120, 0., 3.);
0524   htmp[1] = ibooker.book1D("PixelFast", "Fast Pixel", 120, 0., 3.);
0525   h300.push_back(htmp);
0526 
0527   // TIB1
0528   htmp[0] = ibooker.book1D("TIB1Full", "Full Tracker Inner Barrel 1", 120, 0., 3.);
0529   htmp[1] = ibooker.book1D("TIB1Fast", "Fast Tracker Inner Barrel 1", 120, 0., 3.);
0530   h100.push_back(htmp);
0531 
0532   // TIB2
0533   htmp[0] = ibooker.book1D("TIB2Full", "Full Tracker Inner Barrel 2", 120, 0., 3.);
0534   htmp[1] = ibooker.book1D("TIB2Fast", "Fast Tracker Inner Barrel 2", 120, 0., 3.);
0535   h100.push_back(htmp);
0536 
0537   // TIB3
0538   htmp[0] = ibooker.book1D("TIB3Full", "Full Tracker Inner Barrel 3", 120, 0., 3.);
0539   htmp[1] = ibooker.book1D("TIB3Fast", "Fast Tracker Inner Barrel 3", 120, 0., 3.);
0540   h100.push_back(htmp);
0541 
0542   // TIB4
0543   htmp[0] = ibooker.book1D("TIB4Full", "Full Tracker Inner Barrel 4", 120, 0., 3.);
0544   htmp[1] = ibooker.book1D("TIB4Fast", "Fast Tracker Inner Barrel 4", 120, 0., 3.);
0545   h100.push_back(htmp);
0546 
0547   // TIB Cables
0548   htmp[0] = ibooker.book1D("TIBCFull", "Full Tracker Inner Barrel Cables", 120, 0., 3.);
0549   htmp[1] = ibooker.book1D("TIBCFast", "Fast Tracker Inner Barrel Cables", 120, 0., 3.);
0550   h100.push_back(htmp);
0551 
0552   // All TIB
0553   htmp[0] = ibooker.book1D("TIBFull", "Full Tracker Inner Barrel", 120, 0., 3.);
0554   htmp[1] = ibooker.book1D("TIBFast", "Fast Tracker Inner Barrel", 120, 0., 3.);
0555   h200.push_back(htmp);
0556 
0557   // TID1
0558   htmp[0] = ibooker.book1D("TID1Full", "Full Tracker Inner Disk 1", 120, 0., 3.);
0559   htmp[1] = ibooker.book1D("TID1Fast", "Fast Tracker Inner Disk 1", 120, 0., 3.);
0560   h100.push_back(htmp);
0561 
0562   // TID2
0563   htmp[0] = ibooker.book1D("TID2Full", "Full Tracker Inner Disk 2", 120, 0., 3.);
0564   htmp[1] = ibooker.book1D("TID2Fast", "Fast Tracker Inner Disk 2", 120, 0., 3.);
0565   h100.push_back(htmp);
0566 
0567   // TID3
0568   htmp[0] = ibooker.book1D("TID3Full", "Full Tracker Inner Disk 3", 120, 0., 3.);
0569   htmp[1] = ibooker.book1D("TID3Fast", "Fast Tracker Inner Disk 3", 120, 0., 3.);
0570   h100.push_back(htmp);
0571 
0572   // TID Cables
0573   htmp[0] = ibooker.book1D("TIDCFull", "Full Tracker Inner Disk Cables", 120, 0., 3.);
0574   htmp[1] = ibooker.book1D("TIDCFast", "Fast Tracker Inner Disk Cables", 120, 0., 3.);
0575   h100.push_back(htmp);
0576 
0577   // All TID
0578   htmp[0] = ibooker.book1D("TIDFull", "Full Tracker Inner Disk", 120, 0., 3.);
0579   htmp[1] = ibooker.book1D("TIDFast", "Fast Tracker Inner Disk", 120, 0., 3.);
0580   h200.push_back(htmp);
0581 
0582   // All Inner Tracker
0583   htmp[0] = ibooker.book1D("InnerFull", "Full Inner Tracker", 120, 0., 3.);
0584   htmp[1] = ibooker.book1D("InnerFast", "Fast Inner Tracker", 120, 0., 3.);
0585   h300.push_back(htmp);
0586 
0587   // TOB1
0588   htmp[0] = ibooker.book1D("TOB1Full", "Full Tracker Outer Barrel 1", 120, 0., 3.);
0589   htmp[1] = ibooker.book1D("TOB1Fast", "Fast Tracker Outer Barrel 1", 120, 0., 3.);
0590   h100.push_back(htmp);
0591 
0592   // TOB2
0593   htmp[0] = ibooker.book1D("TOB2Full", "Full Tracker Outer Barrel 2", 120, 0., 3.);
0594   htmp[1] = ibooker.book1D("TOB2Fast", "Fast Tracker Outer Barrel 2", 120, 0., 3.);
0595   h100.push_back(htmp);
0596 
0597   // TOB3
0598   htmp[0] = ibooker.book1D("TOB3Full", "Full Tracker Outer Barrel 3", 120, 0., 3.);
0599   htmp[1] = ibooker.book1D("TOB3Fast", "Fast Tracker Outer Barrel 3", 120, 0., 3.);
0600   h100.push_back(htmp);
0601 
0602   // TOB4
0603   htmp[0] = ibooker.book1D("TOB4Full", "Full Tracker Outer Barrel 4", 120, 0., 3.);
0604   htmp[1] = ibooker.book1D("TOB4Fast", "Fast Tracker Outer Barrel 4", 120, 0., 3.);
0605   h100.push_back(htmp);
0606 
0607   // TOB5
0608   htmp[0] = ibooker.book1D("TOB5Full", "Full Tracker Outer Barrel 5", 120, 0., 3.);
0609   htmp[1] = ibooker.book1D("TOB5Fast", "Fast Tracker Outer Barrel 5", 120, 0., 3.);
0610   h100.push_back(htmp);
0611 
0612   // TOB6
0613   htmp[0] = ibooker.book1D("TOB6Full", "Full Tracker Outer Barrel 6", 120, 0., 3.);
0614   htmp[1] = ibooker.book1D("TOB6Fast", "Fast Tracker Outer Barrel 6", 120, 0., 3.);
0615   h100.push_back(htmp);
0616 
0617   // TOB Cables
0618   htmp[0] = ibooker.book1D("TOBCFull", "Full Tracker Outer Barrel Cables", 120, 0., 3.);
0619   htmp[1] = ibooker.book1D("TOBCFast", "Fast Tracker Outer Barrel Cables", 120, 0., 3.);
0620   h100.push_back(htmp);
0621 
0622   // All TOB
0623   htmp[0] = ibooker.book1D("TOBFull", "Full Tracker Outer Barrel", 120, 0., 3.);
0624   htmp[1] = ibooker.book1D("TOBFast", "Fast Tracker Outer Barrel", 120, 0., 3.);
0625   h200.push_back(htmp);
0626 
0627   // TEC1
0628   htmp[0] = ibooker.book1D("TEC1Full", "Full Tracker EndCap 1", 120, 0., 3.);
0629   htmp[1] = ibooker.book1D("TEC1Fast", "Fast Tracker Endcap 1", 120, 0., 3.);
0630   h100.push_back(htmp);
0631 
0632   // TEC2
0633   htmp[0] = ibooker.book1D("TEC2Full", "Full Tracker EndCap 2", 120, 0., 3.);
0634   htmp[1] = ibooker.book1D("TEC2Fast", "Fast Tracker Endcap 2", 120, 0., 3.);
0635   h100.push_back(htmp);
0636 
0637   // TEC3
0638   htmp[0] = ibooker.book1D("TEC3Full", "Full Tracker EndCap 3", 120, 0., 3.);
0639   htmp[1] = ibooker.book1D("TEC3Fast", "Fast Tracker Endcap 3", 120, 0., 3.);
0640   h100.push_back(htmp);
0641 
0642   // TEC4
0643   htmp[0] = ibooker.book1D("TEC4Full", "Full Tracker EndCap 4", 120, 0., 3.);
0644   htmp[1] = ibooker.book1D("TEC4Fast", "Fast Tracker Endcap 4", 120, 0., 3.);
0645   h100.push_back(htmp);
0646 
0647   // TEC5
0648   htmp[0] = ibooker.book1D("TEC5Full", "Full Tracker EndCap 5", 120, 0., 3.);
0649   htmp[1] = ibooker.book1D("TEC5Fast", "Fast Tracker Endcap 5", 120, 0., 3.);
0650   h100.push_back(htmp);
0651 
0652   // TEC6
0653   htmp[0] = ibooker.book1D("TEC6Full", "Full Tracker EndCap 6", 120, 0., 3.);
0654   htmp[1] = ibooker.book1D("TEC6Fast", "Fast Tracker Endcap 6", 120, 0., 3.);
0655   h100.push_back(htmp);
0656 
0657   // TEC7
0658   htmp[0] = ibooker.book1D("TEC7Full", "Full Tracker EndCap 7", 120, 0., 3.);
0659   htmp[1] = ibooker.book1D("TEC7Fast", "Fast Tracker Endcap 7", 120, 0., 3.);
0660   h100.push_back(htmp);
0661 
0662   // TEC8
0663   htmp[0] = ibooker.book1D("TEC8Full", "Full Tracker EndCap 8", 120, 0., 3.);
0664   htmp[1] = ibooker.book1D("TEC8Fast", "Fast Tracker Endcap 8", 120, 0., 3.);
0665   h100.push_back(htmp);
0666 
0667   // TEC9
0668   htmp[0] = ibooker.book1D("TEC9Full", "Full Tracker EndCap 9", 120, 0., 3.);
0669   htmp[1] = ibooker.book1D("TEC9Fast", "Fast Tracker Endcap 9", 120, 0., 3.);
0670   h100.push_back(htmp);
0671 
0672   // All TEC
0673   htmp[0] = ibooker.book1D("TECFull", "Full Tracker EndCap", 120, 0., 3.);
0674   htmp[1] = ibooker.book1D("TECFast", "Fast Tracker EndCap", 120, 0., 3.);
0675   h200.push_back(htmp);
0676 
0677   // All Outer
0678   htmp[0] = ibooker.book1D("OuterFull", "Full Outer Tracker", 120, 0., 3.);
0679   htmp[1] = ibooker.book1D("OuterFast", "Fast Outer Tracker", 120, 0., 3.);
0680   h300.push_back(htmp);
0681 
0682   // Outer Cables
0683   htmp[0] = ibooker.book1D("TECCFull", "Full Tracker Outer Cables", 120, 0., 3.);
0684   htmp[1] = ibooker.book1D("TECCFast", "Fast Tracker Outer Cables", 120, 0., 3.);
0685   h100.push_back(htmp);
0686 
0687   // All TEC
0688   htmp[0] = ibooker.book1D("CablesFull", "Full Tracker Cables", 120, 0., 3.);
0689   htmp[1] = ibooker.book1D("CablesFast", "Fast Tracker Cables", 120, 0., 3.);
0690   h200.push_back(htmp);
0691 
0692   // All
0693   htmp[0] = ibooker.book1D("TrackerFull", "Full Tracker", 120, 0., 3.);
0694   htmp[1] = ibooker.book1D("TrackerFast", "Fast Tracker", 120, 0., 3.);
0695   h300.push_back(htmp);
0696 }
0697 
0698 testNuclearInteractions::~testNuclearInteractions() {
0699   if (saveNU) {
0700     outFile->cd();
0701     // Fill the last (incomplete) nuEvent
0702     nuTree->Fill();
0703     // Conclude the writing on disk
0704     nuTree->Write();
0705     // Print information
0706     nuTree->Print();
0707     // And tidy up everything!
0708     //  outFile->Close();
0709     delete nuEvent;
0710     delete nuTree;
0711     delete outFile;
0712   }
0713 
0714   //  delete mySimEvent;
0715 }
0716 
0717 void testNuclearInteractions::dqmBeginRun(edm::Run const&, const edm::EventSetup& es) {
0718   // init Particle data table (from Pythia)
0719   const edm::ESHandle<HepPDT::ParticleDataTable> pdt = es.getHandle(tok_pdt_);
0720 
0721   mySimEvent[0]->initializePdt(&(*pdt));
0722   mySimEvent[1]->initializePdt(&(*pdt));
0723 }
0724 
0725 void testNuclearInteractions::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0726   ++totalNEvt;
0727   if (totalNEvt / 1000 * 1000 == totalNEvt)
0728     std::cout << "Number of event analysed/NU " << totalNEvt << " / " << totalNU << std::endl;
0729 
0730   std::unique_ptr<edm::SimTrackContainer> nuclSimTracks(new edm::SimTrackContainer);
0731 
0732   //  std::cout << "Fill full event " << std::endl;
0733   const edm::Handle<std::vector<SimTrack>>& fullSimTracks = iEvent.getHandle(tok_fullSimTk_);
0734   const edm::Handle<std::vector<SimVertex>>& fullSimVertices = iEvent.getHandle(tok_fullSimVx_);
0735   mySimEvent[0]->fill(*fullSimTracks, *fullSimVertices);
0736 
0737   //  std::cout << "Fill fast event " << std::endl;
0738   /* */
0739   //  if ( !saveNU ) {
0740   const edm::Handle<std::vector<SimTrack>>& fastSimTracks = iEvent.getHandle(tok_fastSimTk_);
0741   const edm::Handle<std::vector<SimVertex>>& fastSimVertices = iEvent.getHandle(tok_fastSimVx_);
0742   mySimEvent[1]->fill(*fastSimTracks, *fastSimVertices);
0743   //}
0744   /* */
0745 
0746   //mySimEvent[0]->print();
0747   XYZTLorentzVector theProtonMomentum(0., 0., 0., 0.939);
0748 
0749   // Save the object number count for a new NUevent
0750   if (saveNU) {
0751     if (ObjectNumber == -1 || nuEvent->nInteractions() == 1000) {
0752       ObjectNumber = TProcessID::GetObjectCount();
0753       nuEvent->reset();
0754     }
0755   }
0756 
0757   for (unsigned ievt = 0; ievt < 2; ++ievt) {
0758     //    std::cout << "Event number " << ievt << std::endl;
0759     //    mySimEvent[ievt]->print();
0760 
0761     //    const std::vector<FSimVertex>& fsimVertices = *(mySimEvent[ievt]->vertices() );
0762     //    if ( !fsimVertices.size() ) continue;
0763     if (!mySimEvent[ievt]->nVertices())
0764       continue;
0765     const FSimTrack& thePion = mySimEvent[ievt]->track(0);
0766 
0767     //    h1[ievt]->Fill(fsimVertices.size());
0768     //    if ( fsimVertices.size() == 1 ) continue;
0769     h1[ievt]->Fill(mySimEvent[ievt]->nVertices());
0770     // Count stopping particles
0771     if (mySimEvent[ievt]->nVertices() == 1) {
0772       if (thePion.trackerSurfaceMomentum().e() < 1E-10)
0773         ++stoppedPions[ievt];
0774     }
0775     if (mySimEvent[ievt]->nVertices() == 1)
0776       continue;
0777 
0778     const FSimVertex& thePionVertex = mySimEvent[ievt]->vertex(1);
0779 
0780     double zed = thePionVertex.position().z();
0781     double radius = thePionVertex.position().pt();
0782     double eta = thePionVertex.position().eta();
0783 
0784     h0[ievt]->Fill(fabs(thePionVertex.position().z()), thePionVertex.position().pt());
0785 
0786     // Pion's number of daughters
0787     // FSimTrack& thePion = mySimEvent[ievt]->track(0);
0788 
0789     unsigned ndaugh = thePionVertex.nDaughters();
0790     h2[ievt]->Fill(ndaugh);
0791 
0792     // Check for a second vertex
0793     //    bool theSecondVertex = fsimVertices.size() > 2 ?
0794     //      mySimEvent[ievt]->vertex(2).parent().id() == 0 : false;
0795     //    std::cout << "Plusieurs interactions ? " << theSecondVertex << std::endl;
0796 
0797     // First and last daughters
0798     int firstDaughter = -1;
0799     int lastDaughter = -1;
0800     if (thePionVertex.nDaughters()) {
0801       lastDaughter = thePionVertex.daughters()[thePionVertex.nDaughters() - 1];
0802       firstDaughter = thePionVertex.daughters()[0];
0803     }
0804     // Reject charged pion/kaon leptonic decays (already simulated in FAMOS)
0805     //    if ( thePionVertex.nDaughters() == 1 ) {
0806     //      const FSimTrack& myDaugh = mySimEvent[ievt]->track(firstDaughter);
0807     //      if (abs(myDaugh.type()) == 11 || abs(myDaugh.type()) == 13 ) return;
0808     //    }
0809 
0810     XYZTLorentzVector totMoth = thePion.momentum();
0811     XYZTLorentzVector totDaugh(0., 0., 0., 0.);
0812     // double qMoth = thePion.charge();
0813     // double qDaugh = 0;
0814     unsigned nleptons = 0;
0815     unsigned nothers = 0;
0816     if (!(firstDaughter < 0 || lastDaughter < 0)) {
0817       for (int idaugh = firstDaughter; idaugh <= lastDaughter; ++idaugh) {
0818         const FSimTrack& myDaugh = mySimEvent[ievt]->track(idaugh);
0819         totDaugh += myDaugh.momentum();
0820         //  qDaugh += myDaugh.charge();
0821         // Count the leptons
0822         if (abs(myDaugh.type()) == 11 || abs(myDaugh.type()) == 13)
0823           ++nleptons;
0824         // Count the hadrons
0825         if (abs(myDaugh.type()) != 111 && abs(myDaugh.type()) != 211)
0826           ++nothers;
0827       }
0828     }
0829 
0830     // Reject decays (less than one/four daughters, for pions and kaons)
0831     if ((abs(thePion.type()) == 211 && ndaugh == 1) || (abs(thePion.type()) == 130 && ndaugh < 4) ||
0832         (abs(thePion.type()) == 321 && ndaugh < 4)) {
0833       double diffE = (totMoth - totDaugh).E();
0834       double diffP = std::sqrt((totMoth - totDaugh).Vect().Mag2());
0835       double diffm = totMoth.M2() - totDaugh.M2();
0836       // double diffQ = qMoth-qDaugh;
0837       // Neutral particles (K0L) don't experience dE/dx nor multiple scattering nor B deflection.
0838       // E,p are conserved!
0839       if (abs(thePion.type()) == 130 && fabs(diffE) < 1E-5 && diffP < 1E-5)
0840         return;
0841       // Low-multiplicity final states with one electron or one muon
0842       // usually don't come from an interaction. All pions are taken care
0843       // of by this cut
0844       if (nleptons == 1)
0845         return;
0846       // Reserve of tricks in case it does not work
0847       /*
0848       BaseParticlePropagator pDaugh(totDaugh,thePion.endVertex().position(),qDaugh);
0849       double d0 = pDaugh.xyImpactParameter();
0850       double z0 = pDaugh.zImpactParameter();
0851       pDaugh.propagateToNominalVertex();
0852       diffP = std::sqrt((totMoth-pDaugh.Momentum()).Vect().Mag2());
0853       */
0854       // Charge kaons may experience dE/dx and multiple scattering -> relax the cuts
0855       h7[ievt]->Fill(diffE);
0856       h9[ievt]->Fill(diffm);
0857       if (abs(thePion.type()) != 211 &&    // Ignore pions
0858           diffE > -1E-5 && diffE < 0.1 &&  // Small deltaE - to be checked as f(E)
0859           nothers == 0)
0860         return;  // Only pions in the decays
0861       h6[ievt]->Fill(diffE);
0862       h8[ievt]->Fill(diffm);
0863     }
0864 
0865     if (ndaugh)
0866       ++interactingPions[ievt];
0867     else
0868       ++stoppedPions[ievt];
0869 
0870     // Find the daughters, and boost them.
0871     if (!(firstDaughter < 0 || lastDaughter < 0)) {
0872       // Compute the boost for the cm frame, and the cm energy.
0873       XYZTLorentzVector theBoost = thePion.momentum() + theProtonMomentum;
0874       double ecm = theBoost.mag();
0875       theBoost /= theBoost.e();
0876       if (ievt == 0 && saveNU && totalNEvt < maxNU) {
0877         NUEvent::NUInteraction interaction;
0878         interaction.first = nuEvent->nParticles();
0879         interaction.last = interaction.first + lastDaughter - firstDaughter;
0880         nuEvent->addNUInteraction(interaction);
0881         ++totalNU;
0882       }
0883 
0884       // A few checks
0885       double qTot = 0.;
0886       double eBefore = thePion.momentum().E();
0887       double eAfter = 0.;
0888       XYZTLorentzVector theTotal(0., 0., 0., 0.);
0889 
0890       // Rotation to bring the collision axis around z
0891       XYZVector zAxis(0., 0., 1.);
0892       XYZVector rotationAxis = (theBoost.Vect().Cross(zAxis)).Unit();
0893       double rotationAngle = std::acos(theBoost.Vect().Unit().Z());
0894       RawParticle::Rotation rotation(rotationAxis, rotationAngle);
0895 
0896       for (int idaugh = firstDaughter; idaugh <= lastDaughter; ++idaugh) {
0897         // The track
0898         const FSimTrack& myDaugh = mySimEvent[ievt]->track(idaugh);
0899         qTot += myDaugh.charge();
0900         //  std::cout << "Daughter " << idaugh << " " << myDaugh << std::endl;
0901         RawParticle theMom(myDaugh.momentum());
0902         eAfter += theMom.E();
0903 
0904         // Boost the track
0905         theMom.boost(-theBoost.x(), -theBoost.y(), -theBoost.z());
0906         theTotal = theTotal + theMom.momentum();
0907 
0908         // Rotate ->  along the Z axis
0909         theMom.rotate(rotation);
0910 
0911         // Save the fully simulated tracks
0912         if (ievt == 0 && saveNU && totalNEvt <= maxNU) {
0913           NUEvent::NUParticle particle;
0914           particle.px = theMom.px() / ecm;
0915           particle.py = theMom.py() / ecm;
0916           particle.pz = theMom.pz() / ecm;
0917           particle.mass = theMom.mag();
0918           particle.id = myDaugh.type();
0919           nuEvent->addNUParticle(particle);
0920           //      SimTrack nuclSimTrack(myDaugh.type(),theMom/ecm,-1,-1);
0921           //      nuclSimTracks->push_back(nuclSimTrack);
0922         }
0923       }
0924 
0925       // Save some histograms
0926       h3[ievt]->Fill(ecm);
0927       h4[ievt]->Fill(theTotal.mag() / ecm);
0928       h5[ievt]->Fill(sqrt(theTotal.Vect().mag2()));
0929       h10[ievt]->Fill(eAfter / eBefore);
0930 
0931       // Total charge of daughters
0932       totalCharge[ievt]->Fill(qTot);
0933 
0934       // Fill the individual layer histograms !
0935       bool filled = false;
0936       for (unsigned hist = 0; hist < h100.size() && !filled; ++hist) {
0937         if (radius < trackerRadius[hist][ievt] && fabs(zed) < trackerLength[hist][ievt]) {
0938           h100[hist][ievt]->Fill(eta);
0939           filled = true;
0940         }
0941       }
0942 
0943       // Fill the block histograms !
0944       filled = false;
0945       for (unsigned hist = 0; hist < h200.size() && !filled; ++hist) {
0946         if (radius < blockTrackerRadius[hist][ievt] && fabs(zed) < blockTrackerLength[hist][ievt]) {
0947           h200[hist][ievt]->Fill(eta);
0948           filled = true;
0949         }
0950       }
0951 
0952       // Fill the cumulative histograms !
0953       for (unsigned hist = 0; hist < h300.size(); ++hist) {
0954         if ((radius < subTrackerRadius[hist][ievt] && fabs(zed) < subTrackerLength[hist][ievt]) ||
0955             (hist == 2 && radius < subTrackerRadius[1][ievt] && fabs(zed) < subTrackerLength[1][ievt])) {
0956           h300[hist][ievt]->Fill(eta);
0957         }
0958       }
0959     }
0960 
0961     // Save the fully simulated tracks from the nuclear interaction
0962     if (ievt == 0 && saveNU && totalNEvt <= maxNU) {
0963       if (nuEvent->nInteractions() == 1000) {
0964         // Reset Event object count to avoid memory overflows
0965         TProcessID::SetObjectCount(ObjectNumber);
0966         // Save the nuEvent
0967         std::cout << "Saved " << nuEvent->nInteractions() << " Interaction(s) with " << nuEvent->nParticles()
0968                   << " Particles in the NUEvent " << std::endl;
0969         outFile->cd();
0970         nuTree->Fill();
0971         //  nuTree->Print();
0972       }
0973     }
0974   }
0975 }
0976 
0977 //define this as a plug-in
0978 
0979 DEFINE_FWK_MODULE(testNuclearInteractions);