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

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

bug
File size: 10.7 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 // cout << "fill " << indexDebutPiste << "\n";
119 int* PhDiod=PhDiodArray;
120 int* pDataC=pData;
121 for(int k=0;k<nData;k++,pDataC++, PhDiod++) (*PhDiod)=(*pDataC);
122 return;
123}
124
125void PisteEtoile::push(int* pData, int nbdata) {
126 int* PPhDiod=PhDiodArray;
127 int* PPhDiod2=PhDiodArray+nbdata;
128 int* pDataC=pData;
129
130 for(int i=0; i<PhDiodTabLong-nbdata; i++, PPhDiod++, PPhDiod2++) *PPhDiod=(*PPhDiod2);
131 for(int i=0; i<nbdata; i++, PPhDiod++,pDataC++) *PPhDiod=(*pDataC);
132 indexDebutPiste+=nbdata;
133 // cout << "push " << nbdata << " " << indexDebutPiste << "\n";
134 return;
135}
136
137bool PisteEtoile::traque() {
138 if(!isAlive()) return false; //La diode est-elle morte?
139
140 bool testEtoile=trig();
141 if(testEtoile) {
142 //On a declenché sur une Etoile
143 if(archParam.sst.analFine) {
144 demasque();
145// Un second Filtre est possible
146 // if(!(fabs(AmplRes)>SeuilGachette*Vrms)&&SuccesAFine) testEtoile=false;
147 // Test rajouté pour filtrer le bruit en escalier des voies SST
148 // Le fit de forme tasse le bruit et permet un meilleur filtrage.
149 }
150
151// Troisième filtre possible:
152 // On tue les evenements trop larges
153 if (widthHalfMax>(FPulse.FWHM()/PasTempsEch+2.)) testEtoile=false;
154
155 if(testEtoile) {
156 // On emballe une etoile.
157 // cout << "star " << setprecision(10) << TIndex << "\n";
158 LastEtoile.FitForme=SuccesAFine;
159 LastEtoile.X2Calc=X2Min;
160 LastEtoile.InpCurrNoisRMS=Vrms/FPulse.GainElectrique();
161 LastEtoile.TEchan=TIndex;
162 LastEtoile.InpCurrent= AmplRes/FPulse.GainElectrique();
163 LastEtoile.NoDiode=NoPhDiod;
164 LastEtoile.LargMiHauteur=widthHalfMax;
165 }
166 }
167 // On expedie
168 return testEtoile;
169}
170
171bool PisteEtoile::trig() {
172 MaxIndex=-1;
173 AmplMax=0.;
174 gachette=false;
175
176// Copie les données brutes sur la piste. Les met en Volt ADC 12 Bits-10 Volts
177 for(int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]=PhDiodArray[i]*SSTLsbVal;
178
179// Calcule Valeur Moy du prepulse
180 double Moy=0.;
181 for(int i=0; i<Prepulselong; i++) Moy+=PhDiodPiste[i];
182 Moy/=Prepulselong;
183
184// On soustrait la valeur Moy de la piste
185 for (int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]-=Moy;
186
187 double pente=0.;
188 if (archParam.sst.soustPente==true) {
189 // On calcule la pente du prepulse
190 double milprepulse=(Prepulselong-1.)/2.;
191 double SomCarre=0.;
192 for(int i=0; i<Prepulselong; i++) {
193 SomCarre+= (i-milprepulse)*(i-milprepulse);
194 pente+= (i-milprepulse)*PhDiodPiste[i];
195 }
196 pente=pente/SomCarre;
197
198 // On soustrait la pente sur la piste
199 for (int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]-=(i-milprepulse)*pente;
200 }
201
202// On calcule le Vrms sur le prepulse
203 Vrms=0.;
204 for (int i=0; i<Prepulselong; i++) Vrms+=PhDiodPiste[i]*PhDiodPiste[i];
205 Vrms=sqrt(Vrms)/Prepulselong;
206 if(Vrms<SSTLsbVal) Vrms=SSTLsbVal; // Pour eviter les VRMS=0.
207
208// Declenche-t-on?
209 gachette=false;
210 short iSeuil=0;
211 double* pPiste =PhDiodPiste+Prepulselong-1;
212
213// On demande deux échantillons à seuilgachette Vrms
214 for (int i=0; i<(PhDiodTabLong-Prepulselong); i++, pPiste++) {
215 if (fabs(*pPiste)>(SeuilGachette*Vrms))
216 if (fabs(*(pPiste+1))>SeuilGachette*Vrms) {
217 gachette=true;
218 iSeuil=Prepulselong-1+i;
219 break;
220 }
221 }
222
223 if (gachette==true) {
224// On trouve le maximum
225 MaxIndex=iSeuil;
226 AmplMax=fabs(PhDiodPiste[iSeuil]);
227 for (int i=iSeuil; i<PhDiodTabLong; i++) {
228 if (fabs(PhDiodPiste[i])>AmplMax) {
229 MaxIndex=i;
230 AmplMax=fabs(PhDiodPiste[i]);
231 }
232 }
233 if ((MaxIndex>=PhDiodTabLong-4)||(MaxIndex<Prepulselong+4)) gachette=false;
234// Cela sera traite a l'iteration suivante
235// Ou ca a du etre traite a l'iteration precedente
236 }
237 if (gachette==true) {
238// On calcule la largeur a mi hauteur
239 int HalfMaxInf;
240 for(HalfMaxInf=MaxIndex;HalfMaxInf>Prepulselong;HalfMaxInf--)
241 if(fabs(PhDiodPiste[HalfMaxInf])<(AmplMax/2.)) break;
242 int HalfMaxSup;
243 for(HalfMaxSup=MaxIndex;HalfMaxSup<PhDiodTabLong;HalfMaxSup++)
244 if(fabs(PhDiodPiste[HalfMaxSup])<(AmplMax/2.)) break;
245 widthHalfMax=HalfMaxSup-HalfMaxInf;
246
247/*
248// On integre sur la larg a mi hauteur
249 integ=0.;
250 if((HalfMaxSup<PhDiodTabLong)&&(HalfMaxInf>Prepulselong))
251 for(int i=HalfMaxInf; i<=HalfMaxSup; i++) integ+=PhDiodPiste[i];
252*/
253
254
255// On fitte une parabole autour du maximum
256 double tmaxfit;
257 if (widthHalfMax>1) {
258 double ymax=PhDiodPiste[MaxIndex];
259 double yplus=PhDiodPiste[MaxIndex+1];
260 double ymoins=PhDiodPiste[MaxIndex-1];
261 // Estimation du temps
262 if(!(ymax==yplus)) {
263 double Y=(yplus-ymax)/(ymax-ymoins);
264 tmaxfit=0.5*(1.+Y)/(1-Y);
265 }
266 else tmaxfit=0.5;
267
268 // Estimation de l'amplitude
269 double dtplus2=(1.+tmaxfit)*(1.+tmaxfit);
270 double dtmoins2=(1.-tmaxfit)*(1.-tmaxfit);
271 if(dtplus2==dtmoins2) AmplRes=PhDiodPiste[MaxIndex];
272 else AmplRes= (yplus*dtplus2-ymoins*dtmoins2) / (dtplus2-dtmoins2);
273 }
274
275// On stocke les resultats
276 TIndex=indexDebutPiste+MaxIndex+tmaxfit;
277 // cout << "tindex " << TIndex << " " <<indexDebutPiste<<" " <<Prepulselong << " " <<MaxIndex << " " << tmaxfit << "\n";
278 X2Min=0.;
279 SuccesAFine=false;
280 }
281 return gachette;
282}
283
284void PisteEtoile::demasque() {
285// Maintenant que l'on se doute que il y a une etoile
286// On veut connaitre son flux et son temps de passage
287// Nécessite d'avoir un ordre d'idee de la vitesse de rotation du ballon.
288
289 // BUGG devra etre extrait des headers acquisition. Ca peut changer!!!!!!
290
291
292
293 int IndexDebutEvt=MaxIndex-SSTLargEvtIndex();
294
295 pFitArray0= PhDiodPiste+IndexDebutEvt;
296 // Pointeur origine du tableau d'echantillon sur lequel sera calculé le fit
297 double DeltaTemps=0.;
298// Pour chaque decalage en temps, on estime un amplitude et calcule le X2
299 for(int k=0; k<SSTnbPasFit; k++) {
300 DeltaTemps=TFExcursion-k*TPasTFit;
301 pFitArray=pFitArray0;
302 double Som=0.;
303 for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++) {
304// cout<<(*pFitArray)<<'\t'<<FPulse.PulseShape(i*PasTempsEch-DeltaTemps)<<endl;
305 Som=Som+FPulse.PulseShape(i*PasTempsEch-DeltaTemps)*(*pFitArray);
306 }
307 Ampl[k]=Som/AmplNorm;
308 X2[k]=0.;
309 double X2Som=0.; double dum=0.;pFitArray=pFitArray0;
310 for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++) {
311 dum=(*pFitArray)-Ampl[k]*FPulse.PulseShape(i*PasTempsEch-DeltaTemps); // Différence
312 X2Som+=dum*dum; // Au carré
313 }
314 X2[k]=X2Som/(SSTLongIndexEvt()*Vrms*Vrms);
315 }
316
317// On recherche le minimum du X2
318 X2Min=1.e99;
319 int X2MinIndex=-1;
320 for(int k=0; k<SSTnbPasFit; k++)
321 if (X2[k]<X2Min) {
322 X2Min=X2[k];
323 X2MinIndex=k;
324 }
325
326 if((X2MinIndex>0)&&(X2MinIndex<(SSTnbPasFit-1))) {
327
328// On fitte une parabole sur les trois X2 qui encadrent le min
329 double Di; // Index (non entier) du minimum du X2
330 double X=(X2[X2MinIndex+1]-X2Min)/(X2Min-X2[X2MinIndex-1]);
331 Di=0.5*(X+1)/(X-1);
332
333// On calcule l'amplitude estimee a ce temps
334 double Som1=0.; pFitArray=pFitArray0;
335 DeltaTemps=TFExcursion-(X2MinIndex-Di)*TPasTFit;
336 for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++)
337 Som1+=(*pFitArray)*FPulse.PulseShape(i*PasTempsEch-DeltaTemps);
338 AmplRes=Som1/AmplNorm;
339
340//Le X2 est mimimum pour le n° de pas de fit X2MinIndex+Di
341
342#ifdef PEtoileDebug
343 double* PulseTrak=new double [SSTLongIndexEvt()];
344#endif
345
346 double X2Som1=0.; double dummy=0.;pFitArray=pFitArray0;
347 for(int i=0;i<SSTLongIndexEvt();i++,pFitArray++) {
348#ifdef PEtoileDebug
349 PulseTrak[i]=FPulse.PulseShape(i*PasTempsEch-DeltaTemps);
350#endif
351 dummy=(*pFitArray)-AmplRes*FPulse.PulseShape(i*PasTempsEch-DeltaTemps); // Différence
352 X2Som1+=dummy*dummy; // Au carré
353 }
354
355 X2Min=X2Som1/(SSTLongIndexEvt()*Vrms*Vrms);
356
357 TIndex=indexDebutPiste+IndexDebutEvt+DeltaTemps/PasTempsEch;
358 //Ouff!!!!
359 SuccesAFine=true;
360 }
361 else { // Echec, c'est tres bizarre! On doit etre dans le bruit
362 X2Min=0.;
363 }
364 return;
365}
366
367int PisteEtoile::SSTLargEvtIndex() {
368// Demi-largeur a 3 sigma de l'impulsion
369 int dum=(int) (15.e-3/PasTempsEch+1.);
370 return dum;
371}
372
373int PisteEtoile::SSTLongIndexEvt() {
374// Nombre d'echantillon sur lequel est calculé le fit de forme
375 int temp=2*SSTLargEvtIndex()+1;
376 return temp;
377}
378
Note: See TracBrowser for help on using the repository browser.