Changeset 224 in PSPA for Interface_Web/trunk/pspaWT/src/particleBeam.cc

Timestamp:

Dec 21, 2012, 5:30:31 PM (12 years ago)

Author:

touze

Message:

embryon d'histo... à completer

File:

• Interface_Web/trunk/pspaWT/src/particleBeam.cc (modified) (1 diff)

Interface_Web/trunk/pspaWT/src/particleBeam.cc

@@ -519 +519 @@
     }
 }
+
+void particleBeam::histogramme(vector<double>&xcor,vector<int>& hist,int& cnts)
+{
+  double gammin= GRAND;
+  double gammax= -gammin;
+  double Emoy= 0.0;
+  double ecatyp= 0.0;
+
+  for (unsigned int k = 0; k < goodPartic_.size(); k++)
+    {
+      double gamma = goodPartic_.at(k).getGamma();
+      double EMev = (gamma-1.0)*ERESTMeV;
+      if (gamma < gammin) gammin = gamma;
+      else if (gamma > gammax) gammax = gamma;
+      Emoy += EMev;
+      ecatyp += EMev*EMev;
+    }
+
+  double sum= (float)goodPartic_.size();
+  Emoy /= sum;
+  ecatyp /= sum;
+  ecatyp = sqrt(abs(ecatyp-Emoy*Emoy));
+
+  double Emin = (gammin-1.0)*ERESTMeV;
+  double Emax = (gammax-1.0)*ERESTMeV;
+  cout << "energie cinetique -moyenne " << Emoy << " Mev " << "-mini " << Emin << " Mev " << "-maxi " << Emax << " Mev " << "ecart type " << ecatyp*1000.0 << " Kev" << endl;
+
+  vector<double> Eshf;
+  for (unsigned int k = 0; k < goodPartic_.size(); k++)
+    {
+      double gamma = goodPartic_.at(k).getGamma();
+      double EMev = (gamma-1.0)*ERESTMeV;
+      Eshf.push_back(EMev-Emoy);
+    }
+
+  //////////////////////////////////////////////////////////////////////////////
+  
+  // demi fenetre en energie, et pas de l'histogramme
+  double hfene= max(3.*ecatyp-Emoy,Emoy-3.*ecatyp);
+  double hpas = hfene/50.;
+
+  cout << "demi fenetre " << hfene << ", hpas= " << hpas << endl;
+
+  double vmin = -hfene;
+  double dfen = 2.*hfene;
+  int ihist = dfen/hpas;
+  double phist = ihist*hpas;
+  double dpas = hpas-(dfen-phist);
+  if(dpas <= hpas*1.e-03) {
+    ihist++;
+    phist= ihist*hpas;
+  }
+  double vmax= vmin+hpas*ihist;
+  
+  cout << "'xAxisNumberOfBins= " << ihist <<", xAxisMinimum= " << vmin << ", xAxisMaximum= " << vmax << ", NParticules= " << Eshf.size() << ", phist " << phist << endl;
+
+  xcor= vector<double>(ihist);
+  for (unsigned i = 0; i < ihist; ++i) {
+    xcor[i]= vmin+i*hpas;
+  }
+
+  hist= vector<int>(ihist,0);
+  for (unsigned i = 0; i < Eshf.size(); ++i) {
+    double var= Eshf[i]-vmin;
+    if(var < 0 || var >= phist) cout<<"out of range "<<var<<", ("<< i<<")"<< endl;
+    int k= var/hpas;
+    int kk= (int)floor(var/hpas);
+    if(i%20 == 0) cout<<"v("<<i<<")= " <<var<<" ["<<k<<"-"<<kk<<"], "<<endl; 
+    hist[kk]++;
+  }
+
+  cnts= 0;
+  for (unsigned i = 0; i < ihist; ++i) {
+    if(hist.at(i) > 0) cnts++;
+    cout<<"("<<xcor.at(i)<<","<<hist.at(i)<<") ";
+  }
+  cout<< " ... cnts=  " << cnts << endl;
+}