--- calo/flight/CaloBragg/src/CaloBragg.cpp 2011/01/28 10:40:44 1.8 +++ calo/flight/CaloBragg/src/CaloBragg.cpp 2014/01/23 11:23:43 1.16 @@ -26,11 +26,15 @@ debug = false; usetrack = false; usepl18x = false; + newchi2 = false; + usenewBB = false; + fzeta = -1.; // }; void CaloBragg::Clear(){ // + ndf = 0; tr = 0; sntr = 0; // qtchi2 = 0.; @@ -39,10 +43,12 @@ // qtpskip = 0.; lpchi2 = 0.; lpz = 0.; + lpisotope= 0.; lpetot = 0.; lppskip = 0.; + memset(calorimetro,0,44*2*sizeof(Float_t)); - memset(spessore,0,3*sizeof(Float_t)); + memset(spessore,0,4*sizeof(Float_t)); memset(estremi,0,2*2*sizeof(Float_t)); Integrale=0.; @@ -67,8 +73,10 @@ // 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); printf("========================================================================\n"); // }; @@ -84,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 }; } // @@ -117,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; }; }; @@ -183,7 +210,8 @@ // }; // - this->CleanPlanes(*&epiano); + Bool_t zpiano[22][2]; + this->CleanPlanes(*&epiano, *&zpiano); // PamTrack *ptrack = 0; CaloTrkVar *track = 0; @@ -212,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]); }; }; @@ -280,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; // }; @@ -292,16 +332,25 @@ // 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.; + // estremi[0][1] = 0.; + } /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/ for(Int_t pl=0; pl<(2*NPLA); pl++){ @@ -320,6 +369,7 @@ // mediatroncata(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip /*z ed energia con loop*/ + if ( debug ) printf(" call Zdaloop with integrale %f \n",Integrale); Zdaloop(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip @@ -329,17 +379,91 @@ }; +Float_t CaloBragg::Integral(){ + Process(); + + Float_t dEpianiloop[44]; + Int_t tz1=(Int_t)lpz; + 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.; + for(Int_t i=0;i<=estremi[1][0];i++){ + // integ += dEplan[i]; + //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]); + integ += dEpianiloop[i]; + // printf(" step %i integ %f deplan %f \n",i,integ,dEpianiloop[i]); + } + return integ; +} + +Float_t CaloBragg::LastIntegral(){ + Process(); + + Float_t integ = 0.; + for(Int_t i=0;i<=estremi[1][0];i++){ + integ += dEplan[i]; + //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]); + } + return integ; +} + + void CaloBragg::Draw(){ Process(); + this->Draw(0.,0.,0.); + +} + +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=(Int_t)lpz; - // Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata - Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop + Int_t tz1= Z; + Int_t ti1= isotope; + Float_t enet = enetot; + // Float_t enet = lpetot; + + if ( Z > 0. && enetot > 0. ){ + estremi[0][0] = 0; + estremi[1][0] = 43; + + + Float_t ytgx = 0.; + Float_t ytgy = 0.; + + //lunghezza effettiva di silicio attraversata (mm) + Float_t SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy)); + + spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2 + + /*tungsteno*/ + + //rapporto tra rilasci energetici nei due materiali + Float_t WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad + //gcm2W = WCross/10. * rhoW; + + // (g/cm2W)/(g/cm2Si) + spessore[3] = (WCross/10.) * rhoW; + Float_t a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664); //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664); + spessore[1] = a; + //riscala mip allo spessore attraversato + spessore[2] = MIP*(SiCross/ySi); + + } 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, &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; @@ -370,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); @@ -392,7 +516,10 @@ // tc->cd(2); tc->cd(); // - 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++){ + th3->Fill(Depth[i],dEpianiloop[i]); + // printf(" i %i Depth %f depianiloop %f \n",i,Depth[i],dEpianiloop[i]); + } th3->Draw(); th2->Draw("same"); @@ -409,39 +536,263 @@ void CaloBragg::LoadParam(){ - // - elem[0] = 1.00794; //H 1 - elem[1] = 4.0026; //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.9984; //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 @@ -465,7 +816,10 @@ pigr = 3.1415; Na = 6.02e-23; ZA = 0.49; /*Z/A per Si*/ - ISi =182e-06; /*MeV*/ + // ISi =182e-06; /*MeV*/ + ISi = 171e-06; /*MeV*/ + IW = 735e-06; /*MeV*/ + // ISi =0.0001059994; /*GeV!!*/ no era giusto!! Me = 0.511; /* MeV*/ MassP = 931.27;/*MeV*/ r2 = 7.95e-26; /*ro*ro in cm */ @@ -500,7 +854,7 @@ //lunghezza effettiva di silicio attraversata (mm) SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy)); - spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2 + spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2 /*tungsteno*/ ytgx = yW * L2->GetCaloLevel2()->tanx[0]; @@ -510,21 +864,19 @@ WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad //gcm2W = WCross/10. * rhoW; - a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664); //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664); - // (g/cm2W)/(g/cm2Si) + spessore[3] = (WCross/10.) * rhoW; + a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664); //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664); spessore[1] = a; - //riscala mip allo spessore attraversato spessore[2] = MIP*(SiCross/ySi); - };//end conversione -void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out){ +void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out, Float_t II){ //rilascio energetico con bethe bloch con correzioni //in: x: g/cm2 @@ -542,16 +894,19 @@ Float_t lg =0.; Float_t Energia=0.; Float_t C=0.; + Float_t INo = ISi; + + if ( usenewBB ) INo = II; eta = (*gam)*(*Bet); //Bet=3/gam; SQ(*gam) * SQ(*Bet) Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass)); - lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi); + lg = 2.* Me * SQ(eta) * Wmax / SQ(INo); // Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA* SQ(z)/SQ(Bet) * lg; - C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(ISi,2.) + - (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(ISi,3.); + C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(INo,2.) + + (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(INo,3.); if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653); @@ -564,7 +919,7 @@ -void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out){ +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 @@ -581,20 +936,23 @@ 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)); - v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.)))); + // v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.)))); + v= 121.4139*(Bet*pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet*pow((*Z),(2./3.)))); // EMI AAAAGGH!! //carica effettiva Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v)); //perdita energia per un protone Float_t protone =1.; - Float_t Mass=(elem[*Z-1]*MassP); - BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); + // Float_t Mass=(elem[*Z-1]*MassP); //EMI + // BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); + + BetheBloch(dx, &protone, &MassP, &gam, &Bet, &dEP, II);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); //EMI *out= (SQ(Q)*(dEP));//*dx; @@ -604,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 @@ -622,26 +980,38 @@ //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);//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); }; //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno if(ipla%2 == 0){ /*tungsteno*/ dE=0.; - Float_t sp= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato in g/cm2 - ELOSS(&sp, Z, &Ezero, &dE); + Float_t sp = 0.; + Float_t II = ISi; + if ( usenewBB ){ + sp = spessore[3]; + II = IW; + } else { + sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato in g/cm2 + } + // printf(" sp %f II %f \n",sp,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{ @@ -671,62 +1041,77 @@ Float_t badplane=0.; Float_t badplanetot=0.; Float_t w,wi; - - for(Int_t ipla=0; ipla<2*NPLA; ipla++){ - //tutti i piani attraversati dalla traiettoria - if(calorimetro[ipla][0] != -1.){ // - w=0.; //normalizzazione; - wi=1.;//peso + // + if ( newchi2 ){ + ndf = 0; + sum = 0.; + for( Int_t ipla=((int)(estremi[0][0])); ipla<= ((int)(estremi[1][0])); ipla++){ + sum += pow((dE[ipla] - (calorimetro[ipla][1] * spessore[2]))/(0.05*dE[ipla]),2.); + // printf(" quiqui: dE %f calor %f spessore[2] %f \n",dE[ipla],spessore[2]*calorimetro[ipla][1],spessore[2]); + ndf++; + } + ndf -= 2; + if ( ndf > 0 ) sum /= (float)ndf; + out[0] = sum; + out[1] = 0.; + out[2] = (int)(estremi[1][0])-ndf; + // printf(" sum %f ndf %i \n ",sum,ndf); + } else { + for(Int_t ipla=0; ipla<2*NPLA; ipla++){ + //tutti i piani attraversati dalla traiettoria + if(calorimetro[ipla][0] != -1.){ // + w=0.; //normalizzazione; + wi=1.;//peso - //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso - if (iplaestremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.; - if((ipla>estremi[1][0])) wi=0.; + //tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0 + //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.; + if((ipla>estremi[1][0])) wi=0.; - //normalizzazione - if (calorimetro[ipla][1] != 0.) w=1./(calorimetro[ipla][1]* MIP); // + //normalizzazione + if (calorimetro[ipla][1] != 0.) w=1./(calorimetro[ipla][1]* MIP); // - //tolgo piani con rilasci inferiori al 30% del precedente - if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2 - wi=0.; - //se sono piani intermedi (non si e' fermta) li considero non buoni - if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){// - badplane+=1.; - badplanetot+=1.; - }; - }; - - //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; - - 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)); - if(debug){ - printf("dedx calcolata %f e reale %f \n",dE[ipla],(calorimetro[ipla][1] * MIP)); - } - //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente - if(wi != 0.){// - PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente - badplane = 0.;//azzero contatore piani scartati consecutivi + //tolgo piani con rilasci inferiori al 30% del precedente + if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2 + wi=0.; + //se sono piani intermedi (non si e' fermta) li considero non buoni + if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){// + badplane+=1.; + badplanetot+=1.; + }; + }; + + //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)); + if(debug){ + printf("dedx calcolata %f e reale %f \n",dE[ipla],(calorimetro[ipla][1] * MIP)); + } + //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente + if(wi != 0.){// + PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente + badplane = 0.;//azzero contatore piani scartati consecutivi + }; }; - }; - //da Emi - if(badplane > 2){ - // printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n"); - out[1] =79.; - break; - }; - - };//fine loop piani - //chi2,frammentato,pskip - out[0]=sum; - out[2]=badplanetot; + //da Emi + if(badplane > 2){ + // printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n"); + out[1] =79.; + break; + }; + };//fine loop piani + //chi2,frammentato,pskip + out[0]=sum; + out[2]=badplanetot; + } };//end chiquadro @@ -742,11 +1127,12 @@ //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2 // - - Float_t dEplan[2*NPLA];//energia rilasciata calcolata + //printf("entrato"); memset(dEplan,0,2*NPLA*sizeof(Float_t)); Int_t Z = 0;// z iniziale + + Int_t isotope=0; Float_t Massa = 0.; @@ -755,67 +1141,113 @@ Float_t energia =0.;//energia del loop Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati - + + Int_t zmin = (int)Zstart; 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; + } 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; //loop elementi - for(Int_t inucl=(int)(Zstart); 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 - 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; + 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][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=(int)(Zstart); 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 @@ -939,42 +1371,54 @@ 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.; - Float_t zero=0.; + bestchi2[4]=0.; + // Float_t zero=0.; //------------primo loop ---------------------- // energia ezero, zstart zstop // loopze(Integrale,zero,zmin,zmax); - loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50); + + //-> loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50); + loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200); + // loopze(Integrale*2.,Integrale/100.,zmin,zmax); - // 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 + // 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,step/2.,zmin,zmax,200); // + //loopze(step,zero,zmin,zmax); // + + //-> loopze(step,step/2.,zmin,zmax,200); // + loopze(step,step/2.,zmin,zmax,500); // + + //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