Changeset 362 in PSPA for Interface_Web/trunk/pspaWT/sources/controler/src

Timestamp:

Mar 5, 2013, 11:48:58 AM (12 years ago)

Author:

touze

Message:

histogramme en x, y ou z

File:

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

Interface_Web/trunk/pspaWT/sources/controler/src/particleBeam.cc

@@ -510 +510 @@
 }
 
+void particleBeam::histogramme(int iabs,vector<double>&xcor,vector<int>& hist,int& cnts,double out[3])
+{ 
+  vector<double> vshf;
+  preparationDesHistogrammes(iabs,vshf,out);
+  remplissageDesHistogrammes(out,vshf,xcor,hist,cnts);
+
+  // sortie pour la legende: out[0]= entries, out[1]= mean, out[2]= std deviation
+  if(iabs < 3) {
+    out[1] *= 1.0 ; // en m ?
+    out[2] *= 1000.0; //mm
+  } else {
+    out[1] *= EREST_MeV ; //MeV
+    out[2] *= EREST_keV; //KeV
+  }
+}
+
+void particleBeam::preparationDesHistogrammes(int index,vector<double>& vshf,double out[3])
+{
+  double vmin= GRAND;
+  double vmax= -vmin;
+  double vmoy= 0.0;
+  double ecatyp= 0.0;
+
+  double sum= (float)goodPartic_.size();
+  if(index < 3) 
+    {
+      // collecte d'une coordonnee x(index= 0), y(index= 1) ou z(index= 2) des particules
+      for (unsigned int k = 0; k < goodPartic_.size(); k++) {
+        double val = goodPartic_.at(k).getPosition().getComponent(index);
+        if (val < vmin) vmin = val;
+        else if (val > vmax) vmax = val;
+        vmoy += val;
+        ecatyp += val*val;
+      }
+
+      out[0]= sum;
+      vmoy /= sum;
+      out[1]= vmoy;
+      ecatyp /= sum;
+      ecatyp = sqrt(abs(ecatyp-vmoy*vmoy));
+      out[2]= ecatyp;
+    } 
+  else 
+    {
+      // collecte de l'energie acquise par les particules
+      for (unsigned int k = 0; k < goodPartic_.size(); k++) {
+        double gamma = goodPartic_.at(k).getGamma();
+        if (gamma < vmin) vmin = gamma;
+        else if (gamma > vmax) vmax = gamma;
+        double val = gamma-1.0;
+        vmoy += val;
+        ecatyp += val*val;
+      }
+      
+      out[0]= sum;
+      vmoy /= sum;
+      out[1]= vmoy;
+      ecatyp /= sum;
+      ecatyp = sqrt(abs(ecatyp-vmoy*vmoy));
+      out[2]= ecatyp;
+    }
+  
+  if(index == 0) {
+    cout << "position X -moyenne " << vmoy << " m " << "-mini " << vmin << " m " << "-maxi " << vmax << " m " << "ecart type " << ecatyp*1000.0 << " mm" << endl;
+  } else if(index == 1) {
+    cout << "position Y -moyenne " << vmoy << " m " << "-mini " << vmin << " m " << "-maxi " << vmax << " m " << "ecart type " << ecatyp*1000.0 << " mm" << endl;
+  } else if(index == 2) {
+    cout << "position Z -moyenne " << vmoy << " m " << "-mini " << vmin << " m " << "-maxi " << vmax << " m " << "ecart type " << ecatyp*1000.0 << " mm" << endl;
+  } else {
+    double gmin = vmin-1.0;
+    double gmax = vmax-1.0;
+    cout << "energie cinetique -moyenne " << vmoy*EREST_MeV << " Mev " << "-mini " << gmin*EREST_MeV << " Mev " << "-maxi " << gmax*EREST_MeV << " Mev " << "ecart type " << ecatyp*EREST_MeV << " Kev" << endl;
+  }
+
+  //pouVoir: autre maniere de calculer l'ecart-type
+  double stdev= 0.0;
+
+  if(!vshf.empty()) vshf.clear();
+  vshf.resize(goodPartic_.size());  
+  if(index < 3) 
+    {
+      for (unsigned int k = 0; k < goodPartic_.size(); k++) {
+        double val = goodPartic_.at(k).getPosition().getComponent(index);
+        vshf[k]= val-vmoy;
+        stdev += vshf[k]*vshf[k];
+      }
+    } 
+  else 
+    {
+      for (unsigned int k = 0; k < goodPartic_.size(); k++) {
+        double gamma = goodPartic_.at(k).getGamma();
+        double val = gamma-1.0;
+        vshf[k]= val-vmoy;
+        stdev += vshf[k]*vshf[k];
+      }
+    }
+
+  stdev /= sum;
+  stdev = sqrt(stdev);
+  cout << "stdev " << stdev << ", ecatyp " << ecatyp << ", diff " << stdev-ecatyp << endl;
+}
+
+void particleBeam::remplissageDesHistogrammes(double out[3],vector<double>vshf,vector<double>&xcor,vector<int>& hist,int& cnts)
+{
+  // recupere la moyenne et l'ecart-type
+  double vmoy= out[1];
+  double ecatyp= out[2];
+
+  // demi fenetre hfene= max(3.*ecatyp-vmoy,vmoy-3.*ecatyp);
+  double hfene= 3.*ecatyp;
+  // et pas de l'histogramme
+  double hpas = hfene/25.;
+  
+  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= " << vshf.size() << ", phist " << phist << endl;
+
+  if(!xcor.empty()) xcor.clear();
+  xcor.resize(ihist);  
+  for (int i = 0; i < ihist; ++i) {
+
+    // on gradue l'abscisse en pourcents
+    xcor[i]= 100.*( vmin+i*hpas )/vmoy;
+  }
+
+  if(!hist.empty()) hist.clear();
+  hist.resize(ihist,0);  
+  for (unsigned i = 0; i < vshf.size(); ++i) {
+    double var= vshf[i]-vmin;
+    if(var < 0 ) {
+      cout<<"lesser that the minimum "<<var<<", ("<< i<<")"<< endl;
+      hist[0]++;
+    }
+
+    if(var >= phist) {
+      cout<<"greater that the maximum "<<var<<", ("<< i<<")"<< endl;
+      hist[ihist-1]++;
+    }
+
+    int kk= (int)floor(var/hpas);
+    hist[kk]++;
+  }
+
+  cnts= 0;
+  for (int i = 0; i < ihist; ++i) {
+    if(hist.at(i) > 0) cnts++;
+  }
+  cout<< " ... cnts=  " << cnts << endl;
+}
+
 void particleBeam::histogramme(vector<double>&xcor,vector<int>& hist,int& cnts,double out[3])
 {