52 produces<ValueMap<DeDxData> >();
59 if(!usePixel && !useStrip)
60 edm::LogWarning(
"DeDxHitsProducer") <<
"Pixel Hits AND Strip Hits will not be used to estimate dEdx --> BUG, Please Update the config file";
75 useCalibration = iConfig.
getParameter<
bool>(
"UseCalibration");
76 m_calibrationPath = iConfig.
getParameter<
string>(
"calibrationPath");
78 Prob_ChargePath =
NULL;
88 if( strcmp(Reccord.c_str(),
"SiStripDeDxMip_3D_Rcd")==0){
90 }
else if(strcmp(Reccord.c_str(),
"SiStripDeDxPion_3D_Rcd")==0){
92 }
else if(strcmp(Reccord.c_str(),
"SiStripDeDxKaon_3D_Rcd")==0){
94 }
else if(strcmp(Reccord.c_str(),
"SiStripDeDxProton_3D_Rcd")==0){
96 }
else if(strcmp(Reccord.c_str(),
"SiStripDeDxElectron_3D_Rcd")==0){
103 DeDxMap_ = *DeDxMapHandle_.
product();
105 double xmin = DeDxMap_.rangeX().min;
106 double xmax = DeDxMap_.rangeX().max;
107 double ymin = DeDxMap_.rangeY().min;
108 double ymax = DeDxMap_.rangeY().max;
109 double zmin = DeDxMap_.rangeZ().min;
110 double zmax = DeDxMap_.rangeZ().max;
112 if(Prob_ChargePath)
delete Prob_ChargePath;
113 Prob_ChargePath =
new TH3D (
"Prob_ChargePath" ,
"Prob_ChargePath" , DeDxMap_.numberOfBinsX(),
xmin,
xmax, DeDxMap_.numberOfBinsY() ,
ymin,
ymax, DeDxMap_.numberOfBinsZ(),
zmin,
zmax);
117 if(strcmp(ProbabilityMode.c_str(),
"Accumulation")==0){
119 for(
int i=0;
i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;
i++){
121 for(
int j=0;
j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;
j++){
125 for(
int k=0;
k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;
k++){ Ni+=DeDxMap_.binContent(
i,
j,
k);}
127 for(
int k=0;
k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;
k++){
129 for(
int l=0;
l<=
k;
l++){ tmp+=DeDxMap_.binContent(
i,
j,
l);}
132 Prob_ChargePath->SetBinContent (
i,
j,
k, tmp/Ni);
135 Prob_ChargePath->SetBinContent (
i,
j,
k, 0);
140 }
else if(strcmp(ProbabilityMode.c_str(),
"Integral")==0){
142 for(
int i=0;
i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;
i++){
144 for(
int j=0;
j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;
j++){
148 for(
int k=0;
k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;
k++){ Ni+=DeDxMap_.binContent(
i,
j,
k);}
150 for(
int k=0;
k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;
k++){
151 double tmp = DeDxMap_.binContent(
i,
j,
k);
154 Prob_ChargePath->SetBinContent (
i,
j,
k, tmp/Ni);
157 Prob_ChargePath->SetBinContent (
i,
j,
k, 0);
182 if(MODsColl.size()==0){
187 vector<GeomDet*>
Det = tkGeom->dets();
188 for(
unsigned int i=0;
i<Det.size();
i++){
189 DetId Detid = Det[
i]->geographicalId();
196 if(!DetUnit)
continue;
199 unsigned int NAPV = Topo.
nstrips()/128;
212 MODsColl[MOD->
DetId] = MOD;
216 MakeCalibrationMap();
234 iEvent.
getByLabel(m_trajTrackAssociationTag, trajTrackAssociationHandle);
238 iEvent.
getByLabel(m_tracksTag,trackCollectionHandle);
244 std::vector<DeDxData> dEdxDiscrims( TrajToTrackMap.
size() );
246 unsigned track_index = 0;
248 dEdxDiscrims[track_index] =
DeDxData(-1, -2, 0 );
250 const Track track = *it->val;
253 if(track.
eta() <MinTrackEta || track.
eta()>MaxTrackEta ){
continue;}
254 if(track.
p() <MinTrackMomentum || track.
p() >MaxTrackMomentum){
continue;}
255 if(track.
found()<MinTrackHits ){
continue;}
257 std::vector<double> vect_probs;
258 vector<TrajectoryMeasurement> measurements = traj.
measurements();
260 unsigned int NClusterSaturating = 0;
261 for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
264 if( !trajState.
isValid() )
continue;
267 const SiStripRecHit2D* sistripsimplehit =
dynamic_cast<const SiStripRecHit2D*
>(hit);
268 const SiStripMatchedRecHit2D* sistripmatchedhit =
dynamic_cast<const SiStripMatchedRecHit2D*
>(hit);
269 const SiStripRecHit1D* sistripsimple1dhit =
dynamic_cast<const SiStripRecHit1D*
>(hit);
272 if(sistripsimplehit){
273 Prob = GetProbability(
DeDxTools::GetCluster(sistripsimplehit), trajState,sistripsimplehit->geographicalId());
275 if(Prob>=0) vect_probs.push_back(Prob);
279 }
else if(sistripmatchedhit){
280 Prob = GetProbability(
DeDxTools::GetCluster(sistripmatchedhit->monoHit()), trajState, sistripmatchedhit->monoId());
282 if(Prob>=0) vect_probs.push_back(Prob);
284 Prob = GetProbability(
DeDxTools::GetCluster(sistripmatchedhit->stereoHit()), trajState,sistripmatchedhit->stereoId());
285 if(Prob>=0) vect_probs.push_back(Prob);
287 if(ClusterSaturatingStrip(
DeDxTools::GetCluster(sistripmatchedhit->monoHit()),sistripmatchedhit->monoId()) >0)NClusterSaturating++;
288 if(ClusterSaturatingStrip(
DeDxTools::GetCluster(sistripmatchedhit->stereoHit()),sistripmatchedhit->stereoId())>0)NClusterSaturating++;
289 }
else if(sistripsimple1dhit){
290 Prob = GetProbability(
DeDxTools::GetCluster(sistripsimple1dhit), trajState, sistripsimple1dhit->geographicalId());
292 if(Prob>=0) vect_probs.push_back(Prob);
298 double estimator = ComputeDiscriminator(vect_probs);
299 int size = vect_probs.size();
304 Error = NClusterSaturating;
305 dEdxDiscrims[track_index] =
DeDxData(estimator, Error, size );
310 filler.
insert(trackCollectionHandle, dEdxDiscrims.begin(), dEdxDiscrims.end());
312 iEvent.
put(trackDeDxDiscrimAssociation);
317 const uint32_t & detId){
318 const vector<uint8_t>& ampls = cluster->
amplitudes();
321 if(useCalibration)MOD = MODsColl[detId];
323 int SaturatingStrip = 0;
324 for(
unsigned int i=0;
i<ampls.size();
i++){
325 int StripCharge = ampls[
i];
326 if(MOD){StripCharge = (int)(StripCharge / MOD->
Gain);}
327 if(StripCharge>=254)SaturatingStrip++;
329 return SaturatingStrip;
333 const uint32_t & detId)
337 double cosine = trackDirection.
z()/trackDirection.
mag();
338 const vector<uint8_t>& ampls = cluster->
amplitudes();
345 if( ampls.size()>MaxNrStrips) {
return -1;}
354 for(
unsigned int i=0;
i<ampls.size();
i++){
355 int CalibratedCharge = ampls[
i];
356 CalibratedCharge = (int)(CalibratedCharge / MOD->
Gain);
357 if(CalibratedCharge>=1024){
358 CalibratedCharge = 255;
359 }
else if(CalibratedCharge>254){
360 CalibratedCharge = 254;
362 charge+=CalibratedCharge;
370 int BinY = Prob_ChargePath->GetYaxis()->FindBin(path);
371 int BinZ = Prob_ChargePath->GetZaxis()->FindBin(charge/path);
372 return Prob_ChargePath->GetBinContent(BinX,BinY,BinZ);
378 double estimator = -1;
379 int size = vect_probs.size();
380 if(size<=0)
return estimator;
385 if(vect_probs[
i]<=0.0001){P *=
pow(0.0001 , 1.0/size);}
386 else {P *=
pow(vect_probs[
i], 1.0/size);}
389 }
else if(Formula==1){
390 std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
391 for(
int i=0;
i<
size;
i++){
if(vect_probs[
i]<=0.0001)vect_probs[
i] = 0.0001; }
394 for(
int i=0;
i<
size;
i++){ SumJet+=
log(vect_probs[
i]); }
396 double Loginvlog=
log(-SumJet);
402 Prob+=
exp(
l*Loginvlog-
log(1.*lfact));
405 double LogProb=
log(Prob);
407 estimator = -log10(ProbJet)/4.;
408 estimator = 1-estimator;
409 }
else if(Formula==2){
410 std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
411 double P = 1.0/(12*
size);
413 P +=
pow(vect_probs[
i-1] - ((2.0*
i-1.0)/(2.0*size)),2);
418 std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
419 double P = 1.0/(12*
size);
421 P += vect_probs[
i-1] *
pow(vect_probs[
i-1] - ((2.0*
i-1.0)/(2.0*size)),2);
432 if(!useCalibration)
return;
435 TChain* t1 =
new TChain(
"SiStripCalib/APVGain");
436 t1->Add(m_calibrationPath.c_str());
438 unsigned int tree_DetId;
439 unsigned char tree_APVId;
442 t1->SetBranchAddress(
"DetId" ,&tree_DetId );
443 t1->SetBranchAddress(
"APVId" ,&tree_APVId );
444 t1->SetBranchAddress(
"Gain" ,&tree_Gain );
446 for (
unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
447 t1->GetEntry(ientry);
449 MOD->
Gain = tree_Gain;
double p() const
momentum vector magnitude
T getParameter(std::string const &) const
virtual int nstrips() const =0
T getUntrackedParameter(std::string const &, T const &) const
virtual void produce(edm::Event &, const edm::EventSetup &) override
~DeDxDiscriminatorProducer()
const_iterator end() const
last iterator over the map (read only)
LocalVector localDirection() const
double GetProbability(const SiStripCluster *cluster, TrajectoryStateOnSurface trajState, const uint32_t &)
void insert(const H &h, I begin, I end)
DeDxDiscriminatorProducer(const edm::ParameterSet &)
const Bounds & bounds() const
virtual void beginRun(edm::Run const &run, const edm::EventSetup &) override
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
const Plane & surface() const
The nominal surface of the GeomDet.
uint32_t geographicalId() const
uint32_t rawId() const
get the raw id
LocalVector localMomentum() const
virtual float thickness() const =0
DataContainer const & measurements() const
const Surface::PositionType & position() const
The position (origin of the R.F.)
double eta() const
pseudorapidity of momentum vector
const T & max(const T &a, const T &b)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
void MakeCalibrationMap()
DEFINE_FWK_MODULE(CosmicTrackingParticleSelector)
size_type size() const
map size
T const * product() const
int ClusterSaturatingStrip(const SiStripCluster *cluster, const uint32_t &)
T const * product() const
std::vector< std::vector< double > > tmp
unsigned short found() const
Number of valid hits on track.
T transverse() const
Another name for perp()
double ComputeDiscriminator(std::vector< double > &vect_probs)
const_iterator begin() const
first iterator over the map (read only)
const BasicVectorType & basicVector() const
tuple size
Write out results.
const std::vector< uint8_t > & amplitudes() const
Power< A, B >::type pow(const A &a, const B &b)