source: Sophya/trunk/Poubelle/archTOI.old/pisteetoile.cc@ 555

Last change on this file since 555 was 555, checked in by ansari, 26 years ago

portage cxx en cours
File size: 10.8 KB
RevLine 
[555]1#include <math.h>
[394]2#include <new.h>
[555]3#include <iostream.h>
[394]4#include "archparam.h"
5#include "pisteetoile.h"
6
7TransFuncElec PisteEtoile::Trapani;
8FormePulse PisteEtoile::FPulse(Trapani);
9
[534]10// D. Yvon, CE Saclay, DAPNIA/SPP, 08/99
[394]11
12PisteEtoile::PisteEtoile(short NoPhotDiod) {
[534]13 PasTempsEch=archParam.acq.perEch;
[394]14 //BUGG doit changer !!!!! Doit etre extrait d'un serveur de param d'acquisition
15 NoPhDiod=NoPhotDiod;
16 DiodVivante=true;
[534]17
18 SSTnbPasFit= (int) (2*TFExcursion/TPasTFit);
19
[394]20 AmplNorm=0.;
21 for(int i=0; i<SSTLongIndexEvt()*10; i++){
[534]22 AmplNorm+=FPulse.PulseShape(i*PasTempsEch/10.)*FPulse.PulseShape(i*PasTempsEch/10.);
[394]23 }
24 AmplNorm/=10.;
25 indexDebutPiste=0;
26
27 PhDiodArray=new(nothrow) int[PhDiodTabLong];
28 PhDiodPiste=new(nothrow) double[PhDiodTabLong];
29 Ampl=new(nothrow) double[SSTnbPasFit];
30 X2=new(nothrow) double[SSTnbPasFit];
31
32 if((X2==0)||(Ampl==0)||(PhDiodPiste==0)||(PhDiodArray==0)) {
33 cerr<<"Erreur a la reservation de memoire dans constructeur pisteetoile"<<endl;
34 exit(-1);
35 }
36}
37
38PisteEtoile& PisteEtoile::operator =(PisteEtoile const & x) {
39 PasTempsEch=x.PasTempsEch;
40 NoPhDiod=x.NoPhDiod;
41 DiodVivante=x.DiodVivante;
42 indexDebutPiste=x.indexDebutPiste;
43 DiodVivante=x.DiodVivante;
44 AmplNorm=x.AmplNorm;
45
46
47 for(int j=0; j<PhDiodTabLong;j++) PhDiodArray[j]=x.PhDiodArray[j];
48 for(int j=0; j<PhDiodTabLong;j++) PhDiodPiste[j]=x.PhDiodPiste[j];
49
50 LastEtoile=x.LastEtoile;
51
52 gachette=x.gachette;
53 Vrms=x.Vrms;
54 MaxIndex=x.MaxIndex; // Index de l'amplitude Max en analyse Trig()
55 AmplMax=x.AmplMax;
56
57 AmplRes=x.AmplRes;
58 TIndex=x.TIndex;
59 X2Min=x.X2Min;
60 SuccesAFine=x.SuccesAFine;
61 AmplNorm=x.AmplNorm;
62
63 return *this;
64}
65
66PisteEtoile::PisteEtoile (PisteEtoile const & x) {
67 PasTempsEch=x.PasTempsEch;
68 NoPhDiod=x.NoPhDiod;
69 DiodVivante=x.DiodVivante;
70 indexDebutPiste=x.indexDebutPiste;
71 DiodVivante=x.DiodVivante;
72 AmplNorm=x.AmplNorm;
73
74 PhDiodArray=new(nothrow) int[PhDiodTabLong];
75 PhDiodPiste=new(nothrow) double[PhDiodTabLong];
76 Ampl=new(nothrow) double[SSTnbPasFit];
77 X2=new(nothrow) double[SSTnbPasFit];
78
79
80 if((X2==0)||(Ampl==0)||(PhDiodPiste==0)||(PhDiodArray==0)) {
81 cerr<<"Erreur a la reservation de memoire dans constructeur pisteetoile"<<endl;
82 exit(-1);
83 }
84
85 for(int j=0; j<PhDiodTabLong;j++) PhDiodArray[j]=x.PhDiodArray[j];
86 for(int j=0; j<PhDiodTabLong;j++) PhDiodPiste[j]=x.PhDiodPiste[j];
87
88 LastEtoile=x.LastEtoile;
89
90 gachette=x.gachette;
91 Vrms=x.Vrms;
92 MaxIndex=x.MaxIndex; // Index de l'amplitude Max en analyse Trig()
93 AmplMax=x.AmplMax;
94
95 AmplRes=x.AmplRes;
96 TIndex=x.TIndex;
97 X2Min=x.X2Min;
98 SuccesAFine=x.SuccesAFine;
99 AmplNorm=x.AmplNorm;
100
101}
102
103
104
105PisteEtoile::~PisteEtoile() {
106 delete[] X2;
107 delete[] Ampl;
108 delete[] PhDiodPiste;
109 delete[] PhDiodArray;
110
111}
112
113void PisteEtoile::fill(int* pData, int InDebutPiste, int nData) {
114// if(!isAlive()) return;
[534]115 indexDebutPiste=InDebutPiste - PhDiodTabLong + nData;
116 int* PhDiod=PhDiodArray + PhDiodTabLong - nData;
[394]117 int* pDataC=pData;
118 for(int k=0;k<nData;k++,pDataC++, PhDiod++) (*PhDiod)=(*pDataC);
119 return;
120}
121
122void PisteEtoile::push(int* pData, int nbdata) {
123 int* PPhDiod=PhDiodArray;
124 int* PPhDiod2=PhDiodArray+nbdata;
125 int* pDataC=pData;
126
127 for(int i=0; i<PhDiodTabLong-nbdata; i++, PPhDiod++, PPhDiod2++) *PPhDiod=(*PPhDiod2);
128 for(int i=0; i<nbdata; i++, PPhDiod++,pDataC++) *PPhDiod=(*pDataC);
129 indexDebutPiste+=nbdata;
130 return;
131}
132
133bool PisteEtoile::traque() {
134 if(!isAlive()) return false; //La diode est-elle morte?
135
136 bool testEtoile=trig();
137 if(testEtoile) {
138 //On a declenché sur une Etoile
139 if(archParam.sst.analFine) {
140 demasque();
141// Un second Filtre est possible
142 // if(!(fabs(AmplRes)>SeuilGachette*Vrms)&&SuccesAFine) testEtoile=false;
143 // Test rajouté pour filtrer le bruit en escalier des voies SST
144 // Le fit de forme tasse le bruit et permet un meilleur filtrage.
145 }
146
147// Troisième filtre possible:
148 // On tue les evenements trop larges
149 if (widthHalfMax>(FPulse.FWHM()/PasTempsEch+2.)) testEtoile=false;
150
151 if(testEtoile) {
[534]152 // On emballe une etoile
153 LastEtoile.InpCurrent= AmplRes/FPulse.GainElectrique();
154 LastEtoile.InpCurFin= AmplFin/FPulse.GainElectrique();
155 LastEtoile.TEchan=TIndex;
156 LastEtoile.TFin=Tfin;
157 LastEtoile.InpCurrNoisRMS=Vrms/FPulse.GainElectrique();
[394]158 LastEtoile.X2Calc=X2Min;
159 LastEtoile.NoDiode=NoPhDiod;
[534]160 LastEtoile.FitForme=SuccesAFine;
[394]161 LastEtoile.LargMiHauteur=widthHalfMax;
[534]162
163/* cout<<setprecision(9);
164 cout<<TIndex<<" "<<Tfin<<" "<<AmplRes<<" "<<AmplFin<<endl;
165 cout<<TIndex-Tfin<<" "<<(AmplRes-AmplFin)/AmplFin<<endl;
166*/
[394]167 }
168 }
169 // On expedie
170 return testEtoile;
171}
172
173bool PisteEtoile::trig() {
174 MaxIndex=-1;
175 AmplMax=0.;
176 gachette=false;
177
178// Copie les données brutes sur la piste. Les met en Volt ADC 12 Bits-10 Volts
179 for(int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]=PhDiodArray[i]*SSTLsbVal;
180
181// Calcule Valeur Moy du prepulse
182 double Moy=0.;
183 for(int i=0; i<Prepulselong; i++) Moy+=PhDiodPiste[i];
184 Moy/=Prepulselong;
185
186// On soustrait la valeur Moy de la piste
187 for (int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]-=Moy;
188
189 double pente=0.;
190 if (archParam.sst.soustPente==true) {
191 // On calcule la pente du prepulse
192 double milprepulse=(Prepulselong-1.)/2.;
193 double SomCarre=0.;
194 for(int i=0; i<Prepulselong; i++) {
195 SomCarre+= (i-milprepulse)*(i-milprepulse);
196 pente+= (i-milprepulse)*PhDiodPiste[i];
197 }
198 pente=pente/SomCarre;
199
200 // On soustrait la pente sur la piste
201 for (int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]-=(i-milprepulse)*pente;
202 }
203
204// On calcule le Vrms sur le prepulse
205 Vrms=0.;
[534]206 for (int i=0; i<Prepulselong; i++) {
207 Vrms+=PhDiodPiste[i]*PhDiodPiste[i];
208 }
209 Vrms=sqrt(Vrms/Prepulselong);
[394]210 if(Vrms<SSTLsbVal) Vrms=SSTLsbVal; // Pour eviter les VRMS=0.
211
212// Declenche-t-on?
213 gachette=false;
214 short iSeuil=0;
[397]215 double* pPiste =PhDiodPiste+Prepulselong-1;
[394]216
[399]217// On demande deux échantillons à seuilgachette Vrms
218 for (int i=0; i<(PhDiodTabLong-Prepulselong); i++, pPiste++) {
[394]219 if (fabs(*pPiste)>(SeuilGachette*Vrms))
220 if (fabs(*(pPiste+1))>SeuilGachette*Vrms) {
221 gachette=true;
[399]222 iSeuil=Prepulselong-1+i;
[394]223 break;
224 }
225 }
226
227 if (gachette==true) {
228// On trouve le maximum
229 MaxIndex=iSeuil;
230 AmplMax=fabs(PhDiodPiste[iSeuil]);
231 for (int i=iSeuil; i<PhDiodTabLong; i++) {
232 if (fabs(PhDiodPiste[i])>AmplMax) {
233 MaxIndex=i;
234 AmplMax=fabs(PhDiodPiste[i]);
235 }
236 }
237 if ((MaxIndex>=PhDiodTabLong-4)||(MaxIndex<Prepulselong+4)) gachette=false;
238// Cela sera traite a l'iteration suivante
239// Ou ca a du etre traite a l'iteration precedente
240 }
241 if (gachette==true) {
242// On calcule la largeur a mi hauteur
243 int HalfMaxInf;
244 for(HalfMaxInf=MaxIndex;HalfMaxInf>Prepulselong;HalfMaxInf--)
245 if(fabs(PhDiodPiste[HalfMaxInf])<(AmplMax/2.)) break;
246 int HalfMaxSup;
247 for(HalfMaxSup=MaxIndex;HalfMaxSup<PhDiodTabLong;HalfMaxSup++)
248 if(fabs(PhDiodPiste[HalfMaxSup])<(AmplMax/2.)) break;
249 widthHalfMax=HalfMaxSup-HalfMaxInf;
250
251/*
252// On integre sur la larg a mi hauteur
253 integ=0.;
254 if((HalfMaxSup<PhDiodTabLong)&&(HalfMaxInf>Prepulselong))
255 for(int i=HalfMaxInf; i<=HalfMaxSup; i++) integ+=PhDiodPiste[i];
256*/
257
258
259// On fitte une parabole autour du maximum
260 double tmaxfit;
261 if (widthHalfMax>1) {
262 double ymax=PhDiodPiste[MaxIndex];
263 double yplus=PhDiodPiste[MaxIndex+1];
264 double ymoins=PhDiodPiste[MaxIndex-1];
265 // Estimation du temps
266 if(!(ymax==yplus)) {
267 double Y=(yplus-ymax)/(ymax-ymoins);
268 tmaxfit=0.5*(1.+Y)/(1-Y);
269 }
270 else tmaxfit=0.5;
271
272 // Estimation de l'amplitude
273 double dtplus2=(1.+tmaxfit)*(1.+tmaxfit);
274 double dtmoins2=(1.-tmaxfit)*(1.-tmaxfit);
275 if(dtplus2==dtmoins2) AmplRes=PhDiodPiste[MaxIndex];
276 else AmplRes= (yplus*dtplus2-ymoins*dtmoins2) / (dtplus2-dtmoins2);
277 }
278
279// On stocke les resultats
[399]280 TIndex=indexDebutPiste+MaxIndex+tmaxfit;
[394]281 X2Min=0.;
282 SuccesAFine=false;
283 }
284 return gachette;
285}
286
287void PisteEtoile::demasque() {
288// Maintenant que l'on se doute que il y a une etoile
289// On veut connaitre son flux et son temps de passage
[534]290// Nécessite de connaitre la vitesse de rotation du ballon.
[394]291
292 // BUGG devra etre extrait des headers acquisition. Ca peut changer!!!!!!
293
294 int IndexDebutEvt=MaxIndex-SSTLargEvtIndex();
295
296 pFitArray0= PhDiodPiste+IndexDebutEvt;
297 // Pointeur origine du tableau d'echantillon sur lequel sera calculé le fit
298 double DeltaTemps=0.;
[534]299
300// Pour chaque decalage en temps, on estime une amplitude et calcule le X2
[394]301 for(int k=0; k<SSTnbPasFit; k++) {
302 DeltaTemps=TFExcursion-k*TPasTFit;
303 pFitArray=pFitArray0;
304 double Som=0.;
305 for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++) {
306// cout<<(*pFitArray)<<'\t'<<FPulse.PulseShape(i*PasTempsEch-DeltaTemps)<<endl;
307 Som=Som+FPulse.PulseShape(i*PasTempsEch-DeltaTemps)*(*pFitArray);
308 }
309 Ampl[k]=Som/AmplNorm;
310 X2[k]=0.;
311 double X2Som=0.; double dum=0.;pFitArray=pFitArray0;
312 for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++) {
313 dum=(*pFitArray)-Ampl[k]*FPulse.PulseShape(i*PasTempsEch-DeltaTemps); // Différence
314 X2Som+=dum*dum; // Au carré
315 }
316 X2[k]=X2Som/(SSTLongIndexEvt()*Vrms*Vrms);
317 }
318
319// On recherche le minimum du X2
320 X2Min=1.e99;
321 int X2MinIndex=-1;
[534]322 for(int k=0; k<SSTnbPasFit; k++) {
[394]323 if (X2[k]<X2Min) {
324 X2Min=X2[k];
325 X2MinIndex=k;
326 }
[534]327 }
328
[394]329 if((X2MinIndex>0)&&(X2MinIndex<(SSTnbPasFit-1))) {
330
331// On fitte une parabole sur les trois X2 qui encadrent le min
[534]332 double Di; // Index (non entier) du minimum du X2
[394]333 double X=(X2[X2MinIndex+1]-X2Min)/(X2Min-X2[X2MinIndex-1]);
334 Di=0.5*(X+1)/(X-1);
335
[534]336// Le X2 est mimimum pour le n° de pas de fit X2MinIndex+Di
[394]337// On calcule l'amplitude estimee a ce temps
338 double Som1=0.; pFitArray=pFitArray0;
339 DeltaTemps=TFExcursion-(X2MinIndex-Di)*TPasTFit;
340 for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++)
[534]341 Som1+=(*pFitArray)*FPulse.PulseShape(i*PasTempsEch-DeltaTemps);
342 AmplFin=Som1/AmplNorm;
343
[394]344
[534]345
[394]346#ifdef PEtoileDebug
[534]347 double *PulseTrak=new double [SSTLongIndexEvt()];
[394]348#endif
[534]349// On calcule le X2 correspondant
[394]350 double X2Som1=0.; double dummy=0.;pFitArray=pFitArray0;
351 for(int i=0;i<SSTLongIndexEvt();i++,pFitArray++) {
[534]352
[394]353#ifdef PEtoileDebug
354 PulseTrak[i]=FPulse.PulseShape(i*PasTempsEch-DeltaTemps);
355#endif
[534]356 dummy=(*pFitArray)-AmplFin*FPulse.PulseShape(i*PasTempsEch-DeltaTemps); // Différence
[394]357 X2Som1+=dummy*dummy; // Au carré
358 }
[534]359
[394]360 X2Min=X2Som1/(SSTLongIndexEvt()*Vrms*Vrms);
[534]361
362 Tfin= indexDebutPiste+IndexDebutEvt+DeltaTemps/PasTempsEch;
363 //Ouff!!!!
364/*
365 cout<<setprecision(9)<<endl;
366 cout<<Tfin<<" "<<indexDebutPiste<<" "<<IndexDebutEvt<<" "<<DeltaTemps/PasTempsEch<<endl;
367*/
[394]368 SuccesAFine=true;
[534]369
370#ifdef PEtoileDebug
371 delete[] PulseTrak;
372#endif
[394]373 }
374 else { // Echec, c'est tres bizarre! On doit etre dans le bruit
375 X2Min=0.;
376 }
377 return;
378}
379
380int PisteEtoile::SSTLargEvtIndex() {
381// Demi-largeur a 3 sigma de l'impulsion
382 int dum=(int) (15.e-3/PasTempsEch+1.);
383 return dum;
384}
385
386int PisteEtoile::SSTLongIndexEvt() {
387// Nombre d'echantillon sur lequel est calculé le fit de forme
388 int temp=2*SSTLargEvtIndex()+1;
389 return temp;
390}
391
Note: See TracBrowser for help on using the repository browser.