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

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

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