source: PSPA/Interface_Web/trunk/pspaWT/sources/controler/src/softwareParmela.cc @ 359

Last change on this file since 359 was 359, checked in by garnier, 11 years ago

Changement de la couleur des sections selectionnes + plein dautres choses
File size: 11.5 KB
Line 
1#include "softwareParmela.h"
2#include "abstractElement.h"
3#include "parmelaParticle.h"
4#include "mathematicalConstants.h"
5#include "PhysicalConstants.h"
6#include "dataManager.h"
7#include "mixedTools.h"
8
9softwareParmela::softwareParmela() : abstractSoftware()
10{;}
11
12softwareParmela::softwareParmela(string inputFileName, globalParameters* globals, dataManager* dt) : abstractSoftware(inputFileName, globals, dt)
13{;}
14
15bool softwareParmela::createInputFile(particleBeam* beamBefore,unsigned int numeroDeb,unsigned int numeroFin,string workingDir)
16{
17    unsigned int k;
18   
19    if ( numeroDeb < 1 || numeroFin > dataManager_->getBeamLineSize() ) {
20        dataManager_->consoleMessage(" softwareParmela::createInputFile : numero of element out of limits " );
21        cerr << " numero of element out of limits " << endl;
22        return false;
23    }
24   
25    string name = workingDir + inputFileName_;
26    ofstream outfile;
27    outfile.open(name.c_str(), ios::out);
28    if (!outfile) {
29        dataManager_->consoleMessage(" softwareParmela::createInputFile : error opening output stream " );
30        cerr << " softwareParmela::createInputFile : error opening output stream " << name << endl;
31        return false;
32    }
33   
34    abstractElement* elPtr;
35    double initalKineticEnergy = 0.0;
36    elPtr = dataManager_->getElementPointerFromNumero(numeroDeb);
37    bool there_is_rfgun = ( elPtr->getNomdElement().getElementType() == RFgun );
38   
39    if ( !there_is_rfgun ) {
40        string nameOfInput = workingDir + "parin.input0";
41        if ( !beamToParmela(nameOfInput,beamBefore) ) return false;
42        initalKineticEnergy = beamBefore->referenceKineticEnergyMeV();
43        // les elements de parmela sont indexes de 1 à max, s'il n'y a pas de rfgun
44        offsetNumElem_ = numeroDeb -1;
45    } else {
46        elPtr->setPhaseStep( globParamPtr_->getIntegrationStep() );
47        initalKineticEnergy = elPtr->getInitialKineticEnergy();
48        // les elements de parmela sont indexes de 0 à max, s'il y a un rfgun
49        offsetNumElem_ = numeroDeb;
50    }
51   
52    outfile << "TITLE" << endl;
53    outfile << " titre provisoire " << endl;
54    outfile << "RUN /n0=1 /ip=999 /freq=" << globParamPtr_->getFrequency() << "  /z0=0.0 /W0=" << initalKineticEnergy << "  /itype=1" << endl;
55    outfile << "OUTPUT 0" << endl;
56    unsigned int premier = numeroDeb ;
57    if ( there_is_rfgun ) {
58        outfile << dataManager_->getElementPointerFromNumero(numeroDeb)->parmelaOutputFlow() << endl;
59        premier++;
60    } else {
61        outfile << "INPUT 0 /NP=" << beamBefore->getNbParticles() << endl;
62    }
63   
64    // commentaire :  si l'element est un snapshot ne rien ecrire dans inputFileName_ (= parmin) un saut de ligne perturbe l'execution de parmela
65    for ( k = premier; k <= numeroFin; k++)
66    {
67        elPtr = dataManager_->getElementPointerFromNumero(k);
68        if (elPtr) {
69            if(elPtr->getNomdElement().getElementType() == snapshot) continue;
70            outfile << elPtr->parmelaOutputFlow() << endl;
71        }
72    }
73   
74    outfile << "ZOUT" << endl;
75    outfile << "START /wt=0.0 /dwt=" << globParamPtr_->getIntegrationStep() << "  /nsteps=" << globParamPtr_->getNbSteps() << "  nsc=" << globParamPtr_->getScPeriod() << "  /nout=10" << endl;
76    outfile << "END" << endl;
77    outfile.close();
78    dataManager_->consoleMessage("fichier input termine pour PARMELA");
79    return true;
80}
81
82bool softwareParmela::execute(unsigned int numeroDeb,unsigned int numeroFin,string workingDir)
83{
84    bool ExecuteStatus = true;
85   
86    ostringstream sortie;
87    sortie << " EXECUTION DE PARMELA DE l'ELEMENT " << numeroDeb << " A L'ELEMENT " << numeroFin << endl;
88   
89    string parmelaJob = workingDir + "parmela";
90    parmelaJob += string("   ");
91    parmelaJob += workingDir;
92   
93    string resultOfRun;
94    bool success = launchJob(parmelaJob,resultOfRun);
95    sortie << resultOfRun << endl;
96    if ( !success ) {
97        sortie << " launching of parmela failed " << endl;
98        ExecuteStatus = false;
99    } else {
100        sortie << " successful launching of parmela " << endl;
101        cout << " execution parmela MARCHE " << endl;
102        string::size_type nn = (resultOfRun).find("normal");
103        if ( nn == string::npos )
104        {
105            sortie << " abnormal exit of parmela " << endl;
106            ExecuteStatus = false;
107        }
108    }
109   
110    dataManager_->consoleMessage(sortie.str());
111    return ExecuteStatus;
112}
113
114bool softwareParmela::buildBeamAfterElements(unsigned int numeroDeb,unsigned int numeroFin, vector<particleBeam>& beams, string workingDir)
115{
116    bool result = true;
117    cout << "debug:: debut " << numeroDeb << ", fin " << numeroFin << endl;
118   
119    // index du premier element de parmela
120    int id= numeroDeb-offsetNumElem_;
121   
122    for(unsigned k = numeroDeb; k <= numeroFin; k++)
123    {
124        abstractElement* elem = dataManager_->getElementPointerFromNumero(k);
125       
126        // si l'element est un snapshot, recuperer la sortie precedente
127        if(elem->getNomdElement().getElementType() == snapshot) {
128            int avantDernier = beams.size() - 1;
129            string* param = elem->getParametersString();
130            string cliche = workingDir + param[2].c_str();
131            if( beamToParmela(cliche,&beams.at(avantDernier)) ) {
132                cout<<"["<<k<<"] : ecrit sur fichier " << cliche << " le contenu de beam["<<avantDernier<<"]"<<endl;
133            }
134            continue;
135        }
136       
137        // sinon c'est un element de parmela d'index id
138        beams.push_back(particleBeam());
139        cout << " creation diagn. no " << beams.size() << " par PARMELA " << endl;
140        vector<double> centroid;
141        bareParticle refPart;
142        vector<bareParticle> particles;
143        double freq= globParamPtr_->getFrequency();
144       
145        if(!beamFromParmela(workingDir,id,freq,centroid,refPart,particles))
146        {
147            if(elem->is_accepted_by_software(nomDeLogiciel::parmela) == TBoolIgnore) {
148                int avantDernier = beams.size() -2;
149                beams.back() = beams.at(avantDernier);
150            } else {
151                dataManager_->consoleMessage("softwareParmela::buildBeamAfterElements : reading parmdesz failed");
152                result = false;
153                break;
154            }
155           
156        } else {
157            beams.back().setWithParticles(centroid,refPart,particles);
158        }
159       
160        // l'element de parmela suivant
161        id++;
162    }
163   
164    return result;
165}
166
167bool softwareParmela::beamFromParmela(string workingDir,unsigned numeroParmel, double referencefrequency, vector<double>& centroid, bareParticle& refPart,vector<bareParticle>& particles )
168{
169   
170    string nomfilefais = workingDir + "parmdesz";
171    cout << " nom fichier desz : " << nomfilefais << endl;
172    FILE *filefais = fopen(nomfilefais.c_str(), "r");
173   
174    if ( filefais == (FILE*)0 ) {
175        dataManager_->consoleMessage(" beamFromParmela() erreur a l'ouverture du fichier 'parmdesz'");
176        cerr << " beamFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;
177        return false;
178    } else
179        cout << " beamFromParmela() : ouverture du fichier " << nomfilefais << endl;
180   
181    parmelaParticle partic;
182    std::vector<parmelaParticle> faisceau;
183    int testNombrePartRef =0;
184    double phaseRef = 0.0;
185   
186    while( partic.readFromParmelaFile(filefais) > 0 )
187    {
188       
189        if ( partic.ne == (int)numeroParmel ) {
190            if ( partic.np == 1 ) {
191                // en principe on est sur la particule de reference
192                if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON) {
193                    printf(" ATTENTION part. reference douteuse  \n");
194                    partic.imprim();
195                }
196                phaseRef = partic.phi;
197               
198                // le 'z' est 'absolu' (le long de la trajectoire)
199                TRIDVECTOR  posRef(partic.xx,partic.yy,partic.z);
200                TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
201                refPart = bareParticle(posRef, betagammaRef);
202                testNombrePartRef++;
203                if ( testNombrePartRef != 1 ) {
204                    dataManager_->consoleMessage(" beamFromParmela : nombre de part. de ref different de 1 :");
205                    cerr << " nombre de part. de ref different de 1 : " << testNombrePartRef << " !! " << endl;
206                    return false;
207                }
208            }
209            faisceau.push_back(partic);
210        }
211    } //while
212   
213    if ( faisceau.size() == 0)
214    {
215        string stringNum = mixedTools::intToString( (int)numeroParmel );
216        dataManager_->consoleMessage("beamFromParmela echec lecture  element numero relatif  parmela : " + stringNum);
217        cerr << " beamFromParmela echec lecture  element numero relatif  parmela " << numeroParmel << endl;
218        return false;
219    }
220   
221    // facteur  c/ 360. pour calculer (c dphi) / (360.freq)
222    // avec freq en Mhz et dphi en degres et résultat en cm:
223    double FACTEUR =  83.3333;  // ameliorer la precision
224   
225    // pour l'instant on choisit un centroid nul;
226    centroid.clear();
227    centroid = vector<double>(6,0.0);
228   
229    particles.clear();
230    particles.resize(faisceau.size(), bareParticle());
231    double x,xp,y,yp;
232    double betagammaz;
233    double betaz;
234    double deltaz;
235    double dephas;
236    double g;
237    TRIDVECTOR  pos;
238    TRIDVECTOR betagamma;
239    // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp
240    // sont donnes en radians
241    for (unsigned k = 0; k < faisceau.size(); k++) {
242        x= faisceau.at(k).xx;
243        xp=faisceau.at(k).xxp;
244        y= faisceau.at(k).yy;
245        yp=faisceau.at(k).yyp;
246        // dephasage par rapport a la reference
247        dephas = faisceau.at(k).phi - phaseRef; // degrés
248        g = faisceau.at(k).wz/EREST_MeV;
249        betagammaz = faisceau.at(k).begamz;
250        betaz = betagammaz/(g+1.0);
251        deltaz = FACTEUR * betaz * dephas / referencefrequency;
252        x += xp * deltaz;
253        y += yp * deltaz;
254        pos.setComponents(x,y,deltaz);
255        betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
256        particles.at(k) = bareParticle(pos,betagamma);
257    }
258   
259    return true;
260}
261
262// sauvegarde d'un 'particleBeam' sur un fichier parmela, en guise d'INPUT
263// pour l'instant de nom standard 'parin.input0'
264bool softwareParmela::beamToParmela(string nameOfFile,particleBeam* beam)
265{
266    if ( !beam->particleRepresentationOk() ) {
267        dataManager_->consoleMessage("softwareParmela::beamToParmela : beam not in particles form : not yet programmed");
268        cout << " softwareParmela::beamToParmela : beam not in particles form : not yet programmed " << endl;
269        return false;
270    }
271   
272    ofstream outfile;
273    outfile.open(nameOfFile.c_str(),ios::out);
274    if (!outfile) {
275        dataManager_->consoleMessage(" softwareParmela::beamToParmela : error opening output stream ");
276        cerr << " softwareParmela::beamToParmela : error opening output stream " << nameOfFile << endl;
277        return false;
278    }
279   
280    const vector<bareParticle>& partic = beam->getParticleVector();
281    double weight = 1.0;
282    double xx,yy,zz;
283    double begamx, begamy, begamz;
284    for (unsigned k = 0; k < partic.size(); k++) {
285        partic.at(k).getPosition().getComponents(xx,yy,zz);
286        partic.at(k).getBetaGamma().getComponents(begamx,begamy,begamz);
287        outfile << xx << " " << begamx << " " <<  yy << " " << begamy << " " << zz << " " << begamz  << " " << weight << endl;
288    }
289    outfile.close();
290    return true;
291}
Note: See TracBrowser for help on using the repository browser.