void drawRaw(){ TELoss fAlphaSi; fAlphaSi.SetEL(1, 2.321); // density in g/cm3 fAlphaSi.AddEL(14., 28.086, 1); //Z, mass // mSi.SetZP(1., 1.); //protons fAlphaSi.SetZP(2., 4.); //alphas, Z, A fAlphaSi.SetEtab(100000, 200.); // ?, MeV calculate ranges fAlphaSi.SetDeltaEtab(300); TString cut[16]; TString cut_0[16]; TString cut1(""); TString cut2(""); TString cut3(""); TString cut4(""); TString hDraw,hname,hcut; TH1F *h0,*h1,*h2,*h[300]; TString space(" || "); for(Int_t i =0;i<16;i++) { cut[i].Form("SQY_L[%d]>130",i); cut_0[i].Form("SQY_L[%d]<130",i); if(i==0) { cut1 = cut[0]; cut3 = cut[0]; cut4 = cut_0[0]; } if(i<8) cut1 = cut1 + " || " + cut[i]; if(i==8) cut2 = cut[8]; if(i>7 && i<16) cut2 = cut2 + " || " + cut[i]; cut3 = cut3 + " || " + cut[i]; cut4 = cut4 + " && " + cut_0[i]; } //cout << "cut1 "<< cut1.Data() << endl; //cout << "cut2 "<< cut2.Data() << endl; //cout << "cut3 "<< cut3.Data() << endl; //cout << "cut4 "<< cut4.Data() << endl; TChain *ch = new TChain("AnalysisxTree"); //ch->Add("/media/analysis_nas/exp201804/calib/si_after/si-1_si-20_35cm_45deg_new1_00*.root"); //ch->Add("/media/analysis_nas/exp201804/calib/si_after/si-1_si-20_35cm_60deg_00*.root"); TChain *ch0 = new TChain("AnalysisxTree"); ch0->Add("/media/analysis_nas/exp201804/calib/si_after/si-1_si-20_35cm_0deg_new1_*.root"); // -0.4293 0.006653 0 /* TCanvas *c1 = new TCanvas("c1","raw data",1800,1000); c1->Divide(3,2); c1->cd(1); ch->Draw("NeEvent.SQY_L[7]*0.006653-0.4293 >> hy_0(200,0.435590,8)","","",1000000,1); h[0] = (TH1F*)gPad->GetPrimitive("hy_0"); h[0]->GetXaxis()->SetTitle("Edep in 1mm middle Y strip"); h[0]->GetXaxis()->SetTitleColor(2); h[0]->GetXaxis()->SetTitleOffset(0.75); h[0]->GetXaxis()->SetTitleSize(0.047); h[0]->GetXaxis()->CenterTitle(); c1->Update(); c1->cd(4); ch->Draw("NeEvent.SQY_L[15]*0.006501-0.4096 >> hy_15(200,0.435590,8)","","",10000000,1); h[1] = (TH1F*)gPad->GetPrimitive("hy_15"); h[1]->GetXaxis()->SetTitle("Edep in 1mm last Y strip"); h[1]->GetXaxis()->SetTitleColor(2); h[1]->GetXaxis()->SetTitleOffset(0.75); h[1]->GetXaxis()->SetTitleSize(0.047); h[1]->GetXaxis()->CenterTitle(); c1->Update(); c1->cd(2); ch->Draw("NeEvent.SQX_L[15]*0.006454+0.09241 >> hx_15(200,1.383210,8)","","",10000000,1); ch->Draw("NeEvent.SQX_L[0]*0.006705-0.001022 >> hx_0(200,1.339978,8)","","same",10000000,1); h[2] = (TH1F*)gPad->GetPrimitive("hx_15"); h[2]->GetXaxis()->SetTitle("Edep in 1mm middle(blue) and first/left(green) X strip"); h[2]->GetXaxis()->SetTitleColor(2); h[2]->GetXaxis()->SetTitleOffset(0.75); h[2]->GetXaxis()->SetTitleSize(0.047); h[2]->GetXaxis()->CenterTitle(); h[3] = (TH1F*)gPad->GetPrimitive("hx_0"); h[3]->SetLineColor(3); c1->Update(); c1->cd(5); ch->Draw("NeEvent.SQX_L[31]*0.006666+0.1119 >> hx_31(200,1.445100,8)","","",10000000,1); h[4] = (TH1F*)gPad->GetPrimitive("hx_31"); h[4]->GetXaxis()->SetTitle("Edep in 1mm last/right X strip"); h[4]->GetXaxis()->SetTitleColor(2); h[4]->GetXaxis()->SetTitleOffset(0.75); h[4]->GetXaxis()->SetTitleSize(0.047); h[4]->GetXaxis()->CenterTitle(); c1->Update(); c1->cd(3); ch->Draw("NeEvent.SQY_R[0] >> h20_0(1300,200,1500)","","",10000000,1); h[5] = (TH1F*)gPad->GetPrimitive("h20_0"); h[5]->GetXaxis()->SetTitle("Edep of 20mkm first/left strip"); h[5]->GetXaxis()->SetTitleColor(2); h[5]->GetXaxis()->SetTitleOffset(0.75); h[5]->GetXaxis()->SetTitleSize(0.047); h[5]->GetXaxis()->CenterTitle(); c1->Update(); c1->cd(6); ch->Draw("NeEvent.SQY_R[15] >> h20_15(1300,200,1500)","","",10000000,1); h[6] = (TH1F*)gPad->GetPrimitive("h20_15"); h[6]->GetXaxis()->SetTitle("Edep in 20mkm last/right strip"); h[6]->GetXaxis()->SetTitleColor(2); h[6]->GetXaxis()->SetTitleOffset(0.75); h[6]->GetXaxis()->SetTitleSize(0.047); h[6]->GetXaxis()->CenterTitle(); c1->Update(); TCanvas *c2 = new TCanvas("c2","selection of SQ20 data",1800,1000); c2->Divide(2,2); c2->cd(1); ch->Draw("NeEvent.SQY_R[15] >> h22_0-7(1300,200,1500)",cut1.Data(),"",10000000,1); h[7] = (TH1F*)gPad->GetPrimitive("h22_0-7"); h[7]->GetXaxis()->SetTitle("Edep in first half of 20mkm last/right strip (strips 0-7 of SQY_L were fired)"); h[7]->GetXaxis()->SetTitleColor(2); h[7]->GetXaxis()->SetTitleOffset(0.75); h[7]->GetXaxis()->SetTitleSize(0.047); h[7]->GetXaxis()->CenterTitle(); c2->Update(); c2->cd(2); ch->Draw("NeEvent.SQY_R[15] >> h22_8-15(1300,200,1500)",cut2.Data(),"",10000000,1); h[8] = (TH1F*)gPad->GetPrimitive("h22_8-15"); h[8]->GetXaxis()->SetTitle("Edep in second half of 20mkm last/right strip (strips 8-15 of SQY_L were fired)"); h[8]->GetXaxis()->SetTitleColor(2); h[8]->GetXaxis()->SetTitleOffset(0.75); h[8]->GetXaxis()->SetTitleSize(0.047); h[8]->GetXaxis()->CenterTitle(); c2->Update(); c2->cd(3); ch->Draw("NeEvent.SQY_R[15] >> h22_0-15(1300,200,1500)",cut3.Data(),"",10000000,1); h[9] = (TH1F*)gPad->GetPrimitive("h22_0-15"); h[9]->GetXaxis()->SetTitle("Edep in 20mkm last/right strip (strips 0-15 of SQY_L were fired)"); h[9]->GetXaxis()->SetTitleColor(2); h[9]->GetXaxis()->SetTitleOffset(0.75); h[9]->GetXaxis()->SetTitleSize(0.047); h[9]->GetXaxis()->CenterTitle(); c2->Update(); c2->cd(4); ch->Draw("NeEvent.SQY_R[15] >> h22_underThr(1300,200,1500)",cut4.Data(),"",10000000,1); h[10] = (TH1F*)gPad->GetPrimitive("h22_underThr"); h[10]->GetXaxis()->SetTitle("Edep in 20mkm last/right strip (strips 0-15 of SQY_L were NOT fired)"); h[10]->GetXaxis()->SetTitleColor(2); h[10]->GetXaxis()->SetTitleOffset(0.75); h[10]->GetXaxis()->SetTitleSize(0.047); h[10]->GetXaxis()->CenterTitle(); c2->Update(); */ /* /// calibration TF1* g1 = new TF1("g1", "gaus", 1040, 1125); g1->SetParLimits(0,1.,5000.); g1->SetParLimits(1,1040.,1125.); g1->SetParLimits(2,5,100.); TF1* g2 = new TF1("g2", "gaus", 1200, 1270); g2->SetParLimits(0,1.,5000.); g2->SetParLimits(1,1200.,1270.); g2->SetParLimits(2,1320.,1400.); TF1* g3 = new TF1("g3", "gaus", 1320, 1400); g3->SetParLimits(0,1.,5000.); g3->SetParLimits(1,1320.,1400.); g3->SetParLimits(2,5,100.); Double_t peakPosition,peakHight,x[3],y[3],e[3],temp; // energies after passing through 0.75 mkm dl y[0] = 4.7844; y[1] = 5.4895; y[2] = 6.0024; TList *functions; TPolyMarker *pm; TGraph *g[50]; TF1 *lin = new TF1("lin","[0]*x+[1]",0,1500); ofstream outcalfile; outcalfile.open("par_20_60.cal"); if (!outcalfile.is_open()) { Error("drawRaw", "Output file %s was not opened", "prepar_20.cal"); return 0; }//if ofstream outXfile; outXfile.open("positions_20_60.cal"); if (!outXfile.is_open()) { Error("drawRaw", "Output file %s was not opened", "positions_20.cal"); return 0; }//if Float_t dl[16]; ifstream myfile; TString line; Int_t count=0; myfile.open("/home/muzalevsky/work/temp/deadLayers.txt"); while (! myfile.eof() ){ line.ReadLine(myfile); if(line.IsNull()) break; sscanf(line.Data(),"%g", dl+count); count++; } for(Int_t i=0;i<16;i++) { cout << dl[i]<Divide(4,4); TCanvas *c5 = new TCanvas("c5","linfits",1800,1000); c5->Divide(4,4); for(Int_t i=0;i<16;i++) { cout << " new pad " << i+1 << endl; c4->cd(i+1); hDraw.Form("NeEvent.SQY_R[%d] >> h20_60_%d(1100,500,1600)",i,i); hname.Form("h20_60_%d",i); ch->Draw(hDraw.Data(),"NeEvent.trigger==3","",10000000,1); h[i] = (TH1F*)gPad->GetPrimitive(hname.Data()); h[i]->ShowPeaks(20,"", 0.5); c4->Update(); functions = h[i]->GetListOfFunctions(); pm = (TPolyMarker*)functions->FindObject("TPolyMarker"); //cout << *pm->GetY() << endl; for(Int_t j=0;j<3;j++) { x[j] = *(pm->GetX()+j); cout << x[j] << " "; } cout < x[n+1]) { temp = x[n]; x[n] = x[n+1]; x[n+1] = temp; } } } for(Int_t j=0;j<3;j++) { cout << x[j] << " "; } cout << endl; // recalculation energies using different effective dead layers for different strips cout << " shifted energies " << endl; for(Int_t j =0;j<16;j++) { for(Int_t k=0;k<3;k++) { e[k] = fAlphaSi.GetE(y[k],dl[j]); //cout << e[k] << " "; } //cout << endl; } g[i] = new TGraph(3,x,e); c5->cd(i+1); g[i]->SetMarkerSize(3); g[i]->Draw("AL*"); g[i]->Fit("lin"); c5->Update(); outcalfile << setprecision(4) << lin->GetParameter(1) << "\t" << lin->GetParameter(0) << "\t" << endl; outXfile << setprecision(4) << x[0] << "\t" << x[1] << "\t" << x[2] << endl; //for (Int_t m = 0; m < 3; m++) { // cout << x[m] << " "; //} //cout << endl; } */ /* cout << ch0->GetEntries() << " entr " << endl; Double_t peakPosition,peakHight,x[3],y[3],temp; y[0] = 4.7844; y[1] = 5.4895; y[2] = 6.0024; TList *functions; TPolyMarker *pm; TGraph *g[50]; TF1 *lin = new TF1("lin","[0]*x+[1]",0,1500); ofstream outcalfile; outcalfile.open("prepar_20.cal"); if (!outcalfile.is_open()) { Error("drawRaw", "Output file %s was not opened", "prepar_20.cal"); return 0; }//if ofstream outXfile; outXfile.open("positions_20.cal"); if (!outXfile.is_open()) { Error("drawRaw", "Output file %s was not opened", "positions_20.cal"); return 0; }//if TCanvas *c4 = new TCanvas("c4","raw 20 mkm data trigger=3",1800,1000); c4->Divide(4,4); //TCanvas *c5 = new TCanvas("c5","linfits",1800,1000); //c5->Divide(4,4); Double_t x0; for(Int_t i=0;i<16;i++) { cout << " new pad " << i+1 << endl; c4->cd(i+1); hDraw.Form("NeEvent.SQY_R[%d] >> h20_%d(1400,200,1600)",i,i); hname.Form("h20_%d",i); ch0->Draw(hDraw.Data(),"NeEvent.trigger==3","",7876702,1); h[i] = (TH1F*)gPad->GetPrimitive(hname.Data()); //x0 = h[i]->GetBinCenter(h[i]->GetMaximumBin()); gStyle->SetStatW(0.3); // Set width of stat-box (fraction of pad size) gStyle->SetStatH(0.4); c4->Update(); //functions = h[i]->GetListOfFunctions(); //pm = (TPolyMarker*)functions->FindObject("TPolyMarker"); //cout << *pm->GetY() << endl; //for(Int_t j=0;j<3;j++) { //x[j] = *(pm->GetX()+j); // cout << x[j] << " "; //} //cout < x[n+1]) { temp = x[n]; x[n] = x[n+1]; x[n+1] = temp; } } } //g[i] = new TGraph(3,x,y); //c5->cd(i+1); //g[i]->SetMarkerSize(3); //g[i]->Draw("AL*"); // g[i]->Fit("lin"); //c5->Update(); //outcalfile << setprecision(4) << lin->GetParameter(1) << "\t" << lin->GetParameter(0) << "\t" << endl; outXfile << setprecision(4) << x0 << endl; //for (Int_t m = 0; m < 3; m++) { // cout << x[m] << " "; //} //cout << endl; } */ //thickness determination /*Int_t SearchPeaks(const TH1 *hin, Double_t sigma, Option_t *option, const Int_t searchedpeaks) { //Function searching peaks in inputed TH1 spectrum and selects the peaks in the histogram. // // hin: // sigma: // option: // threshold: // searchedpeaks: TSpectrum sc; //by default for 100 peaks Int_t nopeaks = sc.Search(hin, sigma, "goff", fFitPeakThreshold); TString opt = option; opt.ToLower(); const Double_t tStep = 0.05; if (!nopeaks) { return 0; } if (opt.Contains("goff")) { return nopeaks; } TPolyMarker *pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker"); if (pm) { hin->GetListOfFunctions()->Remove(pm); delete pm; } pm = new TPolyMarker(nopeaks, sc.GetPositionX(), sc.GetPositionY()); hin->GetListOfFunctions()->Add(pm); pm->SetMarkerStyle(23); pm->SetMarkerColor(kRed); pm->SetMarkerSize(1.3); return nopeaks; } Int_t AculCalibration::PeaksFitting(TH1* hSpectrum, Option_t* option, Double_t sigmamin) { if (!hSpectrum) { cout<< "raw spec was not found " << endl; return 1; } Int_t dimension = hSpectrum->GetDimension(); if (dimension > 1) { Error("PeaksFitting", "Only implemented for 1-d histograms"); return 1; } TString opt = option; opt.ToLower(); if (!kRaNOPEAKS) { Error("PeaksFitting", "kRaNOPEAKS is set to zero; calibration spectrum must be set"); return 1; } const Int_t peaksNumber = SearchPeaks(hSpectrum, sigmamin, "", kRaNOPEAKS); if (peaksNumber != kRaNOPEAKS) { Info("PeaksFitting", "In histogram %s was found %d peaks", hSpectrum->GetName(), peaksNumber); return 1; } //should be optional output Info("PeaksFitting", "Number of peaks in %s: %d", hSpectrum->GetName(), peaksNumber); //working array for peaks, there are founded in accidental order Double_t peak[peaksNumber]; Double_t *peakPosition; Double_t *peakHight; TList *functions = hSpectrum->GetListOfFunctions(); TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker"); peakPosition = pm->GetX(); peakHight = pm->GetY(); for (Int_t i = 0; i < peaksNumber; i++) { Double_t fitMin = 0; Double_t fitMax = 0; Double_t fitStep = hSpectrum->GetXaxis()->GetBinWidth(0); // cout << fitStep << endl; // cout << fLowerPeakRelativeHight << "\t" << fUpperPeakRelativeHight << endl; //fitting range: //shift a range of fit and search for raw boarder of peak determined by fUpperPeakRelativeHight //maximum Int_t j = 0; Double_t currentHight = peakHight[i]; while ( currentHight > (peakHight[i]*fUpperPeakRelativeHight) ) { j++; fitMax = static_cast(peakPosition[i]) + j*fitStep; currentHight = hSpectrum->GetBinContent(hSpectrum->GetXaxis()->FindBin(fitMax)); } //minimum j = 0; currentHight = peakHight[i]; while ( currentHight > (peakHight[i]*fLowerPeakRelativeHight) ) { j++; fitMin = static_cast(peakPosition[i]) - j*fitStep; currentHight = hSpectrum->GetBinContent(hSpectrum->GetXaxis()->FindBin(fitMin)); } //fitting if (opt.Contains("gp")) { Info("PeaksFitting", "Option containing gp"); char fncname[20]; sprintf(fncname, "gaus_aux_%d", i); TF1 *gausAux = new TF1(fncname, "gaus", fitMin - 10, fitMax + 10); //pomocny gaus hSpectrum->Fit(fncname, "0 Q", "", fitMin - 15, fitMax + 15); //prvotni fitovani sprintf(fncname, "auto_gp_%d", i); TF1 *fitAuto = new TF1(fncname, "gaus(0) + pol0(3)", fitMin - 15, fitMax + 15); //fce pro automaticke fitovani fitAuto->SetParameter(0, gausAux->GetParameter(0)); //nastavovani parametru fitovaci fce fitAuto->SetParameter(1, gausAux->GetParameter(1)); fitAuto->SetParameter(2, gausAux->GetParameter(2)); hSpectrum->Fit(fncname, "0 R Q +", "", fitMin - 15, fitMax + 15); //dodelat zapis vsech fci hSpectrum->GetFunction(fncname)->ResetBit(TF1::kNotDraw); peak[i] = fitAuto->GetParameter(1); //zapis asi pozice v kanalech do pomocneho pole if (opt.Contains("V")) { Info("PeaksFitting", "Peak position is\t %4.2f \tresolution is \t %2.1f %%", fitAuto->GetParameter(1), 235*(fitAuto->GetParameter(2))/(fitAuto->GetParameter(1))); } } else { char fncname[20]; sprintf(fncname, "auto_g%d", i); TF1 *fitAuto = new TF1(fncname, "gaus", fitMin, fitMax); //fce pro automaticke fitovani // cout << fitMin << "\t" << fitMax << endl; // fitAuto->SetParameter(2, fitMax-fitMin); fitAuto->SetLineWidth(fFitFuncLineWidth); hSpectrum->Fit(fncname, "+ 0 R Q", "", fitMin - 1, fitMax + 1); // hSpectrum->GetFunction(fncname)->ResetBit(TF1::kNotDraw); hSpectrum->GetFunction(fncname)->InvertBit(TF1::kNotDraw); peak[i] = fitAuto->GetParameter(1); //zapis asi pozice v kanalech do pomocneho pole if (opt.Contains("v")) { Info("PeaksFitting", "Peak position is\t%4.2f\tresolution is \t%2.1f %%", fitAuto->GetParameter(1), 235*(fitAuto->GetParameter(2))/(fitAuto->GetParameter(1))); } }//else //end of fitting }//for over all analyzed peaks //peaks sorting Int_t j[peaksNumber]; TMath::Sort(peaksNumber, peak, j, kFALSE); fPeak.Set(peaksNumber); for (Int_t i = 0; i < peaksNumber; i++) { fPeak[i] = peak[j[i]]; //printf("\tPeak peak\t%f\n", fPeak[i]); } if (!opt.Contains("q") || opt.Contains("v")) { Info("PeaksFitting", "Control output:"); for (Int_t i = 0; i < peaksNumber; i++) { printf("\tPeak position is\t%f\n", fPeak[i]); } } // provest kontrolu pomerne polohy piku, // jestli jsou spatne, provest urcita opatreni, // napr. zapis daneho histogramu do souboru, // zapis do souboru s chybama, vypis na obrazovku, ... for (Int_t i = 0; i < peaksNumber; i++) { if ( !( (((1-fPeakPositionTolerance)*(fEnergy[0]/fEnergy[i])) < (fPeak[0]/fPeak[i])) && (((1+fPeakPositionTolerance)*(fEnergy[0]/fEnergy[i])) > (fPeak[0]/fPeak[i])) ) ) { //printf("\tPeaksFitt fEnergy\t%f\n", fEnergy[i]); if (fCalInformation && opt.Contains("writebad")) { fCalInformation->cd(); hSpectrum->Write(); } //return 2; } }//for return 0; } */ // reading input pars Float_t parYL1[16],parYL2[16]; ifstream myfile4; TString line4; Int_t count=-2; // myfile4.open("/home/muzalevsky/work/exp1803/data/siCal/SQY_L.cal"); myfile4.open("/home/muzalevsky/work/exp1803/data/siCal/presentPars/SQY_L.cal"); while (! myfile4.eof() ){ line4.ReadLine(myfile4); if(count < 0){ count++; continue; } if(line4.IsNull()) break; sscanf(line4.Data(),"%g %g", parYL1+count,parYL2+count); count++; } cout << endl << " pars for YL strips" << endl; for(Int_t i=0;i<16;i++) cout << parYL1[i] << " " << parYL2[i] << endl; // angles for DL in strips SQY_L Double_t alpha[16]; for(Int_t i=0;i<16;i++) { if(i<8) alpha[i] = atan(3.625*(7.5-i)/350.); if(i>7) alpha[i] = atan(3.625*(0.5+i-8)/350.); } cout << "angles for DL in strips SQY_L [Rad]" << endl; for(Int_t i=0;i<16;i++) { cout << alpha[i] << " "; } cout << endl; // angles for SQ20 strips Double_t betta[16]; for(Int_t i=0;i<16;i++) { if(i<8) betta[i] = atan(3.125*(7.5-i)/339.5); if(i>7) betta[i] = atan(3.125*(0.5+i-8)/339.5); } cout << "angles for SQ20 strips [Rad]" << endl; for(Int_t i=0;i<16;i++) { cout << betta[i] << " "; } cout << endl; // TH1* hSpectrum; Int_t peaksNumber,nopeaks; Double_t peak[peaksNumber]; Double_t *peakPosition; Double_t *peakHight; Double_t *peakPositionS;//sorted peaks Double_t *peakHightS; Double_t minPos,minH; TList *functions; TPolyMarker *pm; TSpectrum sc; //by default for 100 peaks Double_t fitMin, fitMax, fitStep; TF1 *fitAuto = new TF1("fitAuto", "gaus",0,10); Double_t upFitBorder,lowFitBorder,E[2],E0[2],nCh[2],th[2],thAv; upFitBorder = 0.5; lowFitBorder = 0.5; // initial alpha energies [MeV] E0[0] = 7.6869; E0[1] = 6.0024; Int_t m; TH2F *hTh = new TH2F("hTh","thickness distr",16,0,15,15,0,14); for(Int_t i =0;i<16;i++) { for(Int_t j=0;j<15;j++) { hTh->SetBinContent(i,j,0); } } TCanvas *c[16]; TString cName,cTitle; Int_t nPad,nCan; for(Int_t i=0;i<16;i++) { cName.Form("c%d",i); cTitle.Form("eDepostis in %d strip",i); c[i] = new TCanvas(cName.Data(),cTitle.Data(),1800,1000); c[i]->Divide(4,4); } for(Int_t i=0;i<16;i++) { // strips in SQ20 (0-15) for(Int_t nStr =0;nStr<15;nStr++) { // strips in SQY_L (0-14) cout << "#########" << endl << "SQ20 " << i << " SQY_L " << nStr << endl; //nCan = i/16; //nPad = i - nCan*16 + 1; c[i]->cd(nStr+1); /* } return; TCanvas *c6 = new TCanvas("c6","Edeposit distr",1800,1000); c6->Divide(4,4); for(Int_t i=0;i<4;i++) { c6->cd(i+1);*/ // обнуление fitMin = 0; fitMax = 0; peakPosition = NULL; peakHight = NULL; peakPositionS = NULL; peakHightS = NULL; // заполнение гистограмм через tree->Draw() hDraw.Form("NeEvent.SQY_L[%d]>>edep%d_%d(870,130,1000)",nStr,i,nStr); hname.Form("edep%d_%d",i,nStr); hcut.Form("NeEvent.SQY_R[%d]>200",i); ch0->Draw(hDraw.Data(),hcut.Data(),""); h[i] = (TH1F*)gPad->GetPrimitive(hname.Data()); hSpectrum = h[i]; if (!hSpectrum) { cout<< "raw spec was not found " << endl; return 1; } // подсчёт количества пиков nopeaks = sc.Search(hSpectrum, 8, "goff", 0.35); if (nopeaks<2) { cout << " less than 2 peaks were found " << endl; return; } cout << nopeaks << endl; //functions = hSpectrum->GetListOfFunctions(); //pm = (TPolyMarker*)functions->FindObject("TPolyMarker"); pm = new TPolyMarker(nopeaks, sc.GetPositionX(), sc.GetPositionY()); hSpectrum->GetListOfFunctions()->Add(pm); pm->SetMarkerStyle(23); pm->SetMarkerColor(kRed); pm->SetMarkerSize(1.3); peakPosition = pm->GetX(); peakHight = pm->GetY(); peakPositionS = new Double_t [pm->GetN()]; peakHightS = new Double_t [pm->GetN()]; for(Int_t a=0;aGetN();a++) { peakPositionS[a] = peakPosition[a]; peakHightS[a] = peakHight[a]; } //peakPositionS = pm->GetX(); //peakHightS = pm->GetY(); /* cout << "positions and heights " << endl; for(Int_t a=0;aGetN();a++) { cout << peakPosition[a] << " " << peakHight[a] << endl; }*/ for (Int_t b = 0; b < pm->GetN(); ++b){ for (Int_t a = 0; a < pm->GetN()-1; ++a){ if(peakPositionS[a] < peakPositionS[a+1]) { minPos = peakPositionS[a]; minH = peakHightS[a]; peakPositionS[a] = peakPositionS[a+1]; peakHightS[a] = peakHightS[a+1]; peakPositionS[a+1] = minPos; peakHightS[a+1] = minH; } } } /*cout << "sorted positions and heights " << endl; for(Int_t a=0;aGetN();a++) { cout << peakPositionS[a] << " " << peakHightS[a] << endl; }*/ fitStep = hSpectrum->GetXaxis()->GetBinWidth(0); //фитирование гауссом cout << "thickness!!" << endl; for(Int_t k = 0; k < 2; k++) { //cout << peakPosition[k] << " " << peakHight[k] << endl; //fitting range: //shift a range of fit and search for raw boarder of peak determined by fUpperPeakRelativeHight //maximum m = 0; Double_t currentHight = peakHightS[k]; while ( currentHight > (peakHightS[k]*upFitBorder) ) { m++; fitMax = static_cast(peakPositionS[k]) + m*fitStep; currentHight = hSpectrum->GetBinContent(hSpectrum->GetXaxis()->FindBin(fitMax)); } //minimum m = 0; currentHight = peakHightS[k]; while ( currentHight > (peakHightS[k]*lowFitBorder) ) { m++; fitMin = static_cast(peakPositionS[k]) - m*fitStep; currentHight = hSpectrum->GetBinContent(hSpectrum->GetXaxis()->FindBin(fitMin)); } //fitting fitAuto->SetRange(fitMin,fitMax); hSpectrum->Fit("fitAuto", "R+", "", fitMin, fitMax); //cout << fitAuto->GetParameter(1) << " " ; nCh[k] = fitAuto->GetParameter(1); E[k] = nCh[k]*parYL2[nStr] + parYL1[nStr]; //th[2],thAv th[k] = fAlphaSi.GetR(E0[k],E[k]); cout << nCh[k] << " " << E[k]<< " " << th[k] << endl; } // nopeaks == 2 thAv = (th[0] + th[1])/2.; thAv = ( thAv - 0.2 - (2.44/cos(alpha[nStr])) )/cos(betta[i]); hTh->SetBinContent(i, nStr, thAv); //SetBinContent(Int_t binx, Int_t biny, Double_t content) cout << "THICKNESS"<< endl << thAv << endl; /* peakPosition = pm->GetX(); peakHight = pm->GetY(); for(Int_t j=0;jUpdate(); } // nStr }// i // outfile TFile *fw = new TFile("pdf/thick.root", "RECREATE"); if (fw->IsOpen() == 0) { Error("CalculateCalibParameters", "File %s was not created and won't be processed", "pdf/thick.root"); return 1; } fw->cd(); hTh->Write(); fw->Close(); // printing canvases TString fileName; for(Int_t i=0;i<16;i++){ fileName.Form("pdf/SQ20_%d.png",i); c[i]->Print(fileName.Data()); } return; }