--- calo/flight/CaloBragg/src/CaloBragg.cpp 2008/06/13 08:11:04 1.1.1.1 +++ calo/flight/CaloBragg/src/CaloBragg.cpp 2008/06/18 15:24:12 1.2 @@ -109,7 +109,6 @@ // // - // Always calculate stdedx1 // Int_t view = 0; Int_t plane = 0; @@ -134,7 +133,7 @@ ptrack = L2->GetTrack(ntr); if ( ptrack ) track = ptrack->GetCaloTrack(); } else { - track = L2->GetCaloStoredTrack(ntr); + track = L2->GetCaloStoredTrack(ntr); //al momento e' vera solo questa riga }; // if ( !track && ntr >= 0 ){ @@ -155,26 +154,23 @@ // for(Int_t p=0; p<22; p++){ for(Int_t v=0; v<2; v++){ - cout<< L2->GetCaloLevel2()->cibar[p][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 }; }; - - /*per ogni evento clacolo la conversione mip e w attraversato in equivalente Si*/ - + /*per ogni evento calcolo la conversione mip e w attraversato in equivalente Si*/ conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione + //cout<<"spessore= "<50.)){ if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ //0.7 soglia minima estremi[0][0]=p; @@ -185,51 +181,37 @@ //ultimo parte da 44 e sale p=43; while( (estremi[1][1] == 0) && (p>(int)estremi[0][0]) ){ - cout<<"piano "<

50.)){ - if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ - estremi[1][0]=p; + if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ + estremi[1][0]=p;//era p estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV }; p = p-1; }; // - cout<<"estremi: in "<Print(); - if ( debug ) printf(" esci \n"); + if ( debug ) printf(" fine evento \n"); // }; @@ -256,15 +238,14 @@ Float_t Depth[44]; Int_t tz=(Int_t)qtz; Int_t tz1=(Int_t)lpz; - Enetrack(&tz, &qtetot, &estremi[0][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata - Enetrack(&tz1, &lpetot, &estremi[0][0], dEpianiloop);//calcola rilascio energetico sui piani da loop + 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 Float_t sp= spessore[0]*spessore[1]; for(Int_t i=0;i<44;i++)Depth[i]=i*sp; // gStyle->SetLabelSize(0.04); gStyle->SetNdivisions(510,"XY"); - //gStyle->SetLogy(); // TString hid = Form("cCaloBragg"); @@ -304,14 +285,14 @@ tc->cd(1); // - for(Int_t i=0;i<44;i++)th->Fill(Depth[i],dEpianimean[i]); - for(Int_t i=0;i<44;i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP); + for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]); + for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP); th->Draw(); th2->Draw("same"); tc->cd(2); // - for(Int_t i=0;i<44;i++)th3->Fill(Depth[i],dEpianiloop[i]); + for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]); th3->Draw(); th2->Draw("same"); @@ -325,13 +306,6 @@ }; -void CaloBragg::SWAP( Float_t *A, Float_t *B ){ - Float_t Tmp = *A; - *A = *B; - *B = Tmp; -}; - - void CaloBragg::LoadParam(){ @@ -370,7 +344,6 @@ elem[31] = 72.61; //Ge 32 - //parametri calorimetro NPLA = 22; NCHA = 96; @@ -472,7 +445,6 @@ Float_t lg =0.; Float_t Energia=0.; Float_t C=0.; - eta = (*gam)*(*Bet); @@ -481,8 +453,8 @@ lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi); // 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(ISi,2.) + + (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(ISi,3.); if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653); // spessorecm x ??/massSi x Zsi @@ -534,14 +506,13 @@ *out= (SQ(Q)*(dEP));//*dx; -};//end ELOSS - +};//end ELOSS -void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo, Float_t out[]){ +void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo,Float_t* ultimo, Float_t out[]){ //calcola energia rilasciata sulla traccia (usa ELOSS) // in : Z =>carica @@ -560,16 +531,15 @@ memset(out, 0, 2*NPLA*sizeof(Float_t)); Float_t Massa = (elem[(*Z)-1] * MassP); - //loop piani (dal primo in cui entra) - for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){ + // for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){ + 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 out[ipla] = Ezero - Massa; //MeV return; @@ -585,7 +555,7 @@ dE=0.; Float_t sp= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato in g/cm2 ELOSS(&sp, Z, &Ezero, &dE); - + // cout<<"perdita per piano di W ="<estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=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); //tolgo piani con rilasci inferiori al 30% del precedente if(calorimetro[ipla][1] <= (0.7*PianoPrecedente)){ - wi=0.; //se sono piani intermedi (non si è fermta) li considero non buoni - if( (ipla < estremi[1][0]) && (calorimetro[ipla][1] !=0.)){ + 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; - //do peso maggiore alle ultime 6 misure - // if((ipla >= estr2[1][0]-12)) wi=0.; - 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)); @@ -669,8 +635,9 @@ };//fine loop piani //chi2,frammentato,pskip out[0]=sum; - out[2]=badplanetot; + + if(out[1] ==79.)cout<<"frammentato !!!!!"< energia iniziale (intergale) // Zstart => minimo z da cui patire // Zlimite => z a cui fermarsi (z al minimo di ionizz sul 1o piano) - // spessore => array conversione spessore Si, mip, W - // estr1 => array con primo[0][0] e ultimo[1][0] piano attraversati ed energie([x][1]) - // calo2[44][2]=> [0]strip attraversata [1]energia misurata su ognuno dei 44 piani // //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2 // + + Float_t dEplan[2*NPLA];//energia rilasciata calcolata memset(dEplan,0,2*NPLA*sizeof(Float_t)); @@ -698,7 +664,7 @@ Float_t Massa = 0.; - Float_t Stepint =(step)/1000.;//passo per il calcolo di energia + Float_t Stepint =(step)/100.;//passo per il calcolo di energia //era 1000 Float_t energia =0.;//energia del loop @@ -722,11 +688,11 @@ Massa = elem[inucl-1]*MassP; //loop energia - for(Int_t iene= 0; iene<3000; iene++){ + for(Int_t iene= 0; iene<1000; iene++){ energia= Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica.. - //cout<<"folse"< energia[][1] e strip[][0] passaggio su ogni piano // integrale => energia totale nel calorimetro considerando il W // - // out[6] chi2,z,Etot,Pskip,flagmaxene,step - Float_t ordplane[44];//mi serve per la media troncat + // out[4] chi2,z,Etot,Pskip + + Float_t ordplane[44];//mi serve per la media troncata memset(ordplane,0,44*sizeof(Float_t)); for(Int_t ipla=0; ipla< 2*NPLA; ipla++) ordplane[ipla]=calorimetro[ipla][1]; //energia del piano //ordino tutte le energie dei piani in ordine crescente - Int_t ii=0; -// while( ii < NPLA-1 ){ - -// if(ordplane[ii+1] >= ordplane[ii]){ -// ii++; -// }else{ -// SWAP( &(ordplane[ii]), &(ordplane[ii+1])); -// ii=0; -// }; -// }; - - - //scelgo 4 piani minimo 'sensati' - Float_t sum4=0.; - Int_t pi=0; - cout<<"sum4="<0.7)&& sum4==0.){// (ordplane[pi] >50) - sum4=(ordplane[pi] + ordplane[pi+1] + ordplane[pi+2] + ordplane[pi+3]); - pi=2*NPLA; + //Int_t ii=0; + + Long64_t work[200]; + Int_t ind = 0; + Int_t l = 0; + Int_t RN = 0; + Float_t sum4 = 0.; + Float_t qm = 0.; + // + //Float_t qmt = ethr*0.8; // *0.9 + // + //Int_t uplim = TMath::Max(3,N); + // + while ( l < 4 && ind < 44 ){ + qm = TMath::KOrdStat(44,ordplane,ind,work); + if (qm >= 0.7 ){ + if ( l < 4 ){ + sum4 += qm; + RN++; + }; + l++; + if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4); + }; + ind++; }; - pi++; - }; - cout<<"sum4="< energia[][1] e strip[][0] passaggio su ogni piano // integrale=> energia totale nel calorimetro considerando il W // - // out[6] chi2,z,Etot,Pskip + // out[4] chi2,z,Etot,Pskip /*z se particella fosse al minimo*/ //energia1piano/mip corretta Float_t zmax = (sqrt(estremi[0][1]/spessore[2])); + if(zmax<31)zmax=zmax+1; /*calcolo Z ed E con loop sui vari elementi ed energie*/ - //Int_t step = 0; + Float_t zmin=1.; Float_t bestchitemp[4] = {0,0,0,0}; - // Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip + bestchi2[0]=10000.; bestchi2[1]=0.; bestchi2[2]=0.; @@ -857,7 +845,14 @@ Float_t zero=0.; //primo loop // energia ezero, zstart zstop Si attrav 1 piano energie piani out + + cout<<"inizio primo loop"<