--- calo/flight/CaloBragg/src/CaloBragg.cpp 2009/12/14 14:37:58 1.7 +++ calo/flight/CaloBragg/src/CaloBragg.cpp 2011/01/28 10:40:44 1.8 @@ -33,10 +33,10 @@ // tr = 0; sntr = 0; -// qtchi2 = 0.; -// qtz = 0.; -// qtetot = 0.; -// qtpskip = 0.; + // qtchi2 = 0.; + // qtz = 0.; + // qtetot = 0.; + // qtpskip = 0.; lpchi2 = 0.; lpz = 0.; lpetot = 0.; @@ -61,10 +61,10 @@ 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 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(" energy from loop %f: \n", lpetot); @@ -83,6 +83,56 @@ 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 ){ @@ -133,6 +183,7 @@ // }; // + this->CleanPlanes(*&epiano); // PamTrack *ptrack = 0; CaloTrkVar *track = 0; @@ -164,8 +215,8 @@ 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 + 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 }; }; @@ -181,27 +232,27 @@ //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++; - }; + 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.; + 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*/ @@ -216,7 +267,8 @@ }; p++; }; - //ultimo parte da 44 e sale + + //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)){ @@ -226,9 +278,34 @@ 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); + // + }; + // }; + + /*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 @@ -236,13 +313,14 @@ Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]); }; }; - Integrale=24000;//Integrale*1000; + //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*/ - Zdaloop(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip + Zdaloop(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip if ( debug ) this->Print(); @@ -255,69 +333,71 @@ Process(); -// 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 + // 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 - Float_t sp= spessore[0]*spessore[1]; - for(Int_t i=0;i<44;i++)Depth[i]=i*sp; + 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); -// }; - // - 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); - // - for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]); - th3->Draw(); - th2->Draw("same"); + 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]); + th3->Draw(); + th2->Draw("same"); - tc->Modified(); - tc->Update(); + tc->Modified(); + tc->Update(); // gStyle->SetLabelSize(0); @@ -364,7 +444,7 @@ elem[31] = 72.61; //Ge 32 -//parametri calorimetro + //parametri calorimetro NPLA = 22; NCHA = 96; nView = 2; @@ -430,7 +510,7 @@ 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); + a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664); //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664); // (g/cm2W)/(g/cm2Si) spessore[1] = a; @@ -469,13 +549,13 @@ Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass)); lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi); - // Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA* SQ(z)/SQ(Bet) * lg; + // 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.); 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.); + Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.); *out =Energia;//out @@ -619,7 +699,8 @@ }; //meno peso ai piani con rilasci maggiori di 1000 MIP - if(calorimetro[ipla][1] > 1000) wi=0.5; + // 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)); @@ -636,6 +717,7 @@ //da Emi if(badplane > 2){ + // printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n"); out[1] =79.; break; }; @@ -649,8 +731,8 @@ -void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite){ -// +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) @@ -668,7 +750,7 @@ Float_t Massa = 0.; - Float_t Stepint =(step)/1000.;//passo per il calcolo di energia + Float_t Stepint =(step)/(Float_t)nostep;//passo per il calcolo di energia Float_t energia =0.;//energia del loop @@ -683,6 +765,8 @@ Float_t matrixchi2[colmax][rowmax][3]; memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t)); + Int_t imin = 1-nostep/2; + Int_t imax = nostep/2; //loop elementi for(Int_t inucl=(int)(Zstart); inucl0.)){ bestchi2[0]= matrixchi2[nu][en][0];// chi2 bestchi2[1]= (Float_t)nu; // z @@ -845,12 +932,13 @@ /*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; + // 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[4] = {0,0,0,0}; bestchi2[0]=10000.; @@ -860,8 +948,10 @@ Float_t zero=0.; //------------primo loop ---------------------- // energia ezero, zstart zstop - loopze(Integrale,zero,zmin,zmax); - + // loopze(Integrale,zero,zmin,zmax); + loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50); + // loopze(Integrale*2.,Integrale/100.,zmin,zmax); + // printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]); //------------secondo loop ---------------------- for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i]; @@ -871,9 +961,12 @@ bestchi2[3] = 0.;//riazzero Float_t step = bestchitemp[2];// - zero=0; // qualsiasi altro valore peggiora le cose - zmin=zmax=bestchitemp[1]; - loopze(step,zero,zmin,zmax); // + 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); // //chi2,z,Etot,Pskip