/[PAMELA software]/calo/flight/CaloBragg/src/CaloBragg.cpp
ViewVC logotype

Diff of /calo/flight/CaloBragg/src/CaloBragg.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.7 by mocchiut, Mon Dec 14 14:37:58 2009 UTC revision 1.11 by mocchiut, Wed Jun 15 12:28:45 2011 UTC
# Line 26  CaloBragg::CaloBragg(PamLevel2 *l2p){   Line 26  CaloBragg::CaloBragg(PamLevel2 *l2p){  
26    debug = false;    debug = false;
27    usetrack = false;    usetrack = false;
28    usepl18x = false;    usepl18x = false;
29      newchi2 = false;
30      usenewBB = false;
31      fzeta = -1.;
32    //    //
33  };  };
34    
35  void CaloBragg::Clear(){  void CaloBragg::Clear(){
36    //    //
37      ndf = 0;
38    tr = 0;    tr = 0;
39    sntr = 0;    sntr = 0;
40  //   qtchi2 = 0.;    //   qtchi2 = 0.;
41  //   qtz = 0.;    //   qtz = 0.;
42  //   qtetot = 0.;    //   qtetot = 0.;
43  //   qtpskip = 0.;    //   qtpskip = 0.;
44    lpchi2 = 0.;    lpchi2 = 0.;
45    lpz = 0.;    lpz = 0.;
46    lpetot = 0.;    lpetot = 0.;
47    lppskip = 0.;    lppskip = 0.;
48    memset(calorimetro,0,44*2*sizeof(Float_t));    memset(calorimetro,0,44*2*sizeof(Float_t));
49    memset(spessore,0,3*sizeof(Float_t));    memset(spessore,0,4*sizeof(Float_t));
50    memset(estremi,0,2*2*sizeof(Float_t));    memset(estremi,0,2*2*sizeof(Float_t));
51    Integrale=0.;    Integrale=0.;
52    
# Line 61  void CaloBragg::Print(){ Line 65  void CaloBragg::Print(){
65    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
66    printf(" first plane: %f \n", estremi[0][0]);    printf(" first plane: %f \n", estremi[0][0]);
67    printf(" last plane: %f \n", estremi[1][0]);    printf(" last plane: %f \n", estremi[1][0]);
68  //   printf(" chi 2 from truncated mean: %f \n", qtchi2);    //   printf(" chi 2 from truncated mean: %f \n", qtchi2);
69  //   printf(" Z from truncated mean %f: \n", qtz);    //   printf(" Z from truncated mean %f: \n", qtz);
70  //   printf(" energy from truncated mean %f: \n", qtetot);    //   printf(" energy from truncated mean %f: \n", qtetot);
71  //   printf(" plane not used for truncated mean %f: \n", qtpskip);    //   printf(" plane not used for truncated mean %f: \n", qtpskip);
72    printf(" chi 2 from loop %f:  \n", lpchi2);    printf(" chi 2 from loop %f:  \n", lpchi2);
73    printf(" Z from loop %f: \n", lpz);    printf(" Z from loop %f: \n", lpz);
74    printf(" energy from loop %f: \n", lpetot);    printf(" energy from loop %f: \n", lpetot);
75    printf(" plane not used for loop %f: \n", lppskip);    printf(" plane not used for loop %f: \n", lppskip);
76      printf(" ndf: %i \n",ndf);
77    printf("========================================================================\n");    printf("========================================================================\n");
78    //    //
79  };  };
# Line 83  void CaloBragg::Process(){ Line 88  void CaloBragg::Process(){
88    Process(-1);    Process(-1);
89  };  };
90    
91    
92    void CaloBragg::CleanPlanes(Float_t epiano[22][2]){
93      //  return;
94      Int_t hitplanes = 0;
95      for (Int_t i = 0; i<22; i++){
96        for (Int_t j = 1; j>=0; j--){
97          if ( epiano[i][j] > 0.7 ) hitplanes++;
98        };
99      };
100      Float_t lowlim = 0.85;
101      Float_t dedxone = 0.;
102      Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes;
103      while ( dedxone < step1 ){
104        for (Int_t i = 0; i<22; i++){
105          for (Int_t j = 1; j>=0; j--){
106            if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j];
107          };
108        };    
109      }
110      if ( dedxone < 0.7 ){
111        for (Int_t i = 0; i<22; i++){
112          for (Int_t j = 1; j>=0; j--){
113            if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j];
114          };
115        };    
116      }
117      //
118      //  printf(" dedxone = %f step1 %f  \n",dedxone,step1);
119      Bool_t revulsera = false;
120      Bool_t nullius = false;
121      Int_t nulliferus = 0;
122      for (Int_t i = 0; i<22; i++){
123        for (Int_t j = 1; j>=0; j--){
124          if ( epiano[i][j] < dedxone*lowlim ){
125            //        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
126            epiano[i][j] = 0.;
127          } else {
128            //x        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
129            nulliferus = 0;
130            revulsera = true;
131          };
132          if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++;
133          if ( nulliferus > 10 ) nullius = true;
134          if ( nullius ) epiano[i][j] = 0.;
135        };
136      };
137    
138    }
139    
140    
141  void CaloBragg::Process(Int_t ntr){  void CaloBragg::Process(Int_t ntr){
142    //      //  
143    if ( !L2 ){    if ( !L2 ){
# Line 133  void CaloBragg::Process(Int_t ntr){ Line 188  void CaloBragg::Process(Int_t ntr){
188      //      //
189    };    };
190    //    //
191      this->CleanPlanes(*&epiano);
192    //    //
193    PamTrack *ptrack = 0;    PamTrack *ptrack = 0;
194    CaloTrkVar *track = 0;    CaloTrkVar *track = 0;
# Line 164  void CaloBragg::Process(Int_t ntr){ Line 220  void CaloBragg::Process(Int_t ntr){
220    for(Int_t p=0; p<22; p++){    for(Int_t p=0; p<22; p++){
221      for(Int_t v=0; v<2; v++){      for(Int_t v=0; v<2; v++){
222        /*per usare traccia non del calo camboare cibar*/        /*per usare traccia non del calo camboare cibar*/
223        calorimetro[(2*p)+1-v][0] = L2->GetCaloLevel2()->cibar[p][v];//strip attraversata        calorimetro[(2*p)+1-v][0] = L2->GetCaloLevel2()->cibar[p][v];//strip attraversata      
224        calorimetro[(2*p)+1-v][1] = (epiano[p][v]); //energia del piano //(epiano[p][v])/0.89        calorimetro[(2*p)+1-v][1] = epiano[p][v]; //energia del piano //(epiano[p][v])/0.89
225      };      };
226    };    };
227    
# Line 181  void CaloBragg::Process(Int_t ntr){ Line 237  void CaloBragg::Process(Int_t ntr){
237        
238    //ordino tutte le energie dei piani in ordine crescente    //ordino tutte le energie dei piani in ordine crescente
239    
240      Long64_t work[200];    Long64_t work[200];
241      Int_t ind = 0;    Int_t ind = 0;
242      //Int_t l = 0;    //Int_t l = 0;
243      Int_t RN = 0;    Int_t RN = 0;
244      Float_t sum4 = 0.;    Float_t sum4 = 0.;
245      Float_t qm = 0.;    Float_t qm = 0.;
246      while ( RN < 4 && ind < 44 ){    while ( RN < 4 && ind < 44 ){
247        qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work);      qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work);
248        if (qm >= 0.7 ){        if (qm >= 0.7 ){    
249          if ( RN < 4 ){        if ( RN < 4 ){
250            sum4 += qm;          sum4 += qm;
251            RN++;          RN++;
         };  
252        };        };
       ind++;  
253      };      };
254      //      ind++;
255      //sum4 /= (Float_t)RN;    };
256      Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));    //
257      if(Zmean ==0.) Zmean=1.;    //sum4 /= (Float_t)RN;
258      if ( Zmean < 1. ) Zmean = 1.;    Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));
259      if(Zmean ==0.) Zmean=1.;
260      if ( Zmean < 1. ) Zmean = 1.;
261    
262    
263    /*trova primo e ultimo piano attraversati*/    /*trova primo e ultimo piano attraversati*/
# Line 216  void CaloBragg::Process(Int_t ntr){ Line 272  void CaloBragg::Process(Int_t ntr){
272      };      };
273      p++;      p++;
274    };    };
275   //ultimo parte da 44 e sale  
276      //ultimo parte da 44 e sale
277    p=43;    p=43;
278    while( (estremi[1][1] <= 0.)  &&  (p>(int)estremi[0][0]) ){    while( (estremi[1][1] <= 0.)  &&  (p>(int)estremi[0][0]) ){
279      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){
# Line 226  void CaloBragg::Process(Int_t ntr){ Line 283  void CaloBragg::Process(Int_t ntr){
283      p = p-1;      p = p-1;
284    };    };
285    //    //
286    
287      Float_t lastok = 0.;
288      //  Bool_t goback = false;
289      for ( int o = 0; o < estremi[1][0]; o++ ){
290        //
291        if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
292        if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
293        //    if ( calorimetro[o][1] < 0.7 ) goback = true;
294        //
295      };
296      lastok = 0.;
297      //  if ( goback ){
298      for ( int o = estremi[1][0]; o >= 0;  o-- ){
299        //
300        //    printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
301        if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok;  
302        if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
303        if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
304        //    printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
305        //
306      };
307      //  };
308    
309      if ( startZero ) {
310        estremi[0][0] = 0.;
311        //    estremi[0][1] = 0.;
312      }
313        
314    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
315    for(Int_t pl=0; pl<(2*NPLA); pl++){    for(Int_t pl=0; pl<(2*NPLA); pl++){
316        //    printf(" integrale: calorimetro %f  \n",calorimetro[pl][1]);
317      //calcolo intergale in unita di spessori di silicio      //calcolo intergale in unita di spessori di silicio
318      Integrale += calorimetro[pl][1] * MIP;//piano di silicio      Integrale += calorimetro[pl][1] * MIP;//piano di silicio
319      // se non e'il 1o dopo l'Y (tutti i pari) c'e' il W                // se non e'il 1o dopo l'Y (tutti i pari) c'e' il W          
# Line 236  void CaloBragg::Process(Int_t ntr){ Line 321  void CaloBragg::Process(Int_t ntr){
321        Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]);        Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]);
322      };      };
323    };    };
324    Integrale=24000;//Integrale*1000;    //Integrale=24000;//Integrale*1000;
325      Integrale *= 1000.;
326        
327    /*z ed energia con media troncata*/    /*z ed energia con media troncata*/
328    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
329        
330    /*z ed energia con loop*/    /*z ed energia con loop*/
331      Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip    if ( debug ) printf(" call Zdaloop with integrale %f \n",Integrale);
332      Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
333    
334        
335    if ( debug ) this->Print();    if ( debug ) this->Print();
# Line 251  void CaloBragg::Process(Int_t ntr){ Line 338  void CaloBragg::Process(Int_t ntr){
338  };  };
339    
340    
341    Float_t CaloBragg::Integral(){
342      Process();
343    
344      Float_t dEpianiloop[44];
345      Int_t tz1=(Int_t)lpz;
346      Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
347    
348    
349      Float_t integ = 0.;
350      for(Int_t i=0;i<=estremi[1][0];i++){
351        //    integ += dEplan[i];
352        //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
353        integ += dEpianiloop[i];
354        //    printf(" step %i integ %f deplan %f \n",i,integ,dEpianiloop[i]);
355      }
356      return integ;
357    }
358    
359    Float_t CaloBragg::LastIntegral(){
360      Process();
361    
362      Float_t integ = 0.;
363      for(Int_t i=0;i<=estremi[1][0];i++){
364        integ += dEplan[i];
365        //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
366      }
367      return integ;
368    }
369    
370  void CaloBragg::Draw(){  void CaloBragg::Draw(){
371    
372    Process();    Process();
373    
374  //  Float_t dEpianimean[44];    //  Float_t dEpianimean[44];
375   Float_t dEpianiloop[44];    Float_t dEpianiloop[44];
376   Float_t Depth[44];    Float_t Depth[44];
377  //  Int_t tz=(Int_t)qtz;    //  Int_t tz=(Int_t)qtz;
378   Int_t tz1=(Int_t)lpz;    Int_t tz1=(Int_t)lpz;
379  //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata    //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
380   Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop    Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
381    
382   Float_t sp= spessore[0]*spessore[1];    Float_t sp= spessore[0]*spessore[1];
383   for(Int_t i=0;i<44;i++)Depth[i]=i*sp;    for(Int_t i=0;i<44;i++)Depth[i]=i*sp;
384    //    //
385    gStyle->SetLabelSize(0.04);    gStyle->SetLabelSize(0.04);
386    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
387    //    //
388          TString hid = Form("cCaloBragg");          TString hid = Form("cCaloBragg");    
389          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));    TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
390          if ( tc ){    if ( tc ){
391  //       tc->Clear();      //   tc->Clear();
392          } else {    } else {
393           tc = new TCanvas(hid,hid);      tc = new TCanvas(hid,hid);
394  //       tc->Divide(1,2);      //   tc->Divide(1,2);
395          };    };
396          //    //
397  //      TString thid = Form("hCaloBragg");          //    TString thid = Form("hCaloBragg");      
398  //         TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));    //         TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
399  //      if ( th ) th->Delete();    //    if ( th ) th->Delete();
400  //       th->Clear();    //     th->Clear();
401  //       th->Reset();    //     th->Reset();
402  //      } else {    //    } else {
403  //      th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);    //    th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);
404  //      th->SetMarkerStyle(20);    //    th->SetMarkerStyle(20);
405  //      };    //    };
406          //    //
407          TString thid2 = Form("hCaloBragg2");        tc->cd();
408          TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));    TString thid2 = Form("hCaloBragg2");  
409          if ( th2 ) th2->Delete();    TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));
410          th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);    if ( th2 ) th2->Delete();
411          th2->SetMarkerStyle(20);    th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);
412          th2->SetMarkerColor(kRed);    th2->SetMarkerStyle(20);
413          //    th2->SetMarkerColor(kRed);
414          TString thid3 = Form("hCaloBragg3");        //
415          TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));    TString thid3 = Form("hCaloBragg3");  
416          if ( th3 ) th3->Delete();    TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));
417          th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);    if ( th3 ) th3->Delete();
418          th3->SetMarkerStyle(20);    th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);
419          th3->SetMarkerColor(kBlue);    th3->SetMarkerStyle(20);
420      th3->SetMarkerColor(kBlue);
421    
422          tc->cd(1);  
423          //    //  tc->cd(1);
424  //      for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);    //
425          for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);    //    for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);
426  //      th->Draw();    for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);
427          th2->Draw("same");    //    th->Draw();
428      th2->Draw("same");
429          tc->cd(2);  
430          //    //  tc->cd(2);
431          for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]);    tc->cd();
432          th3->Draw();    //
433          th2->Draw("same");    for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]);
434      th3->Draw();
435      th2->Draw("same");
436    
437          tc->Modified();    tc->Modified();
438          tc->Update();    tc->Update();
439    
440    //    //
441    gStyle->SetLabelSize(0);    gStyle->SetLabelSize(0);
# Line 331  void CaloBragg::LoadParam(){ Line 449  void CaloBragg::LoadParam(){
449    
450    //    //
451    elem[0] = 1.00794; //H  1    elem[0] = 1.00794; //H  1
452    elem[1] = 4.0026;  //He 2    elem[1] = 4.002602;  //He 2
453    elem[2] = 6.941;   //Li 3    elem[2] = 6.941;   //Li 3
454    elem[3] = 9.012182;//Be 4    elem[3] = 9.012182;//Be 4
455    elem[4] = 10.811;  //B  5    elem[4] = 10.811;  //B  5
456    elem[5] = 12.0107; //C  6    elem[5] = 12.0107; //C  6
457    elem[6] = 14.00674;//N  7    elem[6] = 14.00674;//N  7
458    elem[7] = 15.9994; //O  8    elem[7] = 15.9994; //O  8
459    elem[8] = 18.9984; //F  9    elem[8] = 18.9984032; //F  9
460    elem[9] = 20.1797; //Ne 10    elem[9] = 20.1797; //Ne 10
461    elem[10] = 22.98977;//Na 11    elem[10] = 22.98977;//Na 11
462    elem[11] = 24.3050; //Mg 12    elem[11] = 24.3050; //Mg 12
# Line 364  void CaloBragg::LoadParam(){ Line 482  void CaloBragg::LoadParam(){
482    elem[31] = 72.61; //Ge 32    elem[31] = 72.61; //Ge 32
483    
484    
485  //parametri calorimetro    //parametri calorimetro
486    NPLA = 22;    NPLA = 22;
487    NCHA = 96;    NCHA = 96;
488    nView = 2;    nView = 2;
# Line 385  void CaloBragg::LoadParam(){ Line 503  void CaloBragg::LoadParam(){
503    pigr = 3.1415;    pigr = 3.1415;
504    Na = 6.02e-23;    Na = 6.02e-23;
505    ZA = 0.49; /*Z/A per Si*/    ZA = 0.49; /*Z/A per Si*/
506    ISi =182e-06; /*MeV*/    //  ISi =182e-06; /*MeV*/
507      ISi = 171e-06; /*MeV*/
508      IW  = 735e-06; /*MeV*/
509      //  ISi =0.0001059994; /*GeV!!*/ no era giusto!!
510    Me = 0.511; /* MeV*/    Me = 0.511; /* MeV*/
511    MassP = 931.27;/*MeV*/    MassP = 931.27;/*MeV*/
512    r2 = 7.95e-26; /*ro*ro in cm */    r2 = 7.95e-26; /*ro*ro in cm */
# Line 420  void CaloBragg::conversione(){ Line 541  void CaloBragg::conversione(){
541    //lunghezza effettiva di silicio attraversata (mm)    //lunghezza effettiva di silicio attraversata (mm)
542    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
543        
544    spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2    spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2
545        
546    /*tungsteno*/    /*tungsteno*/
547    ytgx = yW * L2->GetCaloLevel2()->tanx[0];    ytgx = yW * L2->GetCaloLevel2()->tanx[0];
# Line 430  void CaloBragg::conversione(){ Line 551  void CaloBragg::conversione(){
551    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
552    //gcm2W = WCross/10. * rhoW;    //gcm2W = WCross/10. * rhoW;
553    
  a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);  
   
554    //       (g/cm2W)/(g/cm2Si)    //       (g/cm2W)/(g/cm2Si)
555      spessore[3] =  (WCross/10.) * rhoW;
556      a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
557    spessore[1] =  a;    spessore[1] =  a;
   
558    //riscala mip allo spessore attraversato    //riscala mip allo spessore attraversato
559    spessore[2] = MIP*(SiCross/ySi);      spessore[2] = MIP*(SiCross/ySi);  
   
560  };//end conversione  };//end conversione
561    
562    
563    
564    
565    
566  void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out){  void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out, Float_t II){
567    
568    //rilascio energetico con bethe bloch con correzioni    //rilascio energetico con bethe bloch con correzioni
569    //in:    x: g/cm2    //in:    x: g/cm2
# Line 462  void CaloBragg::BetheBloch(Float_t *x, F Line 581  void CaloBragg::BetheBloch(Float_t *x, F
581    Float_t lg =0.;    Float_t lg =0.;
582    Float_t Energia=0.;    Float_t Energia=0.;
583    Float_t C=0.;    Float_t C=0.;
584      Float_t INo = ISi;
585    
586      if ( usenewBB ) INo = II;
587    
588    eta = (*gam)*(*Bet);    eta = (*gam)*(*Bet);
589    
590    //Bet=3/gam;  SQ(*gam) * SQ(*Bet)    //Bet=3/gam;  SQ(*gam) * SQ(*Bet)
591    Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass));    Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass));
592        
593    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);    lg = 2.* Me * SQ(eta) * Wmax / SQ(INo);
594                                                                                                 //  Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA*  SQ(z)/SQ(Bet) * lg;    //  Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA*  SQ(z)/SQ(Bet) * lg;
595    C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(ISi,2.) +    C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(INo,2.) +
596      (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(ISi,3.);      (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(INo,3.);
597        
598    if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653);          if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653);      
599    
600     Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.);    Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.);
601    
602    *out =Energia;//out    *out =Energia;//out
603    
# Line 484  void CaloBragg::BetheBloch(Float_t *x, F Line 606  void CaloBragg::BetheBloch(Float_t *x, F
606    
607    
608    
609  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out){  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out, Float_t II){
610        
611    /*perdita di energia per ioni pesanti (come da routine geant)*/    /*perdita di energia per ioni pesanti (come da routine geant)*/
612    //  in : dx    => spessore g/cm2    //  in : dx    => spessore g/cm2
# Line 506  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 628  void CaloBragg::ELOSS(Float_t *dx, Int_t
628    
629    Bet = sqrt((SQ(gam) -1.)/SQ(gam));    Bet = sqrt((SQ(gam) -1.)/SQ(gam));
630        
631    v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.))));    //  v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.))));
632      v= 121.4139*(Bet*pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet*pow((*Z),(2./3.)))); // EMI AAAAGGH!!
633    
634    //carica effettiva    //carica effettiva
635    Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v));    Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v));
636    
637    //perdita energia per un protone    //perdita energia per un protone
638    Float_t protone =1.;    Float_t protone =1.;
639    Float_t Mass=(elem[*Z-1]*MassP);    //  Float_t Mass=(elem[*Z-1]*MassP); //EMI
640    BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP);    //  BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP);
641    
642      BetheBloch(dx, &protone, &MassP, &gam, &Bet, &dEP, II);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); //EMI
643        
644    *out= (SQ(Q)*(dEP));//*dx;    *out= (SQ(Q)*(dEP));//*dx;
645    
# Line 547  void CaloBragg::Enetrack(Int_t* Z, Float Line 672  void CaloBragg::Enetrack(Int_t* Z, Float
672    for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){    for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
673      dE=0.;      dE=0.;
674      //spessore silicio corretto x inclinazione, z, energia, out:rilascio      //spessore silicio corretto x inclinazione, z, energia, out:rilascio
675      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!      ELOSS(&spessore[0], Z, &Ezero, &dE, ISi);//spessore in g/cm2!!
676      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
677        out[ipla] = Ezero - Massa; //MeV        out[ipla] = Ezero - Massa; //MeV
678        return;        return;
# Line 555  void CaloBragg::Enetrack(Int_t* Z, Float Line 680  void CaloBragg::Enetrack(Int_t* Z, Float
680      }else{      }else{
681        out[ipla] = dE; //MeV        out[ipla] = dE; //MeV
682        Ezero = Ezero - dE;//energia residua        Ezero = Ezero - dE;//energia residua
683          if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero);
684      };      };
685      //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno      //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno
686      if(ipla%2 == 0){      if(ipla%2 == 0){
687        /*tungsteno*/        /*tungsteno*/
688        dE=0.;        dE=0.;
689        Float_t sp= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2        Float_t sp = 0.;
690        ELOSS(&sp, Z, &Ezero, &dE);        Float_t II = ISi;
691          if ( usenewBB ){
692            sp = spessore[3];        
693            II = IW;
694          } else {
695            sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2
696          }
697          //      printf(" sp %f II %f \n",sp,II);
698          ELOSS(&sp, Z, &Ezero, &dE,II);
699        if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop        if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
700          return;          return;
701        }else{        }else{
# Line 591  void CaloBragg::chiquadro(Float_t dE[], Line 725  void CaloBragg::chiquadro(Float_t dE[],
725    Float_t badplane=0.;    Float_t badplane=0.;
726    Float_t badplanetot=0.;    Float_t badplanetot=0.;
727    Float_t w,wi;    Float_t w,wi;
728      //
729    for(Int_t ipla=0; ipla<2*NPLA; ipla++){    if ( newchi2 ){
730      //tutti i piani attraversati dalla traiettoria      ndf = 0;
731      if(calorimetro[ipla][0] != -1.){ //      sum = 0.;
732        w=0.; //normalizzazione;      for( Int_t ipla=((int)(estremi[0][0])); ipla<= ((int)(estremi[1][0])); ipla++){
733        wi=1.;//peso        sum += pow((dE[ipla] - (calorimetro[ipla][1] * spessore[2]))/(0.05*dE[ipla]),2.);
734          //      printf(" quiqui: dE %f calor %f spessore[2] %f \n",dE[ipla],spessore[2]*calorimetro[ipla][1],spessore[2]);
735          ndf++;
736        }
737        ndf -= 2;
738        if ( ndf > 0 ) sum /= (float)ndf;
739        out[0] = sum;
740        out[1] = 0.;
741        out[2] = (int)(estremi[1][0])-ndf;
742        //    printf(" sum %f ndf %i \n ",sum,ndf);
743      } else {
744        for(Int_t ipla=0; ipla<2*NPLA; ipla++){
745          //tutti i piani attraversati dalla traiettoria
746          if(calorimetro[ipla][0] != -1.){ //
747            w=0.; //normalizzazione;
748            wi=1.;//peso
749    
750        //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso          //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso
751        if (ipla<estremi[0][0])  wi=0.;          if (ipla<estremi[0][0])  wi=0.;
752                
753        //tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0          //tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0
754        //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;          //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;
755        if((ipla>estremi[1][0])) wi=0.;          if((ipla>estremi[1][0])) wi=0.;
756    
757        //normalizzazione          //normalizzazione      
758        if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);    //            if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);    //        
759                
760        //tolgo piani con rilasci inferiori al 30% del precedente          //tolgo piani con rilasci inferiori al 30% del precedente
761        if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2          if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2
762          wi=0.;            wi=0.;
763          //se sono piani intermedi (non si e' fermta) li considero non buoni            //se sono piani intermedi (non si e' fermta) li considero non buoni
764          if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){//            if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){//
765            badplane+=1.;              badplane+=1.;
766            badplanetot+=1.;              badplanetot+=1.;
767          };            };
768        };          };
769        
770        //meno peso ai piani con rilasci maggiori di 1000 MIP          //meno peso ai piani con rilasci maggiori di 1000 MIP
771          if(calorimetro[ipla][1] > 1000) wi=0.5;          //      if(calorimetro[ipla][1] > 1000) wi=0.5;
772            if(calorimetro[ipla][1] > 1200.) wi=0.5;
773                            
774        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));          Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
775    
776        sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP));          sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP));
777        if(debug){          if(debug){
778          printf("dedx  calcolata  %f e reale  %f  \n",dE[ipla],(calorimetro[ipla][1] * MIP));            printf("dedx  calcolata  %f e reale  %f  \n",dE[ipla],(calorimetro[ipla][1] * MIP));
779        }          }
780        //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente          //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente
781        if(wi != 0.){//          if(wi != 0.){//
782          PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente            PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente
783          badplane = 0.;//azzero contatore piani scartati consecutivi            badplane = 0.;//azzero contatore piani scartati consecutivi
784            };
785        };        };
     };  
786                        
787      //da Emi        //da Emi
788      if(badplane > 2){              if(badplane > 2){      
789        out[1] =79.;          //      printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n");
790        break;          out[1] =79.;
791      };          break;
792            };
   };//fine loop piani  
   //chi2,frammentato,pskip    
   out[0]=sum;  
   out[2]=badplanetot;  
793        
794        };//fine loop piani
795        //chi2,frammentato,pskip  
796        out[0]=sum;
797        out[2]=badplanetot;
798      }  
799  };//end chiquadro  };//end chiquadro
800    
801    
802    
803  void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite){  void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite, Int_t nostep = 1000){
804  //    //
805    //loop su z ed energie per trovare miglior z (ed energia)    //loop su z ed energie per trovare miglior z (ed energia)
806    //in:  nloop     => energia massima da provare (nloop x E0)    //in:  nloop     => energia massima da provare (nloop x E0)
807    //     E0        => energia iniziale (intergale)            //     E0        => energia iniziale (intergale)        
# Line 661  void CaloBragg::loopze( Float_t step, Fl Line 812  void CaloBragg::loopze( Float_t step, Fl
812    //    //
813    
814    
   Float_t dEplan[2*NPLA];//energia rilasciata calcolata  
815    memset(dEplan,0,2*NPLA*sizeof(Float_t));    memset(dEplan,0,2*NPLA*sizeof(Float_t));
816    
817    Int_t Z = 0;// z iniziale    Int_t Z = 0;// z iniziale
818    
819    Float_t Massa = 0.;    Float_t Massa = 0.;
820    
821    Float_t Stepint =(step)/1000.;//passo per il calcolo di energia    Float_t Stepint =(step)/(Float_t)nostep;//passo per il calcolo di energia
822    
823    Float_t energia =0.;//energia del loop    Float_t energia =0.;//energia del loop
824    
825    Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati    Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati
826      
827      Int_t zmin = (int)Zstart;  
828    Int_t max=32;//max z di cui so la massa :P    Int_t max=32;//max z di cui so la massa :P
829    if((Zlimite)<=31) max=(int)(Zlimite) + 1;    if((Zlimite)<=31) max=(int)(Zlimite) + 1;
830    
831      if ( fzeta > 0. ){
832        zmin = fzeta;
833        max = fzeta+1;
834      }
835        
836    Int_t colmax=32;    Int_t colmax=32;
837    Int_t rowmax=3000;    Int_t rowmax=3000;
# Line 683  void CaloBragg::loopze( Float_t step, Fl Line 839  void CaloBragg::loopze( Float_t step, Fl
839    Float_t matrixchi2[colmax][rowmax][3];    Float_t matrixchi2[colmax][rowmax][3];
840    memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));    memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));
841    
842      Int_t imin = 1-nostep/2;
843      Int_t imax = nostep/2;
844    
845    //loop elementi      //loop elementi  
846    for(Int_t inucl=(int)(Zstart); inucl<max; inucl++){    for(Int_t inucl=zmin; inucl<max; inucl++){
847            
848      Z= inucl;      Z= inucl;
849        
850      Massa = elem[inucl-1]*MassP;      Massa = elem[inucl-1]*MassP;
851            
852      //loop energia      //loop energia
853      for(Int_t iene= 0; iene<1000; iene++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=??      Int_t iene2 = 0;
854            //    for(Int_t iene= 0; iene<1000; iene++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=??
855        for(Int_t iene= imin; iene<imax; iene++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=??
856          
857          iene2++;
858        energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..        energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
859    
860        Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani        Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
         
861        //calcolo chi2        //calcolo chi2
862        chiquadro(dEplan,chi2);        chiquadro(dEplan,chi2);
863    
864          //      printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]);
865                
866        if( (chi2[1] != 79.) ){//salto quelli che frammentano        if( (chi2[1] != 79.) ){//salto quelli che frammentano
867          matrixchi2[inucl][iene][0]=chi2[0];//valore chi2 per questo z a questa energia          matrixchi2[inucl][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia
868          matrixchi2[inucl][iene][1]=energia;//energia per questo chi2          matrixchi2[inucl][iene2][1]=energia;//energia per questo chi2
869          matrixchi2[inucl][iene][2]=chi2[2];//piani saltati nel chi2          matrixchi2[inucl][iene2][2]=chi2[2];//piani saltati nel chi2
870        } else {          } else {  
871          matrixchi2[inucl][iene][0]=1000;//valore chi2 per questo z a questa energia          matrixchi2[inucl][iene2][0]=1000.;//valore chi2 per questo z a questa energia
872          matrixchi2[inucl][iene][1]=1000;//energia per questo chi2          matrixchi2[inucl][iene2][1]=1000.;//energia per questo chi2
873          matrixchi2[inucl][iene][2]=1000;//piani saltati nel chi2          matrixchi2[inucl][iene2][2]=1000.;//piani saltati nel chi2
874          break;          break;
875        }        }
876      }//fine loop energia      }//fine loop energia
877    
878    
879   };//fine loop z    };//fine loop z
880    
881    
882    //Emi    //Emi
883    for (Int_t nu=(int)(Zstart); nu<max; nu++){    for (Int_t nu=zmin; nu<max; nu++){
884      for (Int_t en=0; en<1000; en++){      for (Int_t en=0; en<nostep; en++){
885        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){
886          bestchi2[0]= matrixchi2[nu][en][0];// chi2          bestchi2[0]= matrixchi2[nu][en][0];// chi2
887          bestchi2[1]= (Float_t)nu; // z          bestchi2[1]= (Float_t)nu; // z
# Line 845  void CaloBragg::Zdaloop(){ Line 1007  void CaloBragg::Zdaloop(){
1007    
1008    
1009    /*z se particella fosse al minimo*/  //energia1piano/mip corretta    /*z se particella fosse al minimo*/  //energia1piano/mip corretta
1010    Float_t zmax = round(sqrt(estremi[0][1]/spessore[2]));    //  Float_t zmax = round(sqrt(estremi[0][1]/spessore[2]));
1011    if(zmax<31)zmax=zmax+1;    //  if(zmax<31)zmax=zmax+1;
1012        
1013    /*calcolo Z ed E con loop sui vari elementi ed energie*/    /*calcolo Z ed E con loop sui vari elementi ed energie*/
1014    
1015    Float_t zmin=1.;    Float_t zmin=1.;
1016      Float_t zmax=32.;
1017    Float_t bestchitemp[4] = {0,0,0,0};    Float_t bestchitemp[4] = {0,0,0,0};
1018    
1019    bestchi2[0]=10000.;    bestchi2[0]=10000.;
# Line 860  void CaloBragg::Zdaloop(){ Line 1023  void CaloBragg::Zdaloop(){
1023    Float_t zero=0.;    Float_t zero=0.;
1024    //------------primo loop   ----------------------    //------------primo loop   ----------------------
1025    //     energia   ezero, zstart  zstop    //     energia   ezero, zstart  zstop
1026    loopze(Integrale,zero,zmin,zmax);    //  loopze(Integrale,zero,zmin,zmax);
1027    
1028      //->  loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50);
1029      loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200);
1030    
1031      //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);
1032      if ( debug ) printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]);
1033        
1034    //------------secondo loop  ----------------------    //------------secondo loop  ----------------------
1035    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];
# Line 871  void CaloBragg::Zdaloop(){ Line 1039  void CaloBragg::Zdaloop(){
1039    bestchi2[3] = 0.;//riazzero    bestchi2[3] = 0.;//riazzero
1040    
1041    Float_t step = bestchitemp[2];//    Float_t step = bestchitemp[2];//
1042    zero=0;  // qualsiasi altro valore peggiora le cose    zero=0.;  // qualsiasi altro valore peggiora le cose
1043    zmin=zmax=bestchitemp[1];    //  zmin=zmax=bestchitemp[1];
1044    loopze(step,zero,zmin,zmax); //    zmin=bestchitemp[1]-1;
1045      zmax=bestchitemp[1]+1;
1046      //  loopze(step,zero,zmin,zmax); //
1047    
1048      //->  loopze(step,step/2.,zmin,zmax,200); //
1049      loopze(step,step/2.,zmin,zmax,500); //
1050    
1051      if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step);
1052    
1053        
1054    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.23