--- calo/flight/CaloBragg/src/CaloBragg.cpp 2011/06/17 11:56:42 1.12 +++ calo/flight/CaloBragg/src/CaloBragg.cpp 2011/11/23 14:33:32 1.15 @@ -43,8 +43,10 @@ // qtpskip = 0.; lpchi2 = 0.; lpz = 0.; + lpisotope= 0.; lpetot = 0.; lppskip = 0.; + memset(calorimetro,0,44*2*sizeof(Float_t)); memset(spessore,0,4*sizeof(Float_t)); memset(estremi,0,2*2*sizeof(Float_t)); @@ -71,6 +73,7 @@ // printf(" plane not used for truncated mean %f: \n", qtpskip); printf(" chi 2 from loop %f: \n", lpchi2); printf(" Z from loop %f: \n", lpz); + printf(" isotope from loop %f: \n", lpisotope); printf(" energy from loop %f: \n", lpetot); printf(" plane not used for loop %f: \n", lppskip); printf(" ndf: %i \n",ndf); @@ -89,29 +92,44 @@ }; -void CaloBragg::CleanPlanes(Float_t epiano[22][2]){ +void CaloBragg::CleanPlanes(Float_t epiano[22][2], Bool_t zpiano[22][2]){ // return; Int_t hitplanes = 0; + Float_t f5 = 0.; for (Int_t i = 0; i<22; i++){ for (Int_t j = 1; j>=0; j--){ - if ( epiano[i][j] > 0.7 ) hitplanes++; + zpiano[i][j] = false; + if ( epiano[i][j] > 0.7 ){ + if ( hitplanes < 100 ) f5 += epiano[i][j]; + hitplanes++; + }; }; }; + Int_t atl5 = TMath::Min(hitplanes,100); + atl5 = TMath::Max(atl5,1); Float_t lowlim = 0.85; + //Float_t lowlim = 1.; Float_t dedxone = 0.; - Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes; - while ( dedxone < step1 ){ + // Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes; + Float_t step1 = 0.8*f5/atl5; + // while ( dedxone < step1 ){ for (Int_t i = 0; i<22; i++){ for (Int_t j = 1; j>=0; j--){ + if (debug) printf("Acleanplanes: i %i j %i step1 %f dedxone %f epiano[i][j] %f \n",i,j,step1,dedxone,epiano[i][j]); if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j]; + if ( dedxone >= step1 ) break; // new }; + if ( dedxone >= step1 ) break; // new }; - } - if ( dedxone < 0.7 ){ + // } + if ( dedxone < 0.7 ){ // here we could have instead while dedxone == 0. ... perhaps better... for (Int_t i = 0; i<22; i++){ for (Int_t j = 1; j>=0; j--){ + 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]); if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j]; + if ( dedxone >= 0.7 ) break; // new }; + if ( dedxone >= 0.7 ) break; // new }; } // @@ -122,16 +140,20 @@ for (Int_t i = 0; i<22; i++){ for (Int_t j = 1; j>=0; j--){ if ( epiano[i][j] < dedxone*lowlim ){ - // printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius); - epiano[i][j] = 0.; + if ( debug ) printf("Ccleanplanes: %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius); + // epiano[i][j] = 0.; + zpiano[i][j] = true; + if ( epiano[i][j] < dedxone*0.05 ) epiano[i][j] = 0.; } else { - //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); nulliferus = 0; revulsera = true; }; - if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++; + // if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++; + if ( (zpiano[i][j] || epiano[i][j] < 0.7 ) && revulsera ) nulliferus++; if ( nulliferus > 10 ) nullius = true; - if ( nullius ) epiano[i][j] = 0.; + // if ( nullius ) epiano[i][j] = 0.; + if ( nullius ) zpiano[i][j] = true; }; }; @@ -188,7 +210,8 @@ // }; // - this->CleanPlanes(*&epiano); + Bool_t zpiano[22][2]; + this->CleanPlanes(*&epiano, *&zpiano); // PamTrack *ptrack = 0; CaloTrkVar *track = 0; @@ -217,11 +240,14 @@ 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 // + Bool_t zcalo[44]; for(Int_t p=0; p<22; p++){ for(Int_t v=0; v<2; v++){ /*per usare traccia non del calo camboare cibar*/ calorimetro[(2*p)+1-v][0] = L2->GetCaloLevel2()->cibar[p][v];//strip attraversata calorimetro[(2*p)+1-v][1] = epiano[p][v]; //energia del piano //(epiano[p][v])/0.89 + zcalo[(2*p)+1-v] = zpiano[p][v]; + if ( debug ) printf(" idx %i %f %i \n",(2*p)+1-v,epiano[p][v], zpiano[p][v]); }; }; @@ -285,11 +311,20 @@ // Float_t lastok = 0.; + // if ( false ){ // Bool_t goback = false; for ( int o = 0; o < estremi[1][0]; o++ ){ // - if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1]; - if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok; + if (debug) printf(" goforth1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); + if ( calorimetro[o][1] > 0.7 && !zcalo[o] ) lastok = calorimetro[o][1]; + if ( (zcalo[o] || calorimetro[o][1] < 0.7) && lastok > 0. ){ + if ( fabs(calorimetro[o][1]-lastok)/calorimetro[o][1] > 0.5 ) { + if (debug) printf(" goforthchange %f %f \n",calorimetro[o][1],lastok); + calorimetro[o][1] = lastok; + if (debug) printf(" goforthchang+ %f %f \n",calorimetro[o][1],lastok); + } + } + if (debug) printf(" goforth2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); // if ( calorimetro[o][1] < 0.7 ) goback = true; // }; @@ -297,14 +332,20 @@ // if ( goback ){ for ( int o = estremi[1][0]; o >= 0; o-- ){ // - // 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); if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok; - if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1]; - if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok; - // printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); + if ( calorimetro[o][1] > 0.7 && !zcalo[o] ) lastok = calorimetro[o][1]; + if ( (zcalo[o] || calorimetro[o][1] < 0.7) && lastok > 0. ){ + if ( fabs(calorimetro[o][1]-lastok)/calorimetro[o][1] > 0.5 ) { + if ( debug ) printf(" gobackchange \n"); + calorimetro[o][1] = lastok; + } + } + if (debug) printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); // }; // }; + //} if ( startZero ) { estremi[0][0] = 0.; @@ -343,7 +384,9 @@ Float_t dEpianiloop[44]; Int_t tz1=(Int_t)lpz; - Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop + Int_t ti1=(Int_t)lpisotope; + + Enetrack(&tz1, &ti1 , &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop Float_t integ = 0.; @@ -372,17 +415,18 @@ Process(); - this->Draw(0.,0.); + this->Draw(0.,0.,0.); } -void CaloBragg::Draw(Int_t Z, Float_t enetot){ +void CaloBragg::Draw(Int_t Z, Int_t isotope, Float_t enetot){ // Float_t dEpianimean[44]; Float_t dEpianiloop[44]; Float_t Depth[44]; // Int_t tz=(Int_t)qtz; Int_t tz1= Z; + Int_t ti1= isotope; Float_t enet = enetot; // Float_t enet = lpetot; @@ -414,11 +458,12 @@ } else { tz1=(Int_t)lpz; + ti1=(Int_t)lpisotope; enet = lpetot; // Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata } - 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 Float_t sp= spessore[0]*spessore[1]; for(Int_t i=0;i<44;i++)Depth[i]=i*sp; @@ -449,14 +494,14 @@ TString thid2 = Form("hCaloBragg2"); TH2F *th2 = dynamic_cast(gDirectory->FindObject(thid2)); if ( th2 ) th2->Delete(); - 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 th2->SetMarkerStyle(20); th2->SetMarkerColor(kRed); // TString thid3 = Form("hCaloBragg3"); TH2F *th3 = dynamic_cast(gDirectory->FindObject(thid3)); if ( th3 ) th3->Delete(); - 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. th3->SetMarkerStyle(20); th3->SetMarkerColor(kBlue); @@ -491,39 +536,263 @@ void CaloBragg::LoadParam(){ - // - elem[0] = 1.00794; //H 1 - elem[1] = 4.002602; //He 2 - elem[2] = 6.941; //Li 3 - elem[3] = 9.012182;//Be 4 - elem[4] = 10.811; //B 5 - elem[5] = 12.0107; //C 6 - elem[6] = 14.00674;//N 7 - elem[7] = 15.9994; //O 8 - elem[8] = 18.9984032; //F 9 - elem[9] = 20.1797; //Ne 10 - elem[10] = 22.98977;//Na 11 - elem[11] = 24.3050; //Mg 12 - elem[12] = 26.9815; //Al 13 - elem[13] = 28.0855; //Si 14 - elem[14] = 30.974; //P 15 - elem[15] = 32.066; //S 16 - elem[16] = 35.4527; //Cl 17 - elem[17] = 39.948; //Ar 18 - elem[18] = 39.0983; //K 19 - elem[19] = 40.078; //Ca 20 - elem[20] = 44.95591;//Sc 21 - elem[21] = 47.867; //Ti 22 - elem[22] = 50.9415; //V 23 - elem[23] = 51.9961; //Cr 24 - elem[24] = 54.938049;//Mn 25 - elem[25] = 55.845; //Fe 26 - elem[26] = 58.9332; //Co 27 - elem[27] = 58.6934; //Ni 28 - elem[28] = 63.546; //Cu 29 - elem[29] = 65.39; //Zn 30 - elem[30] = 69.723; //Ga 31 - elem[31] = 72.61; //Ge 32 + // elem[Z-1][isotop] 0 is the most common one + // + + elem[0][0] = 1.00782; //H 1 + elem[0][1] = 2.01410; // 2H (Isotope) + elem[0][2] = -1.; + elem[0][3] = -1.; + elem[0][4] = -1.; + elem[0][5] = -1.; + elem[0][6] = -1.; + + elem[1][0] = 4.002602; //He 2 + elem[1][1] = 3.016029; // 3He (Isotope) + elem[1][2] = -1.; + elem[1][3] = -1.; + elem[1][4] = -1.; + elem[1][5] = -1.; + elem[1][6] = -1.; + + elem[2][0] = 7.016004; //Li 3 + elem[2][1] = 6.015123; //6Li (Isotope) + elem[2][2] = -1.; + elem[2][3] = -1.; + elem[2][4] = -1.; + elem[2][5] = -1.; + elem[2][6] = -1.; + + elem[3][0] = 9.012182; //Be 4 + elem[3][1] = 10.01353; //10Be (Isotope) (most stable) + elem[3][2] = 7.01693; //9Be no EC in space? + elem[3][3] = -1.; + elem[3][4] = -1.; + elem[3][5] = -1.; + elem[3][6] = -1.; + + elem[4][0] = 11.00930; //B 5 + elem[4][1] = 10.01294; //10B (Isotope) + elem[4][2] = -1.; + elem[4][3] = -1.; + elem[4][4] = -1.; + elem[4][5] = -1.; + elem[4][5] = -1.; + + elem[5][0] = 12.0107; //C 6 + elem[5][1] = 13.00335; //13C (Isotope) + elem[5][2] = -1.; + elem[5][3] = -1.; + elem[5][4] = -1.; + elem[5][5] = -1.; + elem[5][5] = -1.; + + elem[6][0] = 14.00674; //N 7 + elem[6][1] = 15.00011; //15N (Isotope) + elem[6][2] = -1.; + elem[6][3] = -1.; + elem[6][4] = -1.; + elem[6][5] = -1.; + elem[6][5] = -1.; + + elem[7][0] = 15.99491; //O 8 + elem[7][1] = 17.99916; //18O (Isotope) + elem[7][2] = 16.99916; //17O (Isotope) + elem[7][3] = -1.; + elem[7][4] = -1.; + elem[7][5] = -1.; + elem[7][5] = -1.; + + elem[8][0] = 18.99840; //F 9 + elem[8][1] = -1.; + elem[8][2] = -1.; + elem[8][3] = -1.; + elem[8][4] = -1.; + elem[8][5] = -1.; + elem[8][5] = -1.; + + elem[9][0] = 19.99244; //Ne 10 + elem[9][1] = 21.99138; //22Ne (Isotope) + elem[9][2] = 20.99384; //21Ne 10 + elem[9][3] = -1.; + elem[9][4] = -1.; + elem[9][5] = -1.; + elem[9][6] = -1.; + + elem[10][0] = 22.98977; //Na 11 + elem[10][1] = 21.99444; //22Na (Isotope) (most stable) + elem[10][2] = -1.; + elem[10][3] = -1.; + elem[10][4] = -1.; + elem[10][5] = -1.; + elem[10][6] = -1.; + + elem[11][0] = 23.98504; //Mg 12 + elem[11][1] = 25.98259; //26Mg (Isotope) + elem[11][2] = 24.98504; //25Mg (Isotope) + elem[11][3] = -1.; + elem[11][4] = -1.; + elem[11][5] = -1.; + elem[11][6] = -1.; + + elem[12][0] = 26.98154; //Al 13 + elem[12][1] = 25.98489; //26Al (Isotope) (most stable) + elem[12][2] = -1.; + elem[12][3] = -1.; + elem[12][4] = -1.; + elem[12][5] = -1.; + elem[12][6] = -1.; + + elem[13][0] = 27.97692; //Si 14 + elem[13][1] = 28.97649; //29Si (Isotope) + elem[13][2] = 29.97377; //30Si (Isotope) + elem[13][3] = -1.; + elem[13][4] = -1.; + elem[13][5] = -1.; + elem[13][6] = -1.; + + elem[14][0] = 30.97376; //P 15 + elem[14][1] = -1.; + elem[14][2] = -1.; + elem[14][3] = -1.; + elem[14][4] = -1.; + elem[14][5] = -1.; + elem[14][6] = -1.; + + elem[15][0] = 31.97207; //S 16 + elem[15][1] = 33.96787; //34S (Isotope) + elem[15][2] = 32.97146; //33S (Isotope) + elem[15][3] = 35.96708; //36S (Isotope) + elem[15][4] = -1.; + elem[15][5] = -1.; + elem[15][6] = -1.; + + elem[16][0] = 34.96885; //Cl 17 + elem[16][1] = 36.96831; //37Cl 17 + elem[16][2] = 35.96890; //36Cl (Isotope) + elem[16][3] = -1.; + elem[16][4] = -1.; + elem[16][5] = -1.; + elem[16][6] = -1.; + + elem[17][0] = 39.962383; //Ar 18 + elem[17][1] = 35.967545; //36Ar (Isotope) + elem[17][2] = 37.962732; //38Ar (Isotope) + elem[17][3] = 38.964313; //39Ar (Isotope) + elem[17][4] = -1.; + elem[17][5] = -1.; + elem[17][6] = -1.; + + elem[18][0] = 38.963707; //K 19 + elem[18][1] = 40.961825; //41K (Isotope) + elem[18][2] = 39.963998; //40K (Isotope) + elem[18][3] = -1.; + elem[18][4] = -1.; + elem[18][5] = -1.; + elem[18][6] = -1.; + + elem[19][0] = 39.962590; //Ca 20 + elem[19][1] = 43.955482; //44Ca (Isotope) + elem[19][2] = 41.958618; //42Ca (Isotope) + elem[19][3] = 42.958767; //43Ca (Isotope) + elem[19][4] = 45.953693; //46Ca (Isotope) + elem[19][5] = 40.962278; //41Ca (Isotope) + elem[19][6] = -1.; + + elem[20][0] = 44.955912;//Sc 21 + elem[20][1] = -1.; + elem[20][2] = -1.; + elem[20][3] = -1.; + elem[20][4] = -1.; + elem[20][5] = -1.; + elem[20][6] = -1.; + + elem[21][0] = 47.947946; //Ti 22 + elem[21][1] = 45.952632; //46Ti (Isotope) + elem[21][2] = 46.951763; //47Ti (Isotope) + elem[21][3] = 48.947870; //49Ti (Isotope) + elem[21][4] = 49.944791; //50Ti (Isotope) + elem[21][5] = 43.959690; //44Ti (Isotope) (half life 60y) + elem[21][6] = -1.; + + elem[22][0] = 50.943960; //V 23 + elem[22][1] = 49.947158; //50V (Isotope) + elem[22][2] = -1.; + elem[22][3] = -1.; + elem[22][4] = -1.; + elem[22][5] = -1.; + elem[22][6] = -1.; + + elem[23][0] = 51.940507; //Cr 24 + elem[23][1] = 52.940649; //53Cr (Isotope) + elem[23][2] = 49.946044; //50Cr (Isotope) + elem[23][3] = 53.938880; //54Cr (Isotope) + elem[23][4] = -1.; + elem[23][5] = -1.; + elem[23][6] = -1.; + + elem[24][0] = 54.938049;//Mn 25 + elem[24][1] = 52.941290;//53Mn (Isotope) + elem[24][2] = -1.; + elem[24][3] = -1.; + elem[24][4] = -1.; + elem[24][5] = -1.; + elem[24][6] = -1.; + + elem[25][0] = 55.934937; //Fe 26 + elem[25][1] = 53.939610; //54Fe (Isotope) + elem[25][2] = 56.935394; //57Fe (Isotope) + elem[25][3] = 57.933276; //58Fe (Isotope) + elem[25][4] = 59.934072; //58Fe (Isotope) + + elem[26][0] = 58.933195; //Co 27 + elem[26][1] = 59.933817; //60Co (Isotope) + elem[26][2] = -1.; + elem[26][3] = -1.; + elem[26][4] = -1.; + elem[26][5] = -1.; + elem[26][6] = -1.; + + + elem[27][0] = 57.935343; //Ni 28 + elem[27][1] = 61.928345; //62Ni (Isotope) + elem[27][2] = 59.930786; //60Ni (Isotope) + elem[27][3] = 60.931056; //61Ni (Isotope) + elem[27][4] = 63.927966; //64Ni (Isotope) + elem[27][5] = 58.934346; //59Ni (Isotope) + elem[27][6] = -1.; + + elem[28][0] = 62.929597; //Cu 29 + elem[28][1] = 64.927789; //65Cu (Isotope) + elem[28][2] = -1.; + elem[28][3] = -1.; + elem[28][4] = -1.; + elem[28][5] = -1.; + elem[28][6] = -1.; + + elem[29][0] = 63.929142; //Zn 30 + elem[29][1] = 65.926033; //66Zn (Isotope) + elem[29][2] = 67.924844; //68Zn (Isotope) + elem[29][3] = 66.927127; //67Zn (Isotope) + elem[29][4] = 69.925319; //70Zn (Isotope) + elem[29][5] = -1.; + elem[29][6] = -1.; + + elem[30][0] = 68.925573; //Ga 31 + elem[30][1] = 70.924701; //71Ga (Isotope) + elem[30][2] = -1.; + elem[30][3] = -1.; + elem[30][4] = -1.; + elem[30][5] = -1.; + elem[30][6] = -1.; + + elem[31][0] = 73.921177; //Ge 32 + elem[31][1] = 71.922075; //72Ge (Isotope) + elem[31][2] = 69.924247; //70Ge (Isotope) + elem[31][3] = 75.921403; //76Ge (Isotope) + elem[31][4] = 73.923459; //73Ge (Isotope) + elem[31][5] = -1.; + elem[31][6] = -1.; //parametri calorimetro @@ -650,7 +919,7 @@ -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){ /*perdita di energia per ioni pesanti (come da routine geant)*/ // in : dx => spessore g/cm2 @@ -667,7 +936,7 @@ Float_t dEP=0.; // gamma // Mass = A * MassP; /*in Mev/c2*/ - gam = (*Etot)/(elem[*Z-1]*MassP); // E = gamma M c2 + gam = (*Etot)/(elem[*Z-1][*isotope]*MassP); // E = gamma M c2 Bet = sqrt((SQ(gam) -1.)/SQ(gam)); @@ -693,7 +962,7 @@ -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[]){ //calcola energia rilasciata sulla traccia (usa ELOSS) // in : Z =>carica @@ -711,20 +980,23 @@ //azzero energia rilasciata sui piani memset(out, 0, 2*NPLA*sizeof(Float_t)); - Float_t Massa = (elem[(*Z)-1] * MassP); + Float_t Massa = (elem[(*Z)-1][*isotope] * MassP); for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){ dE=0.; //spessore silicio corretto x inclinazione, z, energia, out:rilascio - ELOSS(&spessore[0], Z, &Ezero, &dE, ISi);//spessore in g/cm2!! - if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop + ELOSS(&spessore[0], Z , isotope , &Ezero, &dE, ISi);//spessore in g/cm2!! + + if(dE!=dE) return; //controlla che non sia un NaN + + if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop out[ipla] = Ezero - Massa; //MeV return; }else{ out[ipla] = dE; //MeV Ezero = Ezero - dE;//energia residua - 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); }; //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno if(ipla%2 == 0){ @@ -739,7 +1011,7 @@ sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato in g/cm2 } // printf(" sp %f II %f \n",sp,II); - ELOSS(&sp, Z, &Ezero, &dE,II); + ELOSS(&sp, Z, isotope , &Ezero, &dE,II); if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop return; }else{ @@ -814,7 +1086,7 @@ //meno peso ai piani con rilasci maggiori di 1000 MIP // if(calorimetro[ipla][1] > 1000) wi=0.5; if(calorimetro[ipla][1] > 1200.) wi=0.5; - + if(debug) printf("chiquadro start \n "); Float_t arg = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP)); sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP)); @@ -855,10 +1127,12 @@ //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2 // - + //printf("entrato"); memset(dEplan,0,2*NPLA*sizeof(Float_t)); Int_t Z = 0;// z iniziale + + Int_t isotope=0; Float_t Massa = 0.; @@ -872,6 +1146,7 @@ Int_t max=32;//max z di cui so la massa :P if((Zlimite)<=31) max=(int)(Zlimite) + 1; +if(debug) printf("loopze inizio max %d \n",max); if ( fzeta > 0. ){ zmin = fzeta; max = fzeta+1; @@ -879,9 +1154,10 @@ Int_t colmax=32; Int_t rowmax=3000; + Int_t isomax=7; - Float_t matrixchi2[colmax][rowmax][3]; - memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t)); + Float_t matrixchi2[colmax][isomax][rowmax][3]; + memset(matrixchi2, 0, colmax*isomax*rowmax*3*sizeof(Float_t)); Int_t imin = 1-nostep/2; Int_t imax = nostep/2; @@ -890,51 +1166,88 @@ for(Int_t inucl=zmin; inucl 0. ){ + + if( fiso != -1 ){ + isotope=fiso; + if(debug) printf("In Loopze - Isotope N %d",isotope); + } + Massa = elem[inucl-1][isotope]*MassP; + //loop energia - Int_t iene2 = 0; + Int_t iene2 = -1; + // 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= imin; iene 0. ) energia=fene; //forza l'energia + if (debug) printf("loopze energia %f, z %d, isotopo %d ,iene %d\n",energia,Z,isotope,iene); +// printf(" energia %f , forzata %f \n",energia,fene); + Enetrack(&Z, &isotope, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani - iene2++; - energia= Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica.. - - Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani - //calcolo chi2 - chiquadro(dEplan,chi2); - + chiquadro(dEplan,chi2); //calcolo chi2 + if (debug) printf("loopze chi %f \n",chi2[0]); + if(debug && TMath::Finite(chi2[0])==1 && (TMath::IsNaN(chi2[0])!=1) ) printf("loopze fin mat %f \n",chi2[0]); // printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]); - - if( (chi2[1] != 79.) ){//salto quelli che frammentano - matrixchi2[inucl][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia - matrixchi2[inucl][iene2][1]=energia;//energia per questo chi2 - matrixchi2[inucl][iene2][2]=chi2[2];//piani saltati nel chi2 - } else { - matrixchi2[inucl][iene2][0]=1000.;//valore chi2 per questo z a questa energia - matrixchi2[inucl][iene2][1]=1000.;//energia per questo chi2 - matrixchi2[inucl][iene2][2]=1000.;//piani saltati nel chi2 - break; + if( (chi2[1] != 79.) ){//salto quelli che frammentano + matrixchi2[inucl][isotope][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia + matrixchi2[inucl][isotope][iene2][1]=energia;//energia per questo chi2 + matrixchi2[inucl][isotope][iene2][2]=chi2[2];//piani saltati nel chi2 + if( fene > 0. ) break; + } else { + matrixchi2[inucl][isotope][iene2][0]=numeric_limits::max();//valore chi2 per questo z a questa energia + matrixchi2[inucl][isotope][iene2][1]=numeric_limits::max();//energia per questo chi2 + matrixchi2[inucl][isotope][iene2][2]=numeric_limits::max();//piani saltati nel chi2 + break; } - }//fine loop energia + }//fine loop energia + + if( fiso != -1 ){ + if(debug) printf("exited form isotopes loop"); + break; + } + + isotope++; //incremento il contatore isotopi + }//fine loop isotopi + isotope=0; //riazzero il contatore isotopi - };//fine loop z + }//fine loop z + isotope=0;//non dovrebbe servire //Emi - for (Int_t nu=zmin; nu0.)){ - bestchi2[0]= matrixchi2[nu][en][0];// chi2 - bestchi2[1]= (Float_t)nu; // z - bestchi2[2]= matrixchi2[nu][en][1];//energia; - bestchi2[3]= matrixchi2[nu][en][2];// totale piani saltati - } - } + for (Int_t nu=zmin; nu 0.){ + for (Int_t en=0; en0.)){ + bestchi2[0]= matrixchi2[nu][isotope][en][0];// chi2 + bestchi2[1]= (Float_t)nu; // z + bestchi2[2]= matrixchi2[nu][isotope][en][1];//energia; + bestchi2[3]= matrixchi2[nu][isotope][en][2];// totale piani saltati + bestchi2[4]= (Float_t)isotope; //isotopo + } + } + + if( fiso != -1 ){ + if(debug) printf("exited form isotopes loop"); + break; } + isotope++; + } + isotope=0; + } };//endloopze @@ -1058,12 +1371,13 @@ Float_t zmin=1.; Float_t zmax=32.; - Float_t bestchitemp[4] = {0,0,0,0}; + Float_t bestchitemp[5] = {0,0,0,0,0}; - bestchi2[0]=10000.; + bestchi2[0]=numeric_limits::max(); bestchi2[1]=0.; bestchi2[2]=0.; bestchi2[3]=0.; + bestchi2[4]=0.; Float_t zero=0.; //------------primo loop ---------------------- // energia ezero, zstart zstop @@ -1073,34 +1387,38 @@ loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200); // loopze(Integrale*2.,Integrale/100.,zmin,zmax); - if ( debug ) printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]); + if ( debug) printf("Zdaloop start Integrale %f , outene %f \n",Integrale,bestchi2[2]); //------------secondo loop ---------------------- - for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i]; - bestchi2[0] = 10000.; + for(Int_t i=0;i<5;i++) bestchitemp[i]=bestchi2[i]; + bestchi2[0]=numeric_limits::max(); bestchi2[1] = 0.; bestchi2[2] = 0.; - bestchi2[3] = 0.;//riazzero + bestchi2[3] = 0.; + bestchi2[4] = 0.;//riazzero Float_t step = bestchitemp[2];// zero=0.; // qualsiasi altro valore peggiora le cose // zmin=zmax=bestchitemp[1]; zmin=bestchitemp[1]-1; zmax=bestchitemp[1]+1; - // loopze(step,zero,zmin,zmax); // + //loopze(step,zero,zmin,zmax); // //-> loopze(step,step/2.,zmin,zmax,200); // loopze(step,step/2.,zmin,zmax,500); // - if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step); + //step = bestchitemp[2];// + + //loopze(step/2,3*step/4.,zmin,zmax,500); // + + if ( debug ) printf("Zdaloop Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step); - //chi2,z,Etot,Pskip lpchi2=bestchi2[0]; lpz=bestchi2[1]; lpetot=bestchi2[2]; lppskip=bestchi2[3]; - + lpisotope=bestchi2[4]; };//endZdaloop