--- calo/flight/CaloBragg/src/CaloBragg.cpp 2011/01/28 10:40:44 1.8 +++ calo/flight/CaloBragg/src/CaloBragg.cpp 2011/06/17 11:56:42 1.12 @@ -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.; @@ -42,7 +46,7 @@ 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.; @@ -69,6 +73,7 @@ printf(" Z from loop %f: \n", lpz); printf(" energy from loop %f: \n", lpetot); printf(" plane not used for loop %f: \n", lppskip); + printf(" ndf: %i \n",ndf); printf("========================================================================\n"); // }; @@ -301,7 +306,10 @@ }; // }; - + 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 +328,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 +338,87 @@ }; +Float_t CaloBragg::Integral(){ + Process(); + + 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 + + + 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.); + +} + +void CaloBragg::Draw(Int_t Z, 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; + 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; + 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 Float_t sp= spessore[0]*spessore[1]; for(Int_t i=0;i<44;i++)Depth[i]=i*sp; @@ -392,7 +471,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"); @@ -411,14 +493,14 @@ // elem[0] = 1.00794; //H 1 - elem[1] = 4.0026; //He 2 + 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.9984; //F 9 + elem[8] = 18.9984032; //F 9 elem[9] = 20.1797; //Ne 10 elem[10] = 22.98977;//Na 11 elem[11] = 24.3050; //Mg 12 @@ -465,7 +547,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 +585,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 +595,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 +625,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 +650,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, Float_t *Etot, Float_t *out, Float_t II){ /*perdita di energia per ioni pesanti (come da routine geant)*/ // in : dx => spessore g/cm2 @@ -586,15 +672,18 @@ 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; @@ -627,7 +716,7 @@ 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!! + 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 out[ipla] = Ezero - Massa; //MeV return; @@ -635,13 +724,22 @@ }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, &Ezero, &dE,II); if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop return; }else{ @@ -671,62 +769,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; + //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)); + 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 + 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 @@ -743,7 +856,6 @@ // - Float_t dEplan[2*NPLA];//energia rilasciata calcolata memset(dEplan,0,2*NPLA*sizeof(Float_t)); Int_t Z = 0;// z iniziale @@ -755,9 +867,15 @@ 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 ( fzeta > 0. ){ + zmin = fzeta; + max = fzeta+1; + } Int_t colmax=32; Int_t rowmax=3000; @@ -769,7 +887,7 @@ Int_t imax = nostep/2; //loop elementi - for(Int_t inucl=(int)(Zstart); inucl0.)){ bestchi2[0]= matrixchi2[nu][en][0];// chi2 @@ -949,9 +1068,12 @@ //------------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(" Integrale %f , outene %f \n",Integrale,bestchi2[2]); //------------secondo loop ---------------------- for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i]; @@ -966,7 +1088,11 @@ zmin=bestchitemp[1]-1; zmax=bestchitemp[1]+1; // loopze(step,zero,zmin,zmax); // - loopze(step,step/2.,zmin,zmax,200); // + + //-> 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); //chi2,z,Etot,Pskip