/[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.10 by mocchiut, Wed Jun 15 10:04:59 2011 UTC revision 1.16 by mocchiut, Thu Jan 23 11:23:43 2014 UTC
# Line 43  void CaloBragg::Clear(){ Line 43  void CaloBragg::Clear(){
43    //   qtpskip = 0.;    //   qtpskip = 0.;
44    lpchi2 = 0.;    lpchi2 = 0.;
45    lpz = 0.;    lpz = 0.;
46      lpisotope= 0.;
47    lpetot = 0.;    lpetot = 0.;
48    lppskip = 0.;    lppskip = 0.;
49    
50    memset(calorimetro,0,44*2*sizeof(Float_t));    memset(calorimetro,0,44*2*sizeof(Float_t));
51    memset(spessore,0,3*sizeof(Float_t));    memset(spessore,0,4*sizeof(Float_t));
52    memset(estremi,0,2*2*sizeof(Float_t));    memset(estremi,0,2*2*sizeof(Float_t));
53    Integrale=0.;    Integrale=0.;
54    
# Line 71  void CaloBragg::Print(){ Line 73  void CaloBragg::Print(){
73    //   printf(" plane not used for truncated mean %f: \n", qtpskip);    //   printf(" plane not used for truncated mean %f: \n", qtpskip);
74    printf(" chi 2 from loop %f:  \n", lpchi2);    printf(" chi 2 from loop %f:  \n", lpchi2);
75    printf(" Z from loop %f: \n", lpz);    printf(" Z from loop %f: \n", lpz);
76      printf(" isotope from loop %f: \n", lpisotope);
77    printf(" energy from loop %f: \n", lpetot);    printf(" energy from loop %f: \n", lpetot);
78    printf(" plane not used for loop %f: \n", lppskip);    printf(" plane not used for loop %f: \n", lppskip);
79    printf(" ndf: %i \n",ndf);    printf(" ndf: %i \n",ndf);
# Line 89  void CaloBragg::Process(){ Line 92  void CaloBragg::Process(){
92  };  };
93    
94    
95  void CaloBragg::CleanPlanes(Float_t epiano[22][2]){  void CaloBragg::CleanPlanes(Float_t epiano[22][2], Bool_t zpiano[22][2]){
96    //  return;    //  return;
97    Int_t hitplanes = 0;    Int_t hitplanes = 0;
98      Float_t f5 = 0.;
99    for (Int_t i = 0; i<22; i++){    for (Int_t i = 0; i<22; i++){
100      for (Int_t j = 1; j>=0; j--){      for (Int_t j = 1; j>=0; j--){
101        if ( epiano[i][j] > 0.7 ) hitplanes++;        zpiano[i][j] = false;
102          if ( epiano[i][j] > 0.7 ){
103            if ( hitplanes < 100 ) f5 += epiano[i][j];
104            hitplanes++;
105          };
106      };      };
107    };    };
108      Int_t atl5 = TMath::Min(hitplanes,100);
109      atl5 = TMath::Max(atl5,1);
110    Float_t lowlim = 0.85;    Float_t lowlim = 0.85;
111      //Float_t lowlim = 1.;
112    Float_t dedxone = 0.;    Float_t dedxone = 0.;
113    Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes;    //  Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes;
114    while ( dedxone < step1 ){    Float_t step1 = 0.8*f5/atl5;
115      //  while ( dedxone < step1 ){
116      for (Int_t i = 0; i<22; i++){      for (Int_t i = 0; i<22; i++){
117        for (Int_t j = 1; j>=0; j--){        for (Int_t j = 1; j>=0; j--){
118            if (debug) printf("Acleanplanes: i %i j %i step1 %f dedxone %f epiano[i][j] %f \n",i,j,step1,dedxone,epiano[i][j]);
119          if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j];          if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j];
120            if ( dedxone >= step1 ) break; // new
121        };        };
122          if ( dedxone >= step1 ) break; // new
123      };          };    
124    }      //  }
125    if ( dedxone < 0.7 ){    if ( dedxone < 0.7 ){   // here we could have instead while dedxone == 0. ... perhaps better...
126      for (Int_t i = 0; i<22; i++){      for (Int_t i = 0; i<22; i++){
127        for (Int_t j = 1; j>=0; j--){        for (Int_t j = 1; j>=0; j--){
128            if (debug) printf("Bcleanplanes dedxone < 0.7: i %i j %i step1 %f dedxone %f epiano[i][j] %f \n",i,j,step1,dedxone,epiano[i][j]);
129          if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j];          if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j];
130            if ( dedxone >= 0.7 ) break; // new
131        };        };
132          if ( dedxone >= 0.7 ) break; // new
133      };          };    
134    }    }
135    //    //
# Line 122  void CaloBragg::CleanPlanes(Float_t epia Line 140  void CaloBragg::CleanPlanes(Float_t epia
140    for (Int_t i = 0; i<22; i++){    for (Int_t i = 0; i<22; i++){
141      for (Int_t j = 1; j>=0; j--){      for (Int_t j = 1; j>=0; j--){
142        if ( epiano[i][j] < dedxone*lowlim ){        if ( epiano[i][j] < dedxone*lowlim ){
143          //        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);          if ( debug ) printf("Ccleanplanes: %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
144          epiano[i][j] = 0.;          //        epiano[i][j] = 0.;
145            zpiano[i][j] = true;
146            if ( epiano[i][j] < dedxone*0.05 ) epiano[i][j] = 0.;
147        } else {        } else {
148          //x        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);          if ( debug ) printf("Dcleanplanes else: %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
149          nulliferus = 0;          nulliferus = 0;
150          revulsera = true;          revulsera = true;
151        };        };
152        if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++;        //      if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++;
153          if ( (zpiano[i][j] || epiano[i][j] < 0.7 ) && revulsera ) nulliferus++;
154        if ( nulliferus > 10 ) nullius = true;        if ( nulliferus > 10 ) nullius = true;
155        if ( nullius ) epiano[i][j] = 0.;        //      if ( nullius ) epiano[i][j] = 0.;
156          if ( nullius ) zpiano[i][j] = true;
157      };      };
158    };    };
159    
# Line 188  void CaloBragg::Process(Int_t ntr){ Line 210  void CaloBragg::Process(Int_t ntr){
210      //      //
211    };    };
212    //    //
213    this->CleanPlanes(*&epiano);    Bool_t zpiano[22][2];
214      this->CleanPlanes(*&epiano, *&zpiano);
215    //    //
216    PamTrack *ptrack = 0;    PamTrack *ptrack = 0;
217    CaloTrkVar *track = 0;    CaloTrkVar *track = 0;
# Line 217  void CaloBragg::Process(Int_t ntr){ Line 240  void CaloBragg::Process(Int_t ntr){
240    if(L2->GetCaloLevel2()->npcfit[0]==0 && L2->GetCaloLevel2()->npcfit[1]==0 && L2->GetCaloLevel2()->npcfit[2]==0 && L2->GetCaloLevel2()->npcfit[3]==0) return;// controllo sulla traccia nel calorimetro    if(L2->GetCaloLevel2()->npcfit[0]==0 && L2->GetCaloLevel2()->npcfit[1]==0 && L2->GetCaloLevel2()->npcfit[2]==0 && L2->GetCaloLevel2()->npcfit[3]==0) return;// controllo sulla traccia nel calorimetro
241    
242    //    //
243      Bool_t zcalo[44];
244    for(Int_t p=0; p<22; p++){    for(Int_t p=0; p<22; p++){
245      for(Int_t v=0; v<2; v++){      for(Int_t v=0; v<2; v++){
246        /*per usare traccia non del calo camboare cibar*/        /*per usare traccia non del calo camboare cibar*/
247        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      
248        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
249          zcalo[(2*p)+1-v] = zpiano[p][v];
250          if ( debug ) printf(" idx %i %f %i \n",(2*p)+1-v,epiano[p][v], zpiano[p][v]);
251      };      };
252    };    };
253    
# Line 285  void CaloBragg::Process(Int_t ntr){ Line 311  void CaloBragg::Process(Int_t ntr){
311    //    //
312    
313    Float_t lastok = 0.;    Float_t lastok = 0.;
314      //  if ( false ){
315    //  Bool_t goback = false;    //  Bool_t goback = false;
316    for ( int o = 0; o < estremi[1][0]; o++ ){    for ( int o = 0; o < estremi[1][0]; o++ ){
317      //      //
318      if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];      if (debug)    printf(" goforth1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
319      if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;          if ( calorimetro[o][1] > 0.7 && !zcalo[o] ) lastok = calorimetro[o][1];
320        if ( (zcalo[o] || calorimetro[o][1] < 0.7) && lastok > 0. ){
321          if ( fabs(calorimetro[o][1]-lastok)/calorimetro[o][1] > 0.5 ) {
322            if (debug) printf(" goforthchange %f %f \n",calorimetro[o][1],lastok);
323            calorimetro[o][1] = lastok;    
324            if (debug)        printf(" goforthchang+ %f %f \n",calorimetro[o][1],lastok);
325          }
326        }
327        if (debug)    printf(" goforth2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
328      //    if ( calorimetro[o][1] < 0.7 ) goback = true;      //    if ( calorimetro[o][1] < 0.7 ) goback = true;
329      //      //
330    };    };
# Line 297  void CaloBragg::Process(Int_t ntr){ Line 332  void CaloBragg::Process(Int_t ntr){
332    //  if ( goback ){    //  if ( goback ){
333    for ( int o = estremi[1][0]; o >= 0;  o-- ){    for ( int o = estremi[1][0]; o >= 0;  o-- ){
334      //      //
335      //    printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);      if (debug)    printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
336      if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok;        if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok;  
337      if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];      if ( calorimetro[o][1] > 0.7 && !zcalo[o] ) lastok = calorimetro[o][1];
338      if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;          if ( (zcalo[o] || calorimetro[o][1] < 0.7) && lastok > 0. ){
339      //    printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);        if ( fabs(calorimetro[o][1]-lastok)/calorimetro[o][1] > 0.5 ) {
340            if ( debug ) printf(" gobackchange \n");
341            calorimetro[o][1] = lastok;    
342          }
343        }
344        if (debug)    printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
345      //      //
346    };    };
347    //  };    //  };
348      //}
349    
350    if ( startZero ) {    if ( startZero ) {
351      estremi[0][0] = 0.;      estremi[0][0] = 0.;
# Line 343  Float_t CaloBragg::Integral(){ Line 384  Float_t CaloBragg::Integral(){
384    
385    Float_t dEpianiloop[44];    Float_t dEpianiloop[44];
386    Int_t tz1=(Int_t)lpz;    Int_t tz1=(Int_t)lpz;
387    Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop    Int_t ti1=(Int_t)lpisotope;
388    
389      Enetrack(&tz1, &ti1 , &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
390    
391    
392    Float_t integ = 0.;    Float_t integ = 0.;
# Line 367  Float_t CaloBragg::LastIntegral(){ Line 410  Float_t CaloBragg::LastIntegral(){
410    return integ;    return integ;
411  }  }
412    
413    
414  void CaloBragg::Draw(){  void CaloBragg::Draw(){
415    
416    Process();    Process();
417    
418      this->Draw(0.,0.,0.);
419    
420    }
421    
422    void CaloBragg::Draw(Int_t Z, Int_t isotope, Float_t enetot){
423    
424    //  Float_t dEpianimean[44];    //  Float_t dEpianimean[44];
425    Float_t dEpianiloop[44];    Float_t dEpianiloop[44];
426    Float_t Depth[44];    Float_t Depth[44];
427    //  Int_t tz=(Int_t)qtz;    //  Int_t tz=(Int_t)qtz;
428    Int_t tz1=(Int_t)lpz;    Int_t tz1= Z;
429    //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata    Int_t ti1= isotope;
430    Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop    Float_t enet = enetot;
431      //  Float_t enet = lpetot;
432      
433      if ( Z > 0. && enetot > 0. ){
434        estremi[0][0] = 0;
435        estremi[1][0] = 43;
436    
437    
438        Float_t    ytgx = 0.;
439        Float_t    ytgy = 0.;
440        
441        //lunghezza effettiva di silicio attraversata (mm)
442        Float_t    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
443        
444        spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2
445        
446        /*tungsteno*/
447        
448        //rapporto tra rilasci energetici nei due materiali
449        Float_t WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
450        //gcm2W = WCross/10. * rhoW;
451        
452        //       (g/cm2W)/(g/cm2Si)
453        spessore[3] =  (WCross/10.) * rhoW;
454        Float_t a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
455        spessore[1] =  a;
456        //riscala mip allo spessore attraversato
457        spessore[2] = MIP*(SiCross/ySi);  
458    
459      } else {
460        tz1=(Int_t)lpz;
461        ti1=(Int_t)lpisotope;
462        enet = lpetot;
463        //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
464      
465      }
466      Enetrack(&tz1, &ti1, &enet, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
467    
468    Float_t sp= spessore[0]*spessore[1];    Float_t sp= spessore[0]*spessore[1];
469    for(Int_t i=0;i<44;i++)Depth[i]=i*sp;    for(Int_t i=0;i<44;i++)Depth[i]=i*sp;
# Line 408  void CaloBragg::Draw(){ Line 494  void CaloBragg::Draw(){
494    TString thid2 = Form("hCaloBragg2");      TString thid2 = Form("hCaloBragg2");  
495    TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));    TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));
496    if ( th2 ) th2->Delete();    if ( th2 ) th2->Delete();
497    th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);    th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,25.); //150
498    th2->SetMarkerStyle(20);    th2->SetMarkerStyle(20);
499    th2->SetMarkerColor(kRed);    th2->SetMarkerColor(kRed);
500    //    //
501    TString thid3 = Form("hCaloBragg3");      TString thid3 = Form("hCaloBragg3");  
502    TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));    TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));
503    if ( th3 ) th3->Delete();    if ( th3 ) th3->Delete();
504    th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);    th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,25.);//150.
505    th3->SetMarkerStyle(20);    th3->SetMarkerStyle(20);
506    th3->SetMarkerColor(kBlue);    th3->SetMarkerColor(kBlue);
507    
# Line 430  void CaloBragg::Draw(){ Line 516  void CaloBragg::Draw(){
516    //  tc->cd(2);    //  tc->cd(2);
517    tc->cd();    tc->cd();
518    //    //
519    for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]);    for(Int_t i=0;i<=estremi[1][0];i++){
520        th3->Fill(Depth[i],dEpianiloop[i]);
521        //    printf(" i %i Depth %f depianiloop %f \n",i,Depth[i],dEpianiloop[i]);
522      }
523    th3->Draw();    th3->Draw();
524    th2->Draw("same");    th2->Draw("same");
525    
# Line 447  void CaloBragg::Draw(){ Line 536  void CaloBragg::Draw(){
536    
537  void CaloBragg::LoadParam(){  void CaloBragg::LoadParam(){
538    
539    //    // elem[Z-1][isotop] 0 is the most common one
540    elem[0] = 1.00794; //H  1    //
541    elem[1] = 4.002602;  //He 2  
542    elem[2] = 6.941;   //Li 3    elem[0][0] = 1.00782;   //H  1
543    elem[3] = 9.012182;//Be 4    elem[0][1] = 2.01410;   // 2H (Isotope)
544    elem[4] = 10.811;  //B  5    elem[0][2] = -1.;
545    elem[5] = 12.0107; //C  6    elem[0][3] = -1.;
546    elem[6] = 14.00674;//N  7    elem[0][4] = -1.;
547    elem[7] = 15.9994; //O  8    elem[0][5] = -1.;
548    elem[8] = 18.9984032; //F  9    elem[0][6] = -1.;
549    elem[9] = 20.1797; //Ne 10  
550    elem[10] = 22.98977;//Na 11    elem[1][0] = 4.002602;  //He 2
551    elem[11] = 24.3050; //Mg 12    elem[1][1] = 3.016029;  // 3He (Isotope)
552    elem[12] = 26.9815; //Al 13    elem[1][2] = -1.;
553    elem[13] = 28.0855; //Si 14    elem[1][3] = -1.;
554    elem[14] = 30.974;  //P  15    elem[1][4] = -1.;
555    elem[15] = 32.066;  //S  16    elem[1][5] = -1.;
556    elem[16] = 35.4527; //Cl 17    elem[1][6] = -1.;
557    elem[17] = 39.948;  //Ar 18  
558    elem[18] = 39.0983; //K  19    elem[2][0] = 7.016004;  //Li 3
559    elem[19] = 40.078;  //Ca 20    elem[2][1] = 6.015123;  //6Li (Isotope)
560    elem[20] = 44.95591;//Sc 21    elem[2][2] = -1.;
561    elem[21] = 47.867;  //Ti 22    elem[2][3] = -1.;
562    elem[22] = 50.9415; //V  23    elem[2][4] = -1.;
563    elem[23] = 51.9961; //Cr 24    elem[2][5] = -1.;
564    elem[24] = 54.938049;//Mn 25    elem[2][6] = -1.;
565    elem[25] = 55.845;  //Fe 26  
566    elem[26] = 58.9332; //Co 27    elem[3][0] = 9.012182;  //Be 4
567    elem[27] = 58.6934; //Ni 28    elem[3][1] = 10.01353;  //10Be (Isotope) (most stable)
568    elem[28] = 63.546; //Cu 29    elem[3][2] = 7.01693;   //9Be no EC in space?
569    elem[29] = 65.39; //Zn 30    elem[3][3] = -1.;
570    elem[30] = 69.723; //Ga 31    elem[3][4] = -1.;
571    elem[31] = 72.61; //Ge 32    elem[3][5] = -1.;
572      elem[3][6] = -1.;
573    
574      elem[4][0] = 11.00930;   //B 5
575      elem[4][1] = 10.01294;  //10B  (Isotope)
576      elem[4][2] = -1.;
577      elem[4][3] = -1.;
578      elem[4][4] = -1.;
579      elem[4][5] = -1.;
580      elem[4][5] = -1.;
581    
582      elem[5][0] = 12.0107;   //C  6
583      elem[5][1] = 13.00335;  //13C (Isotope)
584      elem[5][2] = -1.;
585      elem[5][3] = -1.;
586      elem[5][4] = -1.;
587      elem[5][5] = -1.;
588      elem[5][5] = -1.;
589    
590      elem[6][0] = 14.00674;  //N  7
591      elem[6][1] = 15.00011;  //15N  (Isotope)
592      elem[6][2] = -1.;
593      elem[6][3] = -1.;
594      elem[6][4] = -1.;
595      elem[6][5] = -1.;
596      elem[6][5] = -1.;
597    
598      elem[7][0] = 15.99491;  //O  8
599      elem[7][1] = 17.99916;  //18O  (Isotope)
600      elem[7][2] = 16.99916;  //17O  (Isotope)
601      elem[7][3] = -1.;
602      elem[7][4] = -1.;
603      elem[7][5] = -1.;
604      elem[7][5] = -1.;
605    
606      elem[8][0] = 18.99840;  //F  9
607      elem[8][1] = -1.;
608      elem[8][2] = -1.;
609      elem[8][3] = -1.;
610      elem[8][4] = -1.;
611      elem[8][5] = -1.;
612      elem[8][5] = -1.;
613    
614      elem[9][0] = 19.99244;  //Ne 10
615      elem[9][1] = 21.99138;  //22Ne (Isotope)
616      elem[9][2] = 20.99384;  //21Ne 10
617      elem[9][3] = -1.;
618      elem[9][4] = -1.;
619      elem[9][5] = -1.;
620      elem[9][6] = -1.;
621    
622      elem[10][0] = 22.98977; //Na 11
623      elem[10][1] = 21.99444; //22Na (Isotope) (most stable)
624      elem[10][2] = -1.;
625      elem[10][3] = -1.;
626      elem[10][4] = -1.;
627      elem[10][5] = -1.;
628      elem[10][6] = -1.;
629    
630      elem[11][0] = 23.98504; //Mg 12
631      elem[11][1] = 25.98259; //26Mg (Isotope)
632      elem[11][2] = 24.98504; //25Mg (Isotope)
633      elem[11][3] = -1.;
634      elem[11][4] = -1.;
635      elem[11][5] = -1.;
636      elem[11][6] = -1.;
637    
638      elem[12][0] = 26.98154; //Al 13
639      elem[12][1] = 25.98489; //26Al (Isotope) (most stable)
640      elem[12][2] = -1.;
641      elem[12][3] = -1.;
642      elem[12][4] = -1.;
643      elem[12][5] = -1.;
644      elem[12][6] = -1.;
645    
646      elem[13][0] = 27.97692; //Si 14
647      elem[13][1] = 28.97649; //29Si (Isotope)
648      elem[13][2] = 29.97377; //30Si (Isotope)
649      elem[13][3] = -1.;
650      elem[13][4] = -1.;
651      elem[13][5] = -1.;
652      elem[13][6] = -1.;
653    
654      elem[14][0] = 30.97376;  //P  15
655      elem[14][1] = -1.;
656      elem[14][2] = -1.;
657      elem[14][3] = -1.;
658      elem[14][4] = -1.;
659      elem[14][5] = -1.;
660      elem[14][6] = -1.;
661    
662      elem[15][0] = 31.97207;  //S  16
663      elem[15][1] = 33.96787;  //34S  (Isotope)
664      elem[15][2] = 32.97146;  //33S  (Isotope)
665      elem[15][3] = 35.96708;  //36S  (Isotope)
666      elem[15][4] = -1.;
667      elem[15][5] = -1.;
668      elem[15][6] = -1.;
669    
670      elem[16][0] = 34.96885; //Cl 17
671      elem[16][1] = 36.96831; //37Cl 17
672      elem[16][2] = 35.96890; //36Cl (Isotope)
673      elem[16][3] = -1.;
674      elem[16][4] = -1.;
675      elem[16][5] = -1.;
676      elem[16][6] = -1.;
677    
678      elem[17][0] = 39.962383;  //Ar 18
679      elem[17][1] = 35.967545;  //36Ar (Isotope)
680      elem[17][2] = 37.962732;  //38Ar (Isotope)
681      elem[17][3] = 38.964313;  //39Ar (Isotope)
682      elem[17][4] = -1.;
683      elem[17][5] = -1.;
684      elem[17][6] = -1.;
685    
686      elem[18][0] = 38.963707; //K  19
687      elem[18][1] = 40.961825; //41K  (Isotope)
688      elem[18][2] = 39.963998; //40K  (Isotope)
689      elem[18][3] = -1.;
690      elem[18][4] = -1.;
691      elem[18][5] = -1.;
692      elem[18][6] = -1.;
693    
694      elem[19][0] = 39.962590;  //Ca 20
695      elem[19][1] = 43.955482;  //44Ca (Isotope)
696      elem[19][2] = 41.958618;  //42Ca (Isotope)
697      elem[19][3] = 42.958767;  //43Ca (Isotope)
698      elem[19][4] = 45.953693;  //46Ca (Isotope)
699      elem[19][5] = 40.962278;  //41Ca (Isotope)
700      elem[19][6] = -1.;
701    
702      elem[20][0] = 44.955912;//Sc 21
703      elem[20][1] = -1.;
704      elem[20][2] = -1.;
705      elem[20][3] = -1.;
706      elem[20][4] = -1.;
707      elem[20][5] = -1.;
708      elem[20][6] = -1.;
709    
710      elem[21][0] = 47.947946;  //Ti 22
711      elem[21][1] = 45.952632;  //46Ti (Isotope)
712      elem[21][2] = 46.951763;  //47Ti (Isotope)
713      elem[21][3] = 48.947870;  //49Ti (Isotope)
714      elem[21][4] = 49.944791;  //50Ti (Isotope)
715      elem[21][5] = 43.959690;  //44Ti (Isotope) (half life 60y)
716      elem[21][6] = -1.;
717    
718      elem[22][0] = 50.943960; //V  23
719      elem[22][1] = 49.947158; //50V  (Isotope)
720      elem[22][2] = -1.;
721      elem[22][3] = -1.;
722      elem[22][4] = -1.;
723      elem[22][5] = -1.;
724      elem[22][6] = -1.;
725    
726      elem[23][0] = 51.940507; //Cr 24
727      elem[23][1] = 52.940649; //53Cr (Isotope)
728      elem[23][2] = 49.946044; //50Cr (Isotope)
729      elem[23][3] = 53.938880; //54Cr (Isotope)
730      elem[23][4] = -1.;
731      elem[23][5] = -1.;
732      elem[23][6] = -1.;
733    
734      elem[24][0] = 54.938049;//Mn 25
735      elem[24][1] = 52.941290;//53Mn (Isotope)
736      elem[24][2] = -1.;
737      elem[24][3] = -1.;
738      elem[24][4] = -1.;
739      elem[24][5] = -1.;
740      elem[24][6] = -1.;
741    
742      elem[25][0] = 55.934937;  //Fe 26
743      elem[25][1] = 53.939610;  //54Fe (Isotope)
744      elem[25][2] = 56.935394;  //57Fe (Isotope)
745      elem[25][3] = 57.933276;  //58Fe (Isotope)
746      elem[25][4] = 59.934072;  //58Fe (Isotope)
747    
748      elem[26][0] = 58.933195; //Co 27
749      elem[26][1] = 59.933817; //60Co (Isotope)
750      elem[26][2] = -1.;
751      elem[26][3] = -1.;
752      elem[26][4] = -1.;
753      elem[26][5] = -1.;
754      elem[26][6] = -1.;
755    
756    
757      elem[27][0] = 57.935343; //Ni 28
758      elem[27][1] = 61.928345; //62Ni (Isotope)
759      elem[27][2] = 59.930786; //60Ni (Isotope)
760      elem[27][3] = 60.931056; //61Ni (Isotope)
761      elem[27][4] = 63.927966; //64Ni (Isotope)
762      elem[27][5] = 58.934346; //59Ni (Isotope)
763      elem[27][6] = -1.;
764    
765      elem[28][0] = 62.929597; //Cu 29
766      elem[28][1] = 64.927789; //65Cu (Isotope)
767      elem[28][2] = -1.;
768      elem[28][3] = -1.;
769      elem[28][4] = -1.;
770      elem[28][5] = -1.;
771      elem[28][6] = -1.;
772    
773      elem[29][0] = 63.929142; //Zn 30
774      elem[29][1] = 65.926033; //66Zn (Isotope)
775      elem[29][2] = 67.924844; //68Zn (Isotope)
776      elem[29][3] = 66.927127; //67Zn (Isotope)
777      elem[29][4] = 69.925319; //70Zn (Isotope)
778      elem[29][5] = -1.;
779      elem[29][6] = -1.;
780    
781      elem[30][0] = 68.925573; //Ga 31
782      elem[30][1] = 70.924701; //71Ga (Isotope)
783      elem[30][2] = -1.;
784      elem[30][3] = -1.;
785      elem[30][4] = -1.;
786      elem[30][5] = -1.;
787      elem[30][6] = -1.;
788    
789      elem[31][0] = 73.921177; //Ge 32
790      elem[31][1] = 71.922075; //72Ge (Isotope)
791      elem[31][2] = 69.924247; //70Ge (Isotope)
792      elem[31][3] = 75.921403; //76Ge (Isotope)
793      elem[31][4] = 73.923459; //73Ge (Isotope)
794      elem[31][5] = -1.;
795      elem[31][6] = -1.;
796    
797    
798    //parametri calorimetro    //parametri calorimetro
# Line 552  void CaloBragg::conversione(){ Line 865  void CaloBragg::conversione(){
865    //gcm2W = WCross/10. * rhoW;    //gcm2W = WCross/10. * rhoW;
866    
867    //       (g/cm2W)/(g/cm2Si)    //       (g/cm2W)/(g/cm2Si)
868    if ( usenewBB ){    spessore[3] =  (WCross/10.) * rhoW;
869      spessore[1] =  (WCross/10.) * rhoW;    a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
870    } else {    spessore[1] =  a;
871      //riscala mip allo spessore attraversato
872      a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);    spessore[2] = MIP*(SiCross/ySi);  
     spessore[1] =  a;  
     //riscala mip allo spessore attraversato  
     spessore[2] = MIP*(SiCross/ySi);    
   };  
873  };//end conversione  };//end conversione
874    
875    
# Line 585  void CaloBragg::BetheBloch(Float_t *x, F Line 894  void CaloBragg::BetheBloch(Float_t *x, F
894    Float_t lg =0.;    Float_t lg =0.;
895    Float_t Energia=0.;    Float_t Energia=0.;
896    Float_t C=0.;    Float_t C=0.;
897    Float_t INo = II;    Float_t INo = ISi;
898    
899    if ( !usenewBB ) INo = ISi;    if ( usenewBB ) INo = II;
900    
901    eta = (*gam)*(*Bet);    eta = (*gam)*(*Bet);
902    
# Line 610  void CaloBragg::BetheBloch(Float_t *x, F Line 919  void CaloBragg::BetheBloch(Float_t *x, F
919    
920    
921    
922  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out, Float_t II){  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Int_t *isotope, Float_t *Etot, Float_t *out, Float_t II){
923        
924    /*perdita di energia per ioni pesanti (come da routine geant)*/    /*perdita di energia per ioni pesanti (come da routine geant)*/
925    //  in : dx    => spessore g/cm2    //  in : dx    => spessore g/cm2
# Line 627  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 936  void CaloBragg::ELOSS(Float_t *dx, Int_t
936    Float_t dEP=0.;    Float_t dEP=0.;
937        
938    // gamma  //  Mass = A * MassP; /*in Mev/c2*/    // gamma  //  Mass = A * MassP; /*in Mev/c2*/
939    gam =  (*Etot)/(elem[*Z-1]*MassP); // E = gamma M c2    gam =  (*Etot)/(elem[*Z-1][*isotope]*MassP); // E = gamma M c2
940        
941    
942    Bet = sqrt((SQ(gam) -1.)/SQ(gam));    Bet = sqrt((SQ(gam) -1.)/SQ(gam));
# Line 653  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 962  void CaloBragg::ELOSS(Float_t *dx, Int_t
962    
963    
964    
965  void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo,Float_t* ultimo, Float_t out[]){  void CaloBragg::Enetrack(Int_t* Z, Int_t* isotope, Float_t* E0, Float_t* primo,Float_t* ultimo, Float_t out[]){
966    
967    //calcola energia rilasciata sulla traccia (usa ELOSS)    //calcola energia rilasciata sulla traccia (usa ELOSS)
968    //  in : Z             =>carica    //  in : Z             =>carica
# Line 671  void CaloBragg::Enetrack(Int_t* Z, Float Line 980  void CaloBragg::Enetrack(Int_t* Z, Float
980    //azzero energia rilasciata sui piani    //azzero energia rilasciata sui piani
981    memset(out, 0, 2*NPLA*sizeof(Float_t));    memset(out, 0, 2*NPLA*sizeof(Float_t));
982        
983    Float_t Massa = (elem[(*Z)-1] * MassP);    Float_t Massa = (elem[(*Z)-1][*isotope] * MassP);
984    
985    for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){    for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
986      dE=0.;      dE=0.;
987      //spessore silicio corretto x inclinazione, z, energia, out:rilascio      //spessore silicio corretto x inclinazione, z, energia, out:rilascio
988      ELOSS(&spessore[0], Z, &Ezero, &dE, ISi);//spessore in g/cm2!!      ELOSS(&spessore[0], Z , isotope , &Ezero, &dE, ISi);//spessore in g/cm2!!
989      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop  
990            if(dE!=dE) return; //controlla che non sia un NaN
991        
992            if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
993        out[ipla] = Ezero - Massa; //MeV        out[ipla] = Ezero - Massa; //MeV
994        return;        return;
995                
996      }else{      }else{
997        out[ipla] = dE; //MeV        out[ipla] = dE; //MeV
998        Ezero = Ezero - dE;//energia residua        Ezero = Ezero - dE;//energia residua
999        if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero);        //      if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero);
1000      };      };
1001      //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
1002      if(ipla%2 == 0){      if(ipla%2 == 0){
# Line 693  void CaloBragg::Enetrack(Int_t* Z, Float Line 1005  void CaloBragg::Enetrack(Int_t* Z, Float
1005        Float_t sp = 0.;        Float_t sp = 0.;
1006        Float_t II = ISi;        Float_t II = ISi;
1007        if ( usenewBB ){        if ( usenewBB ){
1008          sp = spessore[1];                  sp = spessore[3];        
1009          II = IW;          II = IW;
1010        } else {        } else {
1011          sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2          sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2
1012        }        }
1013        ELOSS(&sp, Z, &Ezero, &dE,II);        //      printf(" sp %f II %f \n",sp,II);
1014          ELOSS(&sp, Z, isotope , &Ezero, &dE,II);
1015        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
1016          return;          return;
1017        }else{        }else{
# Line 773  void CaloBragg::chiquadro(Float_t dE[], Line 1086  void CaloBragg::chiquadro(Float_t dE[],
1086          //meno peso ai piani con rilasci maggiori di 1000 MIP          //meno peso ai piani con rilasci maggiori di 1000 MIP
1087          //      if(calorimetro[ipla][1] > 1000) wi=0.5;          //      if(calorimetro[ipla][1] > 1000) wi=0.5;
1088          if(calorimetro[ipla][1] > 1200.) wi=0.5;          if(calorimetro[ipla][1] > 1200.) wi=0.5;
1089                           if(debug) printf("chiquadro start \n ");
1090          Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));          Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
1091    
1092          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));
# Line 814  void CaloBragg::loopze( Float_t step, Fl Line 1127  void CaloBragg::loopze( Float_t step, Fl
1127    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2
1128    //    //
1129    
1130     //printf("entrato");
1131    memset(dEplan,0,2*NPLA*sizeof(Float_t));    memset(dEplan,0,2*NPLA*sizeof(Float_t));
1132    
1133    Int_t Z = 0;// z iniziale    Int_t Z = 0;// z iniziale
1134      
1135      Int_t isotope=0;
1136    
1137    Float_t Massa = 0.;    Float_t Massa = 0.;
1138    
# Line 831  void CaloBragg::loopze( Float_t step, Fl Line 1146  void CaloBragg::loopze( Float_t step, Fl
1146    Int_t max=32;//max z di cui so la massa :P    Int_t max=32;//max z di cui so la massa :P
1147    if((Zlimite)<=31) max=(int)(Zlimite) + 1;    if((Zlimite)<=31) max=(int)(Zlimite) + 1;
1148    
1149    if(debug)  printf("loopze inizio max  %d \n",max);
1150    if ( fzeta > 0. ){    if ( fzeta > 0. ){
1151      zmin = fzeta;      zmin = fzeta;
1152      max = fzeta+1;      max = fzeta+1;
# Line 838  void CaloBragg::loopze( Float_t step, Fl Line 1154  void CaloBragg::loopze( Float_t step, Fl
1154        
1155    Int_t colmax=32;    Int_t colmax=32;
1156    Int_t rowmax=3000;    Int_t rowmax=3000;
1157      Int_t isomax=7;
1158    
1159    Float_t matrixchi2[colmax][rowmax][3];    Float_t matrixchi2[colmax][isomax][rowmax][3];
1160    memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));    memset(matrixchi2, 0, colmax*isomax*rowmax*3*sizeof(Float_t));
1161    
1162    Int_t imin = 1-nostep/2;    Int_t imin = 1-nostep/2;
1163    Int_t imax = nostep/2;    Int_t imax = nostep/2;
# Line 849  void CaloBragg::loopze( Float_t step, Fl Line 1166  void CaloBragg::loopze( Float_t step, Fl
1166    for(Int_t inucl=zmin; inucl<max; inucl++){    for(Int_t inucl=zmin; inucl<max; inucl++){
1167            
1168      Z= inucl;      Z= inucl;
1169      
1170      Massa = elem[inucl-1]*MassP;    //loop isotopi
1171          while ( elem[inucl-1][isotope] > 0. ){  
1172    
1173              if( fiso != -1 ){
1174                      isotope=fiso;
1175                      if(debug) printf("In Loopze - Isotope N %d",isotope);
1176              }
1177        Massa = elem[inucl-1][isotope]*MassP;
1178    
1179      //loop energia      //loop energia
1180      Int_t iene2 = 0;      Int_t iene2 = -1;
1181    
1182      //    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 <=??      //    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 <=??
1183    
1184      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 <=??      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 <=??
1185             iene2++;
1186                     energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
1187            
1188            
1189                     if( fene > 0. ) energia=fene;   //forza l'energia
1190                     if (debug) printf("loopze energia %f, z %d, isotopo %d ,iene %d\n",energia,Z,isotope,iene);
1191    //              printf(" energia %f , forzata %f \n",energia,fene);
1192                     Enetrack(&Z, &isotope, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
1193                
1194        iene2++;           chiquadro(dEplan,chi2); //calcolo chi2
1195        energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..                   if (debug) printf("loopze chi %f \n",chi2[0]);
1196         if(debug && TMath::Finite(chi2[0])==1 && (TMath::IsNaN(chi2[0])!=1) )    printf("loopze    fin mat %f \n",chi2[0]);
       Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani  
       //calcolo chi2  
       chiquadro(dEplan,chi2);  
   
1197        //      printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]);        //      printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]);
1198                   if( (chi2[1] != 79.) ){//salto quelli che frammentano
1199        if( (chi2[1] != 79.) ){//salto quelli che frammentano                  matrixchi2[inucl][isotope][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia
1200          matrixchi2[inucl][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia                  matrixchi2[inucl][isotope][iene2][1]=energia;//energia per questo chi2
1201          matrixchi2[inucl][iene2][1]=energia;//energia per questo chi2                          matrixchi2[inucl][isotope][iene2][2]=chi2[2];//piani saltati nel chi2
1202          matrixchi2[inucl][iene2][2]=chi2[2];//piani saltati nel chi2              if( fene > 0. ) break;
1203        } else {             } else {      
1204          matrixchi2[inucl][iene2][0]=1000.;//valore chi2 per questo z a questa energia                          matrixchi2[inucl][isotope][iene2][0]=numeric_limits<Float_t>::max();//valore chi2 per questo z a questa energia
1205          matrixchi2[inucl][iene2][1]=1000.;//energia per questo chi2                          matrixchi2[inucl][isotope][iene2][1]=numeric_limits<Float_t>::max();//energia per questo chi2
1206          matrixchi2[inucl][iene2][2]=1000.;//piani saltati nel chi2                          matrixchi2[inucl][isotope][iene2][2]=numeric_limits<Float_t>::max();//piani saltati nel chi2
1207          break;                  break;
1208        }        }
     }//fine loop energia  
1209    
1210        }//fine loop energia
1211    
1212    };//fine loop z    if( fiso != -1 ){
1213              if(debug) printf("exited form isotopes loop");
1214              break;
1215      }
1216    
1217      isotope++; //incremento il contatore isotopi
1218      }//fine loop isotopi
1219      isotope=0; //riazzero il contatore isotopi
1220    
1221      }//fine loop z
1222    
1223      isotope=0;//non dovrebbe servire
1224    
1225    //Emi    //Emi
1226    for (Int_t nu=zmin; nu<max; nu++){   for (Int_t nu=zmin; nu<max; nu++){
1227      for (Int_t en=0; en<nostep; en++){     if( fiso != -1 ){
1228        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){            isotope=fiso;
1229          bestchi2[0]= matrixchi2[nu][en][0];// chi2            if(debug) printf("In Loopze EMI - Isotope N %d",isotope);
1230          bestchi2[1]= (Float_t)nu; // z     }
1231          bestchi2[2]= matrixchi2[nu][en][1];//energia;       while(elem[nu-1][isotope]> 0.){
1232          bestchi2[3]= matrixchi2[nu][en][2];// totale piani saltati          for (Int_t en=0; en<nostep; en++){
1233        }          if((matrixchi2[nu][isotope][en][0]<bestchi2[0]) && (matrixchi2[nu][isotope][en][0] >0.)){
1234      }            bestchi2[0]= matrixchi2[nu][isotope][en][0];// chi2
1235              bestchi2[1]= (Float_t)nu; // z
1236              bestchi2[2]= matrixchi2[nu][isotope][en][1];//energia;
1237              bestchi2[3]= matrixchi2[nu][isotope][en][2];// totale piani saltati
1238          bestchi2[4]= (Float_t)isotope; //isotopo
1239            }
1240         }
1241    
1242      if( fiso != -1 ){
1243        if(debug) printf("exited form isotopes loop");
1244        break;
1245    }    }
1246    
1247        isotope++;
1248        }
1249       isotope=0;
1250     }
1251    
1252  };//endloopze  };//endloopze
1253    
# Line 1017  void CaloBragg::Zdaloop(){ Line 1371  void CaloBragg::Zdaloop(){
1371    
1372    Float_t zmin=1.;    Float_t zmin=1.;
1373    Float_t zmax=32.;    Float_t zmax=32.;
1374    Float_t bestchitemp[4] = {0,0,0,0};    Float_t bestchitemp[5] = {0,0,0,0,0};
1375    
1376    bestchi2[0]=10000.;    bestchi2[0]=numeric_limits<Float_t>::max();
1377    bestchi2[1]=0.;    bestchi2[1]=0.;
1378    bestchi2[2]=0.;    bestchi2[2]=0.;
1379    bestchi2[3]=0.;    bestchi2[3]=0.;
1380    Float_t zero=0.;    bestchi2[4]=0.;
1381      //  Float_t zero=0.;
1382    //------------primo loop   ----------------------    //------------primo loop   ----------------------
1383    //     energia   ezero, zstart  zstop    //     energia   ezero, zstart  zstop
1384    //  loopze(Integrale,zero,zmin,zmax);    //  loopze(Integrale,zero,zmin,zmax);
# Line 1032  void CaloBragg::Zdaloop(){ Line 1387  void CaloBragg::Zdaloop(){
1387    loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200);    loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200);
1388    
1389    //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);    //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);
1390    if ( debug ) printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]);    if ( debug) printf("Zdaloop start Integrale %f , outene %f \n",Integrale,bestchi2[2]);
1391        
1392    //------------secondo loop  ----------------------    //------------secondo loop  ----------------------
1393    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];    for(Int_t i=0;i<5;i++) bestchitemp[i]=bestchi2[i];
1394    bestchi2[0] = 10000.;    bestchi2[0]=numeric_limits<Float_t>::max();
1395    bestchi2[1] = 0.;    bestchi2[1] = 0.;
1396    bestchi2[2] = 0.;    bestchi2[2] = 0.;
1397    bestchi2[3] = 0.;//riazzero    bestchi2[3] = 0.;
1398      bestchi2[4] = 0.;//riazzero
1399    
1400    Float_t step = bestchitemp[2];//    Float_t step = bestchitemp[2];//
1401    zero=0.;  // qualsiasi altro valore peggiora le cose    //  zero=0.;  // qualsiasi altro valore peggiora le cose
1402    //  zmin=zmax=bestchitemp[1];    //  zmin=zmax=bestchitemp[1];
1403    zmin=bestchitemp[1]-1;    zmin=bestchitemp[1]-1;
1404    zmax=bestchitemp[1]+1;    zmax=bestchitemp[1]+1;
1405    //  loopze(step,zero,zmin,zmax); //    //loopze(step,zero,zmin,zmax); //
1406    
1407    //->  loopze(step,step/2.,zmin,zmax,200); //    //->  loopze(step,step/2.,zmin,zmax,200); //
1408    loopze(step,step/2.,zmin,zmax,500); //    loopze(step,step/2.,zmin,zmax,500); //
1409    
1410    if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step);    //step = bestchitemp[2];//
1411    
1412      //loopze(step/2,3*step/4.,zmin,zmax,500); //
1413    
1414      if ( debug ) printf("Zdaloop Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step);
1415    
     
1416    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
1417    lpchi2=bestchi2[0];    lpchi2=bestchi2[0];
1418    lpz=bestchi2[1];    lpz=bestchi2[1];
1419    lpetot=bestchi2[2];    lpetot=bestchi2[2];
1420    lppskip=bestchi2[3];    lppskip=bestchi2[3];
1421      lpisotope=bestchi2[4];
1422  };//endZdaloop  };//endZdaloop
1423    
1424    
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.23