#include ClassImp(CaloBragg); //-------------------------------------- /* * Default constructor */ CaloBragg::CaloBragg(){ Clear(); }; CaloBragg::CaloBragg(PamLevel2 *l2p){ // Clear(); LoadParam(); // L2 = l2p; // if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n"); // OBT = 0; PKT = 0; atime = 0; // debug = false; usetrack = false; usepl18x = false; newchi2 = false; usenewBB = false; fzeta = -1.; // }; void CaloBragg::Clear(){ // ndf = 0; tr = 0; sntr = 0; // qtchi2 = 0.; // qtz = 0.; // qtetot = 0.; // 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)); Integrale=0.; for(Int_t l=0;l<44;l++){ calorimetro[l][0]=-1.; } }; void CaloBragg::Print(){ // if(!debug) Process(); // printf("========================================================================\n"); printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack); printf(" first plane: %f \n", estremi[0][0]); printf(" last plane: %f \n", estremi[1][0]); // printf(" chi 2 from truncated mean: %f \n", qtchi2); // printf(" Z from truncated mean %f: \n", qtz); // printf(" energy from truncated mean %f: \n", qtetot); // 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"); // }; void CaloBragg::Delete(){ Clear(); //delete this; }; void CaloBragg::Process(){ Process(-1); }; void CaloBragg::CleanPlanes(Float_t epiano[22][2]){ // return; Int_t hitplanes = 0; for (Int_t i = 0; i<22; i++){ for (Int_t j = 1; j>=0; j--){ if ( epiano[i][j] > 0.7 ) hitplanes++; }; }; Float_t lowlim = 0.85; Float_t dedxone = 0.; Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes; while ( dedxone < step1 ){ for (Int_t i = 0; i<22; i++){ for (Int_t j = 1; j>=0; j--){ if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j]; }; }; } if ( dedxone < 0.7 ){ for (Int_t i = 0; i<22; i++){ for (Int_t j = 1; j>=0; j--){ if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j]; }; }; } // // printf(" dedxone = %f step1 %f \n",dedxone,step1); Bool_t revulsera = false; Bool_t nullius = false; Int_t nulliferus = 0; 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.; } else { //x printf(" %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 ( nulliferus > 10 ) nullius = true; if ( nullius ) epiano[i][j] = 0.; }; }; } void CaloBragg::Process(Int_t ntr){ // if ( !L2 ){ printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); printf(" ERROR: CaloBragg variables not filled \n"); return; }; // Bool_t newentry = false; // if ( L2->IsORB() ){ if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){ newentry = true; OBT = L2->GetOrbitalInfo()->OBT; PKT = L2->GetOrbitalInfo()->pkt_num; atime = L2->GetOrbitalInfo()->absTime; sntr = ntr; }; } else { newentry = true; }; // if ( !newentry ) return; // tr = ntr; // if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); // Clear(); // // // Int_t view = 0; Int_t plane = 0; Int_t strip = 0; Float_t mip = 0.; Float_t epiano[22][2]; memset(epiano,0,22*2*sizeof(Float_t)); for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ // mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); // if ( !usepl18x && view==0 && plane==18 ) mip = 0.; // epiano[plane][view]+=mip; // // }; // this->CleanPlanes(*&epiano); // PamTrack *ptrack = 0; CaloTrkVar *track = 0; // if ( usetrack ){ if ( ntr >= 0 ){ ptrack = L2->GetTrack(ntr); if ( ptrack ) track = ptrack->GetCaloTrack(); } else { track = L2->GetCaloStoredTrack(ntr); //al momento e' vera solo questa riga }; // if ( !track && ntr >= 0 ){ printf(" ERROR: cannot find any track!\n"); printf(" ERROR: CaloBragg variables not completely filled \n"); return; }; } else { if ( ntr >= 0 ){ if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr); if ( debug ) printf(" ERROR: CaloBragg variables not completely filled \n"); return; }; }; // 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 // 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 }; }; /*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 /*settaggio della soglia per il loop sulla determinazione del piano di partenza */ 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 Long64_t work[200]; Int_t ind = 0; //Int_t l = 0; Int_t RN = 0; Float_t sum4 = 0.; Float_t qm = 0.; while ( RN < 4 && ind < 44 ){ qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work); if (qm >= 0.7 ){ if ( RN < 4 ){ sum4 += qm; RN++; }; }; ind++; }; // //sum4 /= (Float_t)RN; Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2]))); if(Zmean ==0.) Zmean=1.; if ( Zmean < 1. ) Zmean = 1.; /*trova primo e ultimo piano attraversati*/ Int_t p = 0;//contatore piani //per il primo parte da 0 e va in giu' while( estremi[0][1] <= 0. && p<(2*NPLA) ){ // era ==0 ma ricorda i problemi con Float == !!!!! // if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){ // if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >0.3)){ //0.7 mip = 70MeV soglia minima if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >Zmean*0.7)){ // 70% della MIP estremi[0][0]=p; estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV }; p++; }; //ultimo parte da 44 e sale p=43; while( (estremi[1][1] <= 0.) && (p>(int)estremi[0][0]) ){ if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ estremi[1][0]=p;// estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV }; p = p-1; }; // Float_t lastok = 0.; // 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 ( calorimetro[o][1] < 0.7 ) goback = true; // }; lastok = 0.; // 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 ( 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 ( 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++){ // printf(" integrale: calorimetro %f \n",calorimetro[pl][1]); //calcolo intergale in unita di spessori di silicio Integrale += calorimetro[pl][1] * MIP;//piano di silicio // se non e'il 1o dopo l'Y (tutti i pari) c'e' il W if(pl%2!=0){ //equival W in Si Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]); }; }; //Integrale=24000;//Integrale*1000; Integrale *= 1000.; /*z ed energia con media troncata*/ // 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 if ( debug ) this->Print(); if ( debug ) printf(" fine evento \n"); // }; 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= 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; // gStyle->SetLabelSize(0.04); gStyle->SetNdivisions(510,"XY"); // TString hid = Form("cCaloBragg"); TCanvas *tc = dynamic_cast(gDirectory->FindObject(hid)); if ( tc ){ // tc->Clear(); } else { tc = new TCanvas(hid,hid); // tc->Divide(1,2); }; // // TString thid = Form("hCaloBragg"); // TH2F *th = dynamic_cast(gDirectory->FindObject(thid)); // if ( th ) th->Delete(); // th->Clear(); // th->Reset(); // } else { // th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.); // th->SetMarkerStyle(20); // }; // tc->cd(); 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->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->SetMarkerStyle(20); th3->SetMarkerColor(kBlue); // tc->cd(1); // // 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); tc->cd(); // 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"); tc->Modified(); tc->Update(); // gStyle->SetLabelSize(0); gStyle->SetNdivisions(1,"XY"); // }; void CaloBragg::LoadParam(){ // 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] = -1.; elem[3][3] = -1.; elem[3][4] = -1.; elem[3][5] = -1.; elem[3][6] = -1.; elem[4][0] = 11.0093; //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 NPLA = 22; NCHA = 96; nView = 2; AA = 0.96;//mm larghezza strip ADIST = 80.5;//mm distanza tra pad PIANO = 8.59;//mm distanza ySi = 0.38;//mm spessore silicio yW = 2.66;//mm spessore tungsteno rhoSi = 2.33;//g/cm3 densita' silicio rhoW = 19.3;//g/cm3 densita' tugsteno MIP = 0.106;//Mev g/cm2 energia al minimo nel silicio per 0.38 mm emin = 0.; //parametri bethe-bloch pigr = 3.1415; Na = 6.02e-23; ZA = 0.49; /*Z/A per Si*/ // 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 */ }; // void CaloBragg::conversione(){ // calcolo spessore Si attraverato in funzione dell'inclinazione // e conversione dello spessore di W in Si e correzione del valore // della Mip pe lo spessore effettivo // // in : evento // // out: out[0] = gcm2Si = spessore silicio attraversato nel piano // out[1] = WinSi = spessore equivalente in Si del W attraversato // out[2] = Mip = fattore conversione energia riscalato allo spessore attrversatonel piano Float_t SiCross=0.; Float_t WCross = 0.; Float_t ytgx = 0; Float_t ytgy = 0; Float_t a = 0.; /*silicio*/ ytgx = ySi * L2->GetCaloLevel2()->tanx[0]; ytgy = ySi * L2->GetCaloLevel2()->tany[0]; //lunghezza effettiva di silicio attraversata (mm) SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy)); spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2 /*tungsteno*/ ytgx = yW * L2->GetCaloLevel2()->tanx[0]; ytgy = yW * L2->GetCaloLevel2()->tany[0]; //rapporto tra rilasci energetici nei due materiali 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; 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, Float_t II){ //rilascio energetico con bethe bloch con correzioni //in: x: g/cm2 // z: carica // Mass: Massa uma // Ene: energia particella MeV//tolta // gam: (etot/massa) // Bet: rad((g2-1)/g2) // //out: energia rilasciata MeV Float_t eta =0.; Float_t Wmax =0.; 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(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(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); Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.); *out =Energia;//out };//end Bethebloch 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 // Z => carica // Etot => energia perticella // // out: energia persa Float_t Q=0.; Float_t v=0.; Float_t gam=0.; Float_t Bet=0.; Float_t dEP=0.; // gamma // Mass = A * MassP; /*in Mev/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.)))); // 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); //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; };//end ELOSS 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 // E0 =>energia // spess2[3] => conversione spessore Si, Si in W, mip // primo => posizione primo piano attraversato // // out: array[44] =>rilasci energetici calcolati per ogni piano[44] dopo il primo(estremi[0][0]) Float_t dE=0.; //energia rilasciata Float_t Ezero= *E0;//energia iniziale //azzero energia rilasciata sui piani memset(out, 0, 2*NPLA*sizeof(Float_t)); 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 , 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 = 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{ Ezero = Ezero -dE;//energia residua }; }; };//fine loop piani };//end Enetrack void CaloBragg::chiquadro(Float_t dE[], Float_t out[]){ // calcola chi2 tra energia calcolata e misurata // in : dE[44] =>energia calcolata // calo3[44][2]=> [0]strip attraversata [1]energia misurata per ogni piano // estr2 => array con primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1] // // out: array[3]=> (chi2; piani scartati consecutivi(79= >3 quindi frammentato); piani scartati totale) Float_t sum = 0.; Float_t PianoPrecedente=0.; Float_t badplane=0.; Float_t badplanetot=0.; Float_t w,wi; // 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.; //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; 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; } };//end chiquadro void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite, Int_t nostep = 1000){ // //loop su z ed energie per trovare miglior z (ed energia) //in: nloop => energia massima da provare (nloop x E0) // E0 => energia iniziale (intergale) // Zstart => minimo z da cui patire // Zlimite => z a cui fermarsi (z al minimo di ionizz sul 1o piano) // //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.; Float_t Stepint =(step)/(Float_t)nostep;//passo per il calcolo di energia 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][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=zmin; inucl 0. ){ Massa = elem[inucl-1][isotope]*MassP; //loop energia 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 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]=1000.;//valore chi2 per questo z a questa energia matrixchi2[inucl][isotope][iene2][1]=1000.;//energia per questo chi2 matrixchi2[inucl][isotope][iene2][2]=1000.;//piani saltati nel chi2 break; } }//fine loop energia isotope++; //incremento il contatore isotopi }//fine loop isotopi isotope=0; //riazzero il contatore isotopi }//fine loop z isotope=0;//non dovrebbe servire //Emi 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 } } isotope++; } isotope=0; } };//endloopze // void CaloBragg::mediatroncata(){ // //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore // //in: ordplane[44] => array con energia dei piani // // spess[3] => conversioni spessore di silicio, w, mip // // estr[2][2] => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1] // // calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano // // integrale => energia totale nel calorimetro considerando il W // // // // 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 // 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 ( RN < 4 && ind < 44 ){ // qm = TMath::KOrdStat(44,ordplane,ind,work); // if (qm >= 0.7 ){ // if ( RN < 4 ){ // sum4 += qm; // RN++; // }; // // l++; // // if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4); // }; // ind++; // }; // // // sum4 /= (Float_t)RN; // Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));//ma non e'/1?? // if(Zmean ==0.) Zmean=1.; // if ( Zmean < 1. ) Zmean = 1.; // // Zmean =round(Zmean); // // if(Zmean <1.) Zmean=1.; // // if(Zmean >0.)Zmean =round(Zmean); // //======== per i nuclei======= // if (Zmean >=2.){ // ind = 0; // RN = 0; // sum4 = 0.; // qm = 0.; // while ( RN < 4 && ind < 44 ){ // qm = TMath::KOrdStat(44,ordplane,ind,work); // if (qm >= (Zmean*Zmean)-Zmean*Zmean*0.2 ){ // if ( RN < 4 ){ // sum4 += qm; // RN++; // }; // }; // ind++; // }; // // // sum4 /= (Float_t)RN; // Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));//ma non e' /1?? // } // //calcolo energia migliore per Z trovato con media troncata // // Float_t zmin=Zmean; // Float_t zmin=round(Zmean); // bestchi2[0]=10000.; // bestchi2[1]=0.; // bestchi2[2]=0.; // bestchi2[3]=0.; // Float_t zero=0.; // // step energia zstart zstop // loopze(Integrale,zero,zmin,zmin); // qtchi2=bestchi2[0]; // qtz=bestchi2[1]; // qtetot=bestchi2[2]; // qtpskip=bestchi2[3]; // };//end mediatroncata void CaloBragg::Zdaloop(){ //calcolo Z con un loop su tutti i possibli Z ed energie //in: ordplane[44]=> array con energia dei piani // spess1[3]=> conversioni spessore di silicio, w e mip // estr3[2][2]=> primo[0][0] e ultimo[1][0] piano ed energie // calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano // integrale=> energia totale nel calorimetro considerando il W // // out[4] chi2,z,Etot,Pskip /*z se particella fosse al minimo*/ //energia1piano/mip corretta // Float_t zmax = round(sqrt(estremi[0][1]/spessore[2])); // if(zmax<31)zmax=zmax+1; /*calcolo Z ed E con loop sui vari elementi ed energie*/ Float_t zmin=1.; Float_t zmax=32.; Float_t bestchitemp[5] = {0,0,0,0,0}; 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 // loopze(Integrale,zero,zmin,zmax); //-> 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); if ( debug) printf("Zdaloop start Integrale %f , outene %f \n",Integrale,bestchi2[2]); //------------secondo loop ---------------------- 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.; 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,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