/[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.12 by mocchiut, Fri Jun 17 11:56:42 2011 UTC revision 1.13 by mocchiut, Thu Oct 6 09:27:39 2011 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,4*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));
# 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 343  Float_t CaloBragg::Integral(){ Line 346  Float_t CaloBragg::Integral(){
346    
347    Float_t dEpianiloop[44];    Float_t dEpianiloop[44];
348    Int_t tz1=(Int_t)lpz;    Int_t tz1=(Int_t)lpz;
349    Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop    Int_t ti1=(Int_t)lpisotope;
350    
351      Enetrack(&tz1, &ti1 , &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
352    
353    
354    Float_t integ = 0.;    Float_t integ = 0.;
# Line 372  void CaloBragg::Draw(){ Line 377  void CaloBragg::Draw(){
377    
378    Process();    Process();
379    
380    this->Draw(0.,0.);    this->Draw(0.,0.,0.);
381    
382  }  }
383    
384  void CaloBragg::Draw(Int_t Z, Float_t enetot){  void CaloBragg::Draw(Int_t Z, Int_t isotope, Float_t enetot){
385    
386    //  Float_t dEpianimean[44];    //  Float_t dEpianimean[44];
387    Float_t dEpianiloop[44];    Float_t dEpianiloop[44];
388    Float_t Depth[44];    Float_t Depth[44];
389    //  Int_t tz=(Int_t)qtz;    //  Int_t tz=(Int_t)qtz;
390    Int_t tz1= Z;    Int_t tz1= Z;
391      Int_t ti1= isotope;
392    Float_t enet = enetot;    Float_t enet = enetot;
393    //  Float_t enet = lpetot;    //  Float_t enet = lpetot;
394        
# Line 414  void CaloBragg::Draw(Int_t Z, Float_t en Line 420  void CaloBragg::Draw(Int_t Z, Float_t en
420    
421    } else {    } else {
422      tz1=(Int_t)lpz;      tz1=(Int_t)lpz;
423        ti1=(Int_t)lpisotope;
424      enet = lpetot;      enet = lpetot;
425      //  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
426        
427    }    }
428    Enetrack(&tz1, &enet, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop    Enetrack(&tz1, &ti1, &enet, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
429    
430    Float_t sp= spessore[0]*spessore[1];    Float_t sp= spessore[0]*spessore[1];
431    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 491  void CaloBragg::Draw(Int_t Z, Float_t en Line 498  void CaloBragg::Draw(Int_t Z, Float_t en
498    
499  void CaloBragg::LoadParam(){  void CaloBragg::LoadParam(){
500    
501    //    // elem[Z-1][isotop] 0 is the most common one
502    elem[0] = 1.00794; //H  1    //
503    elem[1] = 4.002602;  //He 2  
504    elem[2] = 6.941;   //Li 3    elem[0][0] = 1.00782;   //H  1
505    elem[3] = 9.012182;//Be 4    elem[0][1] = 2.01410;   // 2H (Isotope)
506    elem[4] = 10.811;  //B  5    elem[0][2] = -1.;
507    elem[5] = 12.0107; //C  6    elem[0][3] = -1.;
508    elem[6] = 14.00674;//N  7    elem[0][4] = -1.;
509    elem[7] = 15.9994; //O  8    elem[0][5] = -1.;
510    elem[8] = 18.9984032; //F  9    elem[0][6] = -1.;
511    elem[9] = 20.1797; //Ne 10  
512    elem[10] = 22.98977;//Na 11    elem[1][0] = 4.002602;  //He 2
513    elem[11] = 24.3050; //Mg 12    elem[1][1] = 3.016029;  // 3He (Isotope)
514    elem[12] = 26.9815; //Al 13    elem[1][2] = -1.;
515    elem[13] = 28.0855; //Si 14    elem[1][3] = -1.;
516    elem[14] = 30.974;  //P  15    elem[1][4] = -1.;
517    elem[15] = 32.066;  //S  16    elem[1][5] = -1.;
518    elem[16] = 35.4527; //Cl 17    elem[1][6] = -1.;
519    elem[17] = 39.948;  //Ar 18  
520    elem[18] = 39.0983; //K  19    elem[2][0] = 7.016004;  //Li 3
521    elem[19] = 40.078;  //Ca 20    elem[2][1] = 6.015123;  //6Li (Isotope)
522    elem[20] = 44.95591;//Sc 21    elem[2][2] = -1.;
523    elem[21] = 47.867;  //Ti 22    elem[2][3] = -1.;
524    elem[22] = 50.9415; //V  23    elem[2][4] = -1.;
525    elem[23] = 51.9961; //Cr 24    elem[2][5] = -1.;
526    elem[24] = 54.938049;//Mn 25    elem[2][6] = -1.;
527    elem[25] = 55.845;  //Fe 26  
528    elem[26] = 58.9332; //Co 27    elem[3][0] = 9.012182;  //Be 4
529    elem[27] = 58.6934; //Ni 28    elem[3][1] = 10.01353;  //10Be (Isotope) (most stable)
530    elem[28] = 63.546; //Cu 29    elem[3][2] = -1.;
531    elem[29] = 65.39; //Zn 30    elem[3][3] = -1.;
532    elem[30] = 69.723; //Ga 31    elem[3][4] = -1.;
533    elem[31] = 72.61; //Ge 32    elem[3][5] = -1.;
534      elem[3][6] = -1.;
535    
536      elem[4][0] = 11.0093;   //B 5
537      elem[4][1] = 10.01294;  //10B  (Isotope)
538      elem[4][2] = -1.;
539      elem[4][3] = -1.;
540      elem[4][4] = -1.;
541      elem[4][5] = -1.;
542      elem[4][5] = -1.;
543    
544      elem[5][0] = 12.0107;   //C  6
545      elem[5][1] = 13.00335;  //13C (Isotope)
546      elem[5][2] = -1.;
547      elem[5][3] = -1.;
548      elem[5][4] = -1.;
549      elem[5][5] = -1.;
550      elem[5][5] = -1.;
551    
552      elem[6][0] = 14.00674;  //N  7
553      elem[6][1] = 15.00011;  //15N  (Isotope)
554      elem[6][2] = -1.;
555      elem[6][3] = -1.;
556      elem[6][4] = -1.;
557      elem[6][5] = -1.;
558      elem[6][5] = -1.;
559    
560      elem[7][0] = 15.99491;  //O  8
561      elem[7][1] = 17.99916;  //18O  (Isotope)
562      elem[7][2] = 16.99916;  //17O  (Isotope)
563      elem[7][3] = -1.;
564      elem[7][4] = -1.;
565      elem[7][5] = -1.;
566      elem[7][5] = -1.;
567    
568      elem[8][0] = 18.99840;  //F  9
569      elem[8][1] = -1.;
570      elem[8][2] = -1.;
571      elem[8][3] = -1.;
572      elem[8][4] = -1.;
573      elem[8][5] = -1.;
574      elem[8][5] = -1.;
575    
576      elem[9][0] = 19.99244;  //Ne 10
577      elem[9][1] = 21.99138;  //22Ne (Isotope)
578      elem[9][2] = 20.99384;  //21Ne 10
579      elem[9][3] = -1.;
580      elem[9][4] = -1.;
581      elem[9][5] = -1.;
582      elem[9][6] = -1.;
583    
584      elem[10][0] = 22.98977; //Na 11
585      elem[10][1] = 21.99444; //22Na (Isotope) (most stable)
586      elem[10][2] = -1.;
587      elem[10][3] = -1.;
588      elem[10][4] = -1.;
589      elem[10][5] = -1.;
590      elem[10][6] = -1.;
591    
592      elem[11][0] = 23.98504; //Mg 12
593      elem[11][1] = 25.98259; //26Mg (Isotope)
594      elem[11][2] = 24.98504; //25Mg (Isotope)
595      elem[11][3] = -1.;
596      elem[11][4] = -1.;
597      elem[11][5] = -1.;
598      elem[11][6] = -1.;
599    
600      elem[12][0] = 26.98154; //Al 13
601      elem[12][1] = 25.98489; //26Al (Isotope) (most stable)
602      elem[12][2] = -1.;
603      elem[12][3] = -1.;
604      elem[12][4] = -1.;
605      elem[12][5] = -1.;
606      elem[12][6] = -1.;
607    
608      elem[13][0] = 27.97692; //Si 14
609      elem[13][1] = 28.97649; //29Si (Isotope)
610      elem[13][2] = 29.97377; //30Si (Isotope)
611      elem[13][3] = -1.;
612      elem[13][4] = -1.;
613      elem[13][5] = -1.;
614      elem[13][6] = -1.;
615    
616      elem[14][0] = 30.97376;  //P  15
617      elem[14][1] = -1.;
618      elem[14][2] = -1.;
619      elem[14][3] = -1.;
620      elem[14][4] = -1.;
621      elem[14][5] = -1.;
622      elem[14][6] = -1.;
623    
624      elem[15][0] = 31.97207;  //S  16
625      elem[15][1] = 33.96787;  //34S  (Isotope)
626      elem[15][2] = 32.97146;  //33S  (Isotope)
627      elem[15][3] = 35.96708;  //36S  (Isotope)
628      elem[15][4] = -1.;
629      elem[15][5] = -1.;
630      elem[15][6] = -1.;
631    
632      elem[16][0] = 34.96885; //Cl 17
633      elem[16][1] = 36.96831; //37Cl 17
634      elem[16][2] = 35.96890; //36Cl (Isotope)
635      elem[16][3] = -1.;
636      elem[16][4] = -1.;
637      elem[16][5] = -1.;
638      elem[16][6] = -1.;
639    
640      elem[17][0] = 39.962383;  //Ar 18
641      elem[17][1] = 35.967545;  //36Ar (Isotope)
642      elem[17][2] = 37.962732;  //38Ar (Isotope)
643      elem[17][3] = 38.964313;  //39Ar (Isotope)
644      elem[17][4] = -1.;
645      elem[17][5] = -1.;
646      elem[17][6] = -1.;
647    
648      elem[18][0] = 38.963707; //K  19
649      elem[18][1] = 40.961825; //41K  (Isotope)
650      elem[18][2] = 39.963998; //40K  (Isotope)
651      elem[18][3] = -1.;
652      elem[18][4] = -1.;
653      elem[18][5] = -1.;
654      elem[18][6] = -1.;
655    
656      elem[19][0] = 39.962590;  //Ca 20
657      elem[19][1] = 43.955482;  //44Ca (Isotope)
658      elem[19][2] = 41.958618;  //42Ca (Isotope)
659      elem[19][3] = 42.958767;  //43Ca (Isotope)
660      elem[19][4] = 45.953693;  //46Ca (Isotope)
661      elem[19][5] = 40.962278;  //41Ca (Isotope)
662      elem[19][6] = -1.;
663    
664      elem[20][0] = 44.955912;//Sc 21
665      elem[20][1] = -1.;
666      elem[20][2] = -1.;
667      elem[20][3] = -1.;
668      elem[20][4] = -1.;
669      elem[20][5] = -1.;
670      elem[20][6] = -1.;
671    
672      elem[21][0] = 47.947946;  //Ti 22
673      elem[21][1] = 45.952632;  //46Ti (Isotope)
674      elem[21][2] = 46.951763;  //47Ti (Isotope)
675      elem[21][3] = 48.947870;  //49Ti (Isotope)
676      elem[21][4] = 49.944791;  //50Ti (Isotope)
677      elem[21][5] = 43.959690;  //44Ti (Isotope) (half life 60y)
678      elem[21][6] = -1.;
679    
680      elem[22][0] = 50.943960; //V  23
681      elem[22][1] = 49.947158; //50V  (Isotope)
682      elem[22][2] = -1.;
683      elem[22][3] = -1.;
684      elem[22][4] = -1.;
685      elem[22][5] = -1.;
686      elem[22][6] = -1.;
687    
688      elem[23][0] = 51.940507; //Cr 24
689      elem[23][1] = 52.940649; //53Cr (Isotope)
690      elem[23][2] = 49.946044; //50Cr (Isotope)
691      elem[23][3] = 53.938880; //54Cr (Isotope)
692      elem[23][4] = -1.;
693      elem[23][5] = -1.;
694      elem[23][6] = -1.;
695    
696      elem[24][0] = 54.938049;//Mn 25
697      elem[24][1] = 52.941290;//53Mn (Isotope)
698      elem[24][2] = -1.;
699      elem[24][3] = -1.;
700      elem[24][4] = -1.;
701      elem[24][5] = -1.;
702      elem[24][6] = -1.;
703    
704      elem[25][0] = 55.934937;  //Fe 26
705      elem[25][1] = 53.939610;  //54Fe (Isotope)
706      elem[25][2] = 56.935394;  //57Fe (Isotope)
707      elem[25][3] = 57.933276;  //58Fe (Isotope)
708      elem[25][4] = 59.934072;  //58Fe (Isotope)
709    
710      elem[26][0] = 58.933195; //Co 27
711      elem[26][1] = 59.933817; //60Co (Isotope)
712      elem[26][2] = -1.;
713      elem[26][3] = -1.;
714      elem[26][4] = -1.;
715      elem[26][5] = -1.;
716      elem[26][6] = -1.;
717    
718    
719      elem[27][0] = 57.935343; //Ni 28
720      elem[27][1] = 61.928345; //62Ni (Isotope)
721      elem[27][2] = 59.930786; //60Ni (Isotope)
722      elem[27][3] = 60.931056; //61Ni (Isotope)
723      elem[27][4] = 63.927966; //64Ni (Isotope)
724      elem[27][5] = 58.934346; //59Ni (Isotope)
725      elem[27][6] = -1.;
726    
727      elem[28][0] = 62.929597; //Cu 29
728      elem[28][1] = 64.927789; //65Cu (Isotope)
729      elem[28][2] = -1.;
730      elem[28][3] = -1.;
731      elem[28][4] = -1.;
732      elem[28][5] = -1.;
733      elem[28][6] = -1.;
734    
735      elem[29][0] = 63.929142; //Zn 30
736      elem[29][1] = 65.926033; //66Zn (Isotope)
737      elem[29][2] = 67.924844; //68Zn (Isotope)
738      elem[29][3] = 66.927127; //67Zn (Isotope)
739      elem[29][4] = 69.925319; //70Zn (Isotope)
740      elem[29][5] = -1.;
741      elem[29][6] = -1.;
742    
743      elem[30][0] = 68.925573; //Ga 31
744      elem[30][1] = 70.924701; //71Ga (Isotope)
745      elem[30][2] = -1.;
746      elem[30][3] = -1.;
747      elem[30][4] = -1.;
748      elem[30][5] = -1.;
749      elem[30][6] = -1.;
750    
751      elem[31][0] = 73.921177; //Ge 32
752      elem[31][1] = 71.922075; //72Ge (Isotope)
753      elem[31][2] = 69.924247; //70Ge (Isotope)
754      elem[31][3] = 75.921403; //76Ge (Isotope)
755      elem[31][4] = 73.923459; //73Ge (Isotope)
756      elem[31][5] = -1.;
757      elem[31][6] = -1.;
758    
759    
760    //parametri calorimetro    //parametri calorimetro
# Line 650  void CaloBragg::BetheBloch(Float_t *x, F Line 881  void CaloBragg::BetheBloch(Float_t *x, F
881    
882    
883    
884  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){
885        
886    /*perdita di energia per ioni pesanti (come da routine geant)*/    /*perdita di energia per ioni pesanti (come da routine geant)*/
887    //  in : dx    => spessore g/cm2    //  in : dx    => spessore g/cm2
# Line 667  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 898  void CaloBragg::ELOSS(Float_t *dx, Int_t
898    Float_t dEP=0.;    Float_t dEP=0.;
899        
900    // gamma  //  Mass = A * MassP; /*in Mev/c2*/    // gamma  //  Mass = A * MassP; /*in Mev/c2*/
901    gam =  (*Etot)/(elem[*Z-1]*MassP); // E = gamma M c2    gam =  (*Etot)/(elem[*Z-1][*isotope]*MassP); // E = gamma M c2
902        
903    
904    Bet = sqrt((SQ(gam) -1.)/SQ(gam));    Bet = sqrt((SQ(gam) -1.)/SQ(gam));
# Line 693  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 924  void CaloBragg::ELOSS(Float_t *dx, Int_t
924    
925    
926    
927  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[]){
928    
929    //calcola energia rilasciata sulla traccia (usa ELOSS)    //calcola energia rilasciata sulla traccia (usa ELOSS)
930    //  in : Z             =>carica    //  in : Z             =>carica
# Line 711  void CaloBragg::Enetrack(Int_t* Z, Float Line 942  void CaloBragg::Enetrack(Int_t* Z, Float
942    //azzero energia rilasciata sui piani    //azzero energia rilasciata sui piani
943    memset(out, 0, 2*NPLA*sizeof(Float_t));    memset(out, 0, 2*NPLA*sizeof(Float_t));
944        
945    Float_t Massa = (elem[(*Z)-1] * MassP);    Float_t Massa = (elem[(*Z)-1][*isotope] * MassP);
946    
947    for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){    for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
948      dE=0.;      dE=0.;
949      //spessore silicio corretto x inclinazione, z, energia, out:rilascio      //spessore silicio corretto x inclinazione, z, energia, out:rilascio
950      ELOSS(&spessore[0], Z, &Ezero, &dE, ISi);//spessore in g/cm2!!      ELOSS(&spessore[0], Z , isotope , &Ezero, &dE, ISi);//spessore in g/cm2!!
951      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop  
952            if(dE!=dE) return; //controlla che non sia un NaN
953        
954            if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
955        out[ipla] = Ezero - Massa; //MeV        out[ipla] = Ezero - Massa; //MeV
956        return;        return;
957                
# Line 739  void CaloBragg::Enetrack(Int_t* Z, Float Line 973  void CaloBragg::Enetrack(Int_t* Z, Float
973          sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2          sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2
974        }        }
975        //      printf(" sp %f II %f \n",sp,II);        //      printf(" sp %f II %f \n",sp,II);
976        ELOSS(&sp, Z, &Ezero, &dE,II);        ELOSS(&sp, Z, isotope , &Ezero, &dE,II);
977        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
978          return;          return;
979        }else{        }else{
# Line 814  void CaloBragg::chiquadro(Float_t dE[], Line 1048  void CaloBragg::chiquadro(Float_t dE[],
1048          //meno peso ai piani con rilasci maggiori di 1000 MIP          //meno peso ai piani con rilasci maggiori di 1000 MIP
1049          //      if(calorimetro[ipla][1] > 1000) wi=0.5;          //      if(calorimetro[ipla][1] > 1000) wi=0.5;
1050          if(calorimetro[ipla][1] > 1200.) wi=0.5;          if(calorimetro[ipla][1] > 1200.) wi=0.5;
1051                           if(debug) printf("chiquadro start \n ");
1052          Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));          Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
1053    
1054          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 855  void CaloBragg::loopze( Float_t step, Fl Line 1089  void CaloBragg::loopze( Float_t step, Fl
1089    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2
1090    //    //
1091    
1092     //printf("entrato");
1093    memset(dEplan,0,2*NPLA*sizeof(Float_t));    memset(dEplan,0,2*NPLA*sizeof(Float_t));
1094    
1095    Int_t Z = 0;// z iniziale    Int_t Z = 0;// z iniziale
1096      
1097      Int_t isotope=0;
1098    
1099    Float_t Massa = 0.;    Float_t Massa = 0.;
1100    
# Line 872  void CaloBragg::loopze( Float_t step, Fl Line 1108  void CaloBragg::loopze( Float_t step, Fl
1108    Int_t max=32;//max z di cui so la massa :P    Int_t max=32;//max z di cui so la massa :P
1109    if((Zlimite)<=31) max=(int)(Zlimite) + 1;    if((Zlimite)<=31) max=(int)(Zlimite) + 1;
1110    
1111    if(debug)  printf("loopze inizio max  %d \n",max);
1112    if ( fzeta > 0. ){    if ( fzeta > 0. ){
1113      zmin = fzeta;      zmin = fzeta;
1114      max = fzeta+1;      max = fzeta+1;
# Line 879  void CaloBragg::loopze( Float_t step, Fl Line 1116  void CaloBragg::loopze( Float_t step, Fl
1116        
1117    Int_t colmax=32;    Int_t colmax=32;
1118    Int_t rowmax=3000;    Int_t rowmax=3000;
1119      Int_t isomax=7;
1120    
1121    Float_t matrixchi2[colmax][rowmax][3];    Float_t matrixchi2[colmax][isomax][rowmax][3];
1122    memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));    memset(matrixchi2, 0, colmax*isomax*rowmax*3*sizeof(Float_t));
1123    
1124    Int_t imin = 1-nostep/2;    Int_t imin = 1-nostep/2;
1125    Int_t imax = nostep/2;    Int_t imax = nostep/2;
# Line 890  void CaloBragg::loopze( Float_t step, Fl Line 1128  void CaloBragg::loopze( Float_t step, Fl
1128    for(Int_t inucl=zmin; inucl<max; inucl++){    for(Int_t inucl=zmin; inucl<max; inucl++){
1129            
1130      Z= inucl;      Z= inucl;
1131      
1132      Massa = elem[inucl-1]*MassP;    //loop isotopi
1133          while ( elem[inucl-1][isotope] > 0. ){  
1134        Massa = elem[inucl-1][isotope]*MassP;
1135    
1136      //loop energia      //loop energia
1137      Int_t iene2 = 0;      Int_t iene2 = -1;
1138    
1139      //    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 <=??
1140    
1141      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 <=??
1142             iene2++;
1143                     energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
1144            
1145            
1146                     if( fene > 0. ) energia=fene;   //forza l'energia
1147                     if (debug) printf("loopze energia %f, z %d, isotopo %d ,iene %d\n",energia,Z,isotope,iene);
1148    //              printf(" energia %f , forzata %f \n",energia,fene);
1149                     Enetrack(&Z, &isotope, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
1150                
1151        iene2++;           chiquadro(dEplan,chi2); //calcolo chi2
1152        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]);
1153         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);  
   
1154        //      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]);
1155                   if( (chi2[1] != 79.) ){//salto quelli che frammentano
1156        if( (chi2[1] != 79.) ){//salto quelli che frammentano                  matrixchi2[inucl][isotope][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia
1157          matrixchi2[inucl][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia                  matrixchi2[inucl][isotope][iene2][1]=energia;//energia per questo chi2
1158          matrixchi2[inucl][iene2][1]=energia;//energia per questo chi2                          matrixchi2[inucl][isotope][iene2][2]=chi2[2];//piani saltati nel chi2
1159          matrixchi2[inucl][iene2][2]=chi2[2];//piani saltati nel chi2              if( fene > 0. ) break;
1160        } else {             } else {      
1161          matrixchi2[inucl][iene2][0]=1000.;//valore chi2 per questo z a questa energia                          matrixchi2[inucl][isotope][iene2][0]=1000.;//valore chi2 per questo z a questa energia
1162          matrixchi2[inucl][iene2][1]=1000.;//energia per questo chi2                          matrixchi2[inucl][isotope][iene2][1]=1000.;//energia per questo chi2
1163          matrixchi2[inucl][iene2][2]=1000.;//piani saltati nel chi2                          matrixchi2[inucl][isotope][iene2][2]=1000.;//piani saltati nel chi2
1164          break;                  break;
1165        }        }
     }//fine loop energia  
1166    
1167        }//fine loop energia
1168    
1169      isotope++; //incremento il contatore isotopi
1170      }//fine loop isotopi
1171      isotope=0; //riazzero il contatore isotopi
1172    
1173    };//fine loop z    }//fine loop z
1174    
1175      isotope=0;//non dovrebbe servire
1176    
1177    //Emi    //Emi
1178    for (Int_t nu=zmin; nu<max; nu++){    for (Int_t nu=zmin; nu<max; nu++){
1179      for (Int_t en=0; en<nostep; en++){      while(elem[nu-1][isotope]> 0.){
1180        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){          for (Int_t en=0; en<nostep; en++){
1181          bestchi2[0]= matrixchi2[nu][en][0];// chi2          if((matrixchi2[nu][isotope][en][0]<bestchi2[0]) && (matrixchi2[nu][isotope][en][0] >0.)){
1182          bestchi2[1]= (Float_t)nu; // z            bestchi2[0]= matrixchi2[nu][isotope][en][0];// chi2
1183          bestchi2[2]= matrixchi2[nu][en][1];//energia;            bestchi2[1]= (Float_t)nu; // z
1184          bestchi2[3]= matrixchi2[nu][en][2];// totale piani saltati            bestchi2[2]= matrixchi2[nu][isotope][en][1];//energia;
1185              bestchi2[3]= matrixchi2[nu][isotope][en][2];// totale piani saltati
1186          bestchi2[4]= (Float_t)isotope; //isotopo
1187            }
1188        }        }
1189        isotope++;
1190      }      }
1191       isotope=0;
1192    }    }
1193    
   
1194  };//endloopze  };//endloopze
1195    
1196    
# Line 1058  void CaloBragg::Zdaloop(){ Line 1313  void CaloBragg::Zdaloop(){
1313    
1314    Float_t zmin=1.;    Float_t zmin=1.;
1315    Float_t zmax=32.;    Float_t zmax=32.;
1316    Float_t bestchitemp[4] = {0,0,0,0};    Float_t bestchitemp[5] = {0,0,0,0,0};
1317    
1318    bestchi2[0]=10000.;    bestchi2[0]=numeric_limits<Float_t>::max();
1319    bestchi2[1]=0.;    bestchi2[1]=0.;
1320    bestchi2[2]=0.;    bestchi2[2]=0.;
1321    bestchi2[3]=0.;    bestchi2[3]=0.;
1322      bestchi2[4]=0.;
1323    Float_t zero=0.;    Float_t zero=0.;
1324    //------------primo loop   ----------------------    //------------primo loop   ----------------------
1325    //     energia   ezero, zstart  zstop    //     energia   ezero, zstart  zstop
# Line 1073  void CaloBragg::Zdaloop(){ Line 1329  void CaloBragg::Zdaloop(){
1329    loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200);    loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200);
1330    
1331    //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);    //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);
1332    if ( debug ) printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]);    if ( debug) printf("Zdaloop start Integrale %f , outene %f \n",Integrale,bestchi2[2]);
1333        
1334    //------------secondo loop  ----------------------    //------------secondo loop  ----------------------
1335    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];    for(Int_t i=0;i<5;i++) bestchitemp[i]=bestchi2[i];
1336    bestchi2[0] = 10000.;    bestchi2[0]=numeric_limits<Float_t>::max();
1337    bestchi2[1] = 0.;    bestchi2[1] = 0.;
1338    bestchi2[2] = 0.;    bestchi2[2] = 0.;
1339    bestchi2[3] = 0.;//riazzero    bestchi2[3] = 0.;
1340      bestchi2[4] = 0.;//riazzero
1341    
1342    Float_t step = bestchitemp[2];//    Float_t step = bestchitemp[2];//
1343    zero=0.;  // qualsiasi altro valore peggiora le cose    zero=0.;  // qualsiasi altro valore peggiora le cose
1344    //  zmin=zmax=bestchitemp[1];    //  zmin=zmax=bestchitemp[1];
1345    zmin=bestchitemp[1]-1;    zmin=bestchitemp[1]-1;
1346    zmax=bestchitemp[1]+1;    zmax=bestchitemp[1]+1;
1347    //  loopze(step,zero,zmin,zmax); //    //loopze(step,zero,zmin,zmax); //
1348    
1349    //->  loopze(step,step/2.,zmin,zmax,200); //    //->  loopze(step,step/2.,zmin,zmax,200); //
1350    loopze(step,step/2.,zmin,zmax,500); //    loopze(step,step/2.,zmin,zmax,500); //
1351    
1352    if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step);    //step = bestchitemp[2];//
1353    
1354      //loopze(step/2,3*step/4.,zmin,zmax,500); //
1355    
1356      if ( debug ) printf("Zdaloop Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step);
1357    
     
1358    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
1359    lpchi2=bestchi2[0];    lpchi2=bestchi2[0];
1360    lpz=bestchi2[1];    lpz=bestchi2[1];
1361    lpetot=bestchi2[2];    lpetot=bestchi2[2];
1362    lppskip=bestchi2[3];    lppskip=bestchi2[3];
1363      lpisotope=bestchi2[4];
1364  };//endZdaloop  };//endZdaloop
1365    
1366    
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.23