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";
74 useCalibration = iConfig.
getParameter<
bool>(
"UseCalibration");
75 m_calibrationPath = iConfig.
getParameter<
string>(
"calibrationPath");
77 Prob_ChargePath =
NULL;
87 if( strcmp(Reccord.c_str(),
"SiStripDeDxMip_3D_Rcd")==0){
89 }
else if(strcmp(Reccord.c_str(),
"SiStripDeDxPion_3D_Rcd")==0){
91 }
else if(strcmp(Reccord.c_str(),
"SiStripDeDxKaon_3D_Rcd")==0){
93 }
else if(strcmp(Reccord.c_str(),
"SiStripDeDxProton_3D_Rcd")==0){
95 }
else if(strcmp(Reccord.c_str(),
"SiStripDeDxElectron_3D_Rcd")==0){
102 DeDxMap_ = *DeDxMapHandle_.
product();
104 double xmin = DeDxMap_.rangeX().min;
105 double xmax = DeDxMap_.rangeX().max;
106 double ymin = DeDxMap_.rangeY().min;
107 double ymax = DeDxMap_.rangeY().max;
108 double zmin = DeDxMap_.rangeZ().min;
109 double zmax = DeDxMap_.rangeZ().max;
111 if(Prob_ChargePath)
delete Prob_ChargePath;
112 Prob_ChargePath =
new TH3D (
"Prob_ChargePath" ,
"Prob_ChargePath" , DeDxMap_.numberOfBinsX(), xmin, xmax, DeDxMap_.numberOfBinsY() , ymin, ymax, DeDxMap_.numberOfBinsZ(), zmin, zmax);
116 if(strcmp(ProbabilityMode.c_str(),
"Accumulation")==0){
118 for(
int i=0;
i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;
i++){
120 for(
int j=0;
j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;
j++){
124 for(
int k=0;
k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;
k++){ Ni+=DeDxMap_.binContent(
i,
j,
k);}
126 for(
int k=0;
k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;
k++){
128 for(
int l=0;
l<=
k;
l++){ tmp+=DeDxMap_.binContent(
i,
j,
l);}
131 Prob_ChargePath->SetBinContent (
i,
j,
k, tmp/Ni);
134 Prob_ChargePath->SetBinContent (
i,
j,
k, 0);
139 }
else if(strcmp(ProbabilityMode.c_str(),
"Integral")==0){
141 for(
int i=0;
i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;
i++){
143 for(
int j=0;
j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;
j++){
147 for(
int k=0;
k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;
k++){ Ni+=DeDxMap_.binContent(
i,
j,
k);}
149 for(
int k=0;
k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;
k++){
150 double tmp = DeDxMap_.binContent(
i,
j,
k);
153 Prob_ChargePath->SetBinContent (
i,
j,
k, tmp/Ni);
156 Prob_ChargePath->SetBinContent (
i,
j,
k, 0);
181 if(MODsColl.size()==0){
186 vector<GeomDet*>
Det = tkGeom->dets();
187 for(
unsigned int i=0;
i<Det.size();
i++){
188 DetId Detid = Det[
i]->geographicalId();
195 if(!DetUnit)
continue;
198 unsigned int NAPV = Topo.
nstrips()/128;
211 MODsColl[MOD->
DetId] = MOD;
215 MakeCalibrationMap();
240 iEvent.
getByLabel(m_trajTrackAssociationTag, trajTrackAssociationHandle);
244 iEvent.
getByLabel(m_tracksTag,trackCollectionHandle);
250 std::vector<DeDxData> dEdxDiscrims( TrajToTrackMap.
size() );
252 unsigned track_index = 0;
254 dEdxDiscrims[track_index] =
DeDxData(-1, -2, 0 );
260 if(track.
p() <MinTrackMomentum || track.
p() >MaxTrackMomentum){
continue;}
261 if(track.
found()<MinTrackHits ){
continue;}
263 std::vector<double> vect_probs;
264 vector<TrajectoryMeasurement> measurements = traj.
measurements();
266 unsigned int NClusterSaturating = 0;
267 for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
270 if( !trajState.
isValid() )
continue;
280 Prob = GetProbability((sistripsimplehit->
cluster()).
get(), trajState);
if(Prob>=0) vect_probs.push_back(Prob);
281 if(ClusterSaturatingStrip((sistripsimplehit->
cluster()).
get())>0)NClusterSaturating++;
282 }
else if(sistripmatchedhit){
283 Prob = GetProbability((sistripmatchedhit->
monoHit()->
cluster()).
get(), trajState);
if(Prob>=0) vect_probs.push_back(Prob);
284 Prob = GetProbability((sistripmatchedhit->
stereoHit()->
cluster()).
get(), trajState);
if(Prob>=0) vect_probs.push_back(Prob);
285 if(ClusterSaturatingStrip((sistripmatchedhit->
monoHit()->
cluster()).
get()) >0)NClusterSaturating++;
286 if(ClusterSaturatingStrip((sistripmatchedhit->
stereoHit()->
cluster()).
get())>0)NClusterSaturating++;
287 }
else if(sistripsimple1dhit){
288 Prob = GetProbability((sistripsimple1dhit->
cluster()).
get(), trajState);
if(Prob>=0) vect_probs.push_back(Prob);
289 if(ClusterSaturatingStrip((sistripsimple1dhit->
cluster()).
get())>0)NClusterSaturating++;
294 double estimator = ComputeDiscriminator(vect_probs);
295 int size = vect_probs.size();
300 Error = NClusterSaturating;
301 dEdxDiscrims[track_index] =
DeDxData(estimator, Error, size );
306 filler.
insert(trackCollectionHandle, dEdxDiscrims.begin(), dEdxDiscrims.end());
308 iEvent.
put(trackDeDxDiscrimAssociation);
313 const vector<uint8_t>& ampls = cluster->
amplitudes();
318 int SaturatingStrip = 0;
319 for(
unsigned int i=0;
i<ampls.size();
i++){
320 int StripCharge = ampls[
i];
321 if(MOD){StripCharge = (int)(StripCharge / MOD->
Gain);}
322 if(StripCharge>=254)SaturatingStrip++;
324 return SaturatingStrip;
331 double cosine = trackDirection.
z()/trackDirection.
mag();
332 const vector<uint8_t>& ampls = cluster->
amplitudes();
339 if( ampls.size()>MaxNrStrips) {
return -1;}
348 for(
unsigned int i=0;
i<ampls.size();
i++){
349 int CalibratedCharge = ampls[
i];
350 CalibratedCharge = (int)(CalibratedCharge / MOD->
Gain);
351 if(CalibratedCharge>=1024){
352 CalibratedCharge = 255;
353 }
else if(CalibratedCharge>254){
354 CalibratedCharge = 254;
356 charge+=CalibratedCharge;
364 int BinY = Prob_ChargePath->GetYaxis()->FindBin(path);
365 int BinZ = Prob_ChargePath->GetZaxis()->FindBin(charge/path);
367 return Prob_ChargePath->GetBinContent(BinX,BinY,BinZ);
374 int size = vect_probs.size();
380 if(vect_probs[
i]<=0.0001){P *=
pow(0.0001 , 1.0/size);}
381 else {P *=
pow(vect_probs[
i], 1.0/size);}
384 }
else if(Formula==1){
385 std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
386 for(
int i=0;
i<
size;
i++){
if(vect_probs[
i]<=0.0001)vect_probs[
i] = 0.0001; }
389 for(
int i=0;
i<
size;
i++){ SumJet+=
log(vect_probs[
i]); }
391 double Loginvlog=
log(-SumJet);
397 Prob+=
exp(
l*Loginvlog-
log(1.*lfact));
400 double LogProb=
log(Prob);
402 estimator = -log10(ProbJet)/4.;
404 }
else if(Formula==2){
405 std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
406 double P = 1.0/(12*
size);
408 P +=
pow(vect_probs[
i-1] - ((2.0*
i-1.0)/(2.0*size)),2);
413 std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
414 double P = 1.0/(12*
size);
416 P += vect_probs[
i-1] *
pow(vect_probs[
i-1] - ((2.0*
i-1.0)/(2.0*size)),2);
427 if(!useCalibration)
return;
430 TChain* t1 =
new TChain(
"SiStripCalib/APVGain");
431 t1->Add(m_calibrationPath.c_str());
433 unsigned int tree_DetId;
434 unsigned char tree_APVId;
437 t1->SetBranchAddress(
"DetId" ,&tree_DetId );
438 t1->SetBranchAddress(
"APVId" ,&tree_APVId );
439 t1->SetBranchAddress(
"Gain" ,&tree_Gain );
441 for (
unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
442 t1->GetEntry(ientry);
444 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
double GetProbability(const SiStripCluster *cluster, TrajectoryStateOnSurface trajState)
~DeDxDiscriminatorProducer()
const_iterator end() const
last iterator over the map (read only)
const SiStripRecHit2D * stereoHit() const
LocalVector localDirection() const
int ClusterSaturatingStrip(const SiStripCluster *cluster)
void insert(const H &h, I begin, I end)
DeDxDiscriminatorProducer(const edm::ParameterSet &)
Exp< T >::type exp(const T &t)
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
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.
virtual void beginRun(edm::Run &run, const edm::EventSetup &)
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)
const Bounds & bounds() const
Log< T >::type log(const T &t)
size_type size() const
map size
virtual void produce(edm::Event &, const edm::EventSetup &)
ClusterRef const & cluster() const
T const * product() const
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)
ClusterRef const & cluster() const
const_iterator begin() const
first iterator over the map (read only)
const SiStripRecHit2D * monoHit() const
const BasicVectorType & basicVector() const
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
tuple size
Write out results.
const std::vector< uint8_t > & amplitudes() const
Power< A, B >::type pow(const A &a, const B &b)