www.eprace.edu.pl » minos-neutrina » Dodatki » Sposób użycia algorytmu - metoda cięć

Sposób użycia algorytmu - metoda cięć

W sposób bardziej szczegółowy, w postaci kodu programu, przedstawiony zostanie algorytm do selekcji przypadków metodą cięć (patrz: ANALIZA). Program selekcja.C został napisany w środowisku root na bazie loon'a (patrz: Dodatek A).

 #include <stdio.h>
 void selekcja() //selekcja.C
 {
 //tutaj następuje deklaracja podstawowych zmiennych oraz używanych
 //funkcji, np.:
 void ustawienia() //funkcja odpowiada za ustawienia wygladu histogramów
 {
 gStyle->SetOptStat(10000011); //bez statystyki
 //pozostałe parametry: styl, rozmiar, kolor, tło etc.
 ...
 }
 //przykład funkcji odpowiedzialnej za rysowanie histogramu dowolnej zmiennej
 //zapisanej jako TEST_NC2 dla przypadków NC i TEST_CC2 dla CC:
 void hist_TEST2()
 {
 TCanvas* TEST2 = new TCanvas (ćanvas-testowy","TEST",70,300,1200,600);
 TEST_NC2->SetFillColor (3); //kolor wypełnienia
 TEST_NC2->SetFillStyle (3001); //styl wypełnienia
 TEST_CC2->SetFillColor (4);
 TEST_CC2->SetFillStyle (3001);
 gStyle->SetOptStat(0); //indywidualne ustawienia wygladu
 gStyle->SetTitleOffset(1.2,"x");
 gStyle->SetTitleOffset(1.2,"y");
 TEST2->Divide (1,1); //podział
 TEST2->cd(1);
 TEST_CC2->Draw();
 TEST_NC2->Draw("SAME"); //rysowanie własciwe
 TEST2->Update(); //dwie niezbędne funkcje
 TEST2->Modified();
 TLegend *legend = new TLegend (0.8, 0.88, 0.88, 0.7); //legenda
 legend->AddEntry (TEST_NC2, "NC", "f");
 legend->AddEntry (TEST_CC2, "CC", "f");
 legend->Draw (); //rysuj legendę
 }
 //przykład funkcji rysujacej czystosc i efektywnosc:
 void Pur_Eff(TH1F* CCjakoCC, TH1F* CCjakoNC, TH1F* NCjakoNC, TH1F* NCjakoCC,
 TH1F* CCall, TH1F* NCall, float itCCjakoCC, float itCCjakoNC,
 float itNCjakoNC, float itNCjakoCC, float iIleCC, float iIleNC)
 {
 double r=0, N=0; //niezbędne do liczenia błędów
 TH1F *hPurCC = new TH1F("PurCC","PurCC",50,0,10);//wlasciwe histogramy
 TH1F *hEffCC = new TH1F("EffCC","EffCC",50,0,10);//czystosci i
 TH1F *hPurNC = new TH1F("PurNC","PurNC",50,0,10);//efektywnosci
 TH1F *hEffNC = new TH1F("EffNC","EffNC",50,0,10);
 TH1F *temp = new TH1F("temp","temp;Energia [GeV];CC_jako_CC + NC_jako_CC",
 50,0,10);
 //obliczenia na podstawie wzorow na Pur i Eff:
 temp -> Add(CCjakoCC,NCjakoCC,1,1);
 hPurCC -> Divide(CCjakoCC,temp,1,1);
 hEffCC -> Divide(CCjakoCC, CCall, 1, 1);
 TH1F *temp1 = new TH1F("temp1","temp1;Energia [GeV];CC_jako_NC + NC_jako_NC",
 50,0,10);
 temp1 -> Add(NCjakoNC,CCjakoNC,1,1);
 hPurNC -> Divide(NCjakoNC,temp1,1,1);
 hEffNC -> Divide(NCjakoNC, NCall, 1, 1);
 //koniec wlasciwych obliczen - mamy histogramy Pur i Eff
 //rysowanie Eff NC:
 TCanvas *c0 = new TCanvas("effNC", "Eff - efektywnosc NC", 0, 10, 800, 600);
 hEffNC->SetMinimum(0); //ustawianie zakresu osi Y
 hEffNC->SetMaximum(1);
 //wyznaczanie slupkow bledow:
 Int_t nbinsx = hEffNC->GetNbinsX(); //liczba binow histogramu (=50)
 r=0; N=0;
 for (Int_t binx=1; binx<=nbinsx; binx++) //petla po binach
 {
 r=hEffNC->GetBinContent(binx); //odczytuje wartosc binu nr 'binx'
 N=NCall->GetBinContent(binx);
 //wstawienie wartosci policzonego bledu do danego binu:
 if (N!=0) hEffNC->SetBinError(binx,sqrt(r*(1.0-r)/N));
 else hEffNC->SetBinError(binx,2);
 }
 hEffNC->Draw(); //rysowanie
 c0->Update();
 c0->Modified();
 hEffNC->SetTitle(""); //ustawienia
 hEffNC->SetXTitle("Energia [GeV]");
 hEffNC->SetYTitle("Eff - efektywnosc NC");
 hEffNC->SetStats(0); //brak legendy statystyki
 //rysowanie Pur CC:
 TCanvas *c1 = new TCanvas("purCC", "Pur - czystosc CC", 100, 110, 800, 600);
 hPurCC->SetMinimum(0);
 hPurCC->SetMaximum(1);
 r=0; N=0;
 nbinsx = hPurCC->GetNbinsX();
 for (Int_t binx=1; binx<=nbinsx; binx++)
 {
 r=hPurCC->GetBinContent(binx);
 //inne N niz w Eff!
 N=(CCjakoCC->GetBinContent(binx))+(NCjakoCC->GetBinContent(binx));
 if (N!=0) hPurCC->SetBinError(binx,sqrt(r*(1.-r)/N));
 else hPurCC->SetBinError(binx,2);
 }
 hPurCC->Draw();
 c1->Update();
 c1->Modified();
 hPurCC->SetTitle("");
 hPurCC->SetXTitle("Energia [GeV]");
 hPurCC->SetYTitle("Pur - czystosc CC");
 hPurCC->SetStats(0);
 //PurNC i EffCC wyznaczamy analogicznie
 ...
 //ustawienia typu punktow na wykresie
 hPurCC->SetMarkerStyle(22);
 hPurNC->SetMarkerStyle(23);
 hEffCC->SetMarkerStyle(22);
 hEffNC->SetMarkerStyle(23);
 hPurCC->SetMarkerColor(kBlue);
 hPurNC->SetMarkerColor(kRed);
 hEffCC->SetMarkerColor(kBlue);
 hEffNC->SetMarkerColor(kRed);
 }
 ...
 //dodawanie bibliotek
 gSystem->Load(łibCandNtupleSR.so");
 gSystem->Load(łibTruthHelperNtuple.so");
 gSystem->Load(łibMCNtuple.so");
 gSystem->Load(łibStandardNtuple.so");
 gROOT->LoadMacro("MyClass.C"); //omówione na końcu dodatku B
 MyClass t; //podstawowa klasa odpowiedzialna za NTuple
 //ustawienia wyswietlania histogramów
 gROOT->SetStyle("Plain");
 gROOT->ForceStyle();
 ustawienia(); //dokładne ustawienia histogramów
 Long64_t nentries = t.fChain->GetEntries();
 //pobiera całkowita liczbę przypadków z
 //wczytanej paczki danych. Inaczej: liczba
 //iteracji w pętli głównej programu
 //potrzebne do obliczeń Pur i Eff
 int itCCjakoCC = 0;
 int itCCjakoNC = 0;
 int itNCjakoNC = 0;
 int itNCjakoCC = 0;
 int iIleCC=0,iIleNC=0;
 //deklaracja dalszych używanych zmiennych pomocniczych
 ...
 //ustawienia zakresów i nazw histogramów w zależnosci
 //od np. analizowanej aktualnie zmiennej
 ...
 t.fChain.SetBranchStatus("*",0);
 //poczatkowo niezbędne jest wyłaczenie wszystkich
 //zmiennych w NTuplu z uwagi na pamięć komputera
 //właczanie wybranych gałęzi czy zmiennych
 t.fChain.SetBranchStatus("mc*",1); //cała gałaz MC
 t.fChain.SetBranchStatus("evthdr*",1);
 t.fChain.SetBranchStatus(śtp*",1);
 ...
 //deklaracje histogramów. Przykład:
 TH1D* TEST_NC = new TH1D ("TEST_NC",napisTEST,50.,zakres_od,zakres_do);
 TH1D* TEST_CC = new TH1D ("TEST_CC",napisTEST,50.,zakres_od,zakres_do);
 ...
 //w przypadku histogramów pomocnych w rysowaniu Pur i Eff należy zwrócić
 //uwagę na wywołanie funkcji 'Sumw2' niezbędnej przy rysowaniu błędów:
 TH1F *NCall = new TH1F("NCall","NCall;Energia [GeV]; wszystkie NC",50,0,10);
 TH1F *CCall = new TH1F("CCall","CCall;Energia [GeV]; wszystkie CC",50,0,10);
 TH1F *CCjakoCC = new TH1F("CCjakoCC","CCjakoCC;Energia [GeV];CC jako CC",50,0,10);
 TH1F *CCjakoNC = new TH1F("CCjakoNC","CCjakoNC;Energia [GeV];CC jako NC",50,0,10);
 TH1F *NCjakoNC = new TH1F("NCjakoNC","NCjakoNC;Energia [GeV];NC jako NC",50,0,10);
 TH1F *NCjakoCC = new TH1F("NCjakoCC","NCjakoCC;Energia [GeV];NC jako CC",50,0,10);
 CCjakoCC->Sumw2();
 CCjakoNC->Sumw2();
 NCjakoCC->Sumw2();
 NCjakoNC->Sumw2();
 NCall->Sumw2();
 CCall->Sumw2();
 ...
 //niezbędny komunikat:
 cout<<endl<<"Calkowita liczba przypadkow="<<nentries<<endl;
 cout<<"Obliczenia... prosze czekac..."<<endl;
 //ROZPOCZYNA SIĘ GŁÓWNA PĘTLA PO WSZYSTKICH PRZYPADKACH:
 for (Int_t jentry = 0; jentry < nentries; jentry++)
 {
 Long64_t ientry = t.LoadTree(jentry);
 if (ientry < 0) break;
 nb = t.fChain->GetEntry(jentry); nbytes += nb;
 int iaction = t.mc_iaction[0]; //= 1 dla CC; =0 dla NC
 float enGeV = (double)t.evthdr_ph_sigcor*((double)0.000082);
 //zamiana energii na GeV
 //deklaracja dodatkowych zmiennych
 ...
 //tworzenie prostych histogramów ogólnych
 ...
 //przykład generowania histogramu rozkładu energii w zależnosci od CC/NC:
 if (enGeV>0)
 {
 if (iaction==1) {jeden++; histogram_CC -> Fill(enGeV);} //dla CC
 if (iaction==0) {zero++; histogram_NC -> Fill(enGeV);} //dla NC
 }
 //rozkład dowolnej zmiennej przed selekcja:
 if (enGeV>0)
 {
 if (iaction==0) TEST_NC2 -> Fill(zmienna);
 if (iaction==1) TEST_CC2 -> Fill(zmienna);
 }
 //obliczanie różnych wielkosci niezbędnych do selekcji, przykładowo:
 for (int i=0; i<t.evthdr_ntrack; i++)
 {
 if (t.trk_range[i]>maks1) maks1=t.trk_range[i]; //zasięg wiodacego toru
 if (t.trk_ph_sigcor[i]>maks2) maks2=t.trk_ph_sigcor[i];
 //sygnał wiodacego toru
 }
 for (int i=0; i<t.evthdr_nshower; i++)
 {
 if (t.shw_ph_sigcor[i]>maks3) maks3=t.shw_ph_sigcor[i];
 //sygnał wiodacej kaskady
 if (t.shw_ndigit[i]>maks4) maks4=t.shw_ndigit[i];
 }
 for (int i=0;i<t.evthdr_nstrip;i++)
 //sygnał z pierwszych 25% zapalonych płaszczyzn
 {
 if ((t.stp_plane[i]>t.evthdr_plane_beg) &&
 (t.stp_plane[i]<=t.evthdr_plane_beg+dlugosc*25./100.))
 sygnal+=t.stp_ph0_sigcor[i]+t.stp_ph1_sigcor[i];
 }
 ...
 //definicje 3 zmiennych do ciec:
 double zmienna1=0; //zabezpieczenie przed dzieleniem przez zero gdy maks3=0;
 if (maks3>0) zmienna1=sygnal/maks3;
 double zmienna2=t.evthdr_plane_n;
 double zmienna3=maks2;
 //SELEKCJA WŁASCIWA:
 if (t.evthdr_ph_sigcor>0) //warunek na niezerowa energię
 {
 if (iaction==1) ileCCwszystkich++; //oblicza ile jest przypadków CC/NC
 if (iaction==0) ileNCwszystkich++;
 ...
 //własciwy algorytm selekcji (patrz: 'Algorytm selekcji')
 if (zmiennaN>=ciecieN) //na podstawie danego algorytmu
 {
 //przykładowe histogramy i zmienne, gdy przypadek
 //jest zaklasyfikowany np. jako NC:
 hist_NC -> Fill(zmiennaX); licznik_NC++; typ_czastki=0;
 //opcjonalnie do wstawienia generacje histogramów 'N-1' lub
 //histogramów po wszystkich poprzednich cięciach
 ...
 //także innych użytecznych zmiennych-liczników, np.:
 ile_przypadkow_NC_po_N_cieciu++;
 ...
 goto wyjscie;
 }
 else
 if (zmiennaM>=ciecieM)
 {
 //analogicznie jak poprzednio
 ...
 }
 else
 //po kolei wszystkie zmienne i cięcia
 ...
 else {typ_czastki=-1; goto wyjscie;} //UNKNOWN
 }
 wyjscie:
 //w tym miejscu jestesmy już po selekcji i rozpatrywany
 //przypadek numer 'jentry' mamy zaklasyfikowany jako
 //NC/CC na podstawie zmiennej 'typ_czastki'
 //rozkład zmiennej po selekcji:
 if (enGeV>0)
 {
 if (typ_czastki==0) {TEST_NC -> Fill(zmienna);}
 if (typ_czastki==1) {TEST_CC -> Fill(zmienna);}
 }
 //tworzenie pozostałych histogramów - 'N-1',
 //po poprzednich cięciach i innych
 ...
 //liczenie wartosci do Pur i Eff
 //- czystosci i efektywnosci:
 if ( iaction == 1 )
 {
 iIleCC++;
 CCall->Fill(enGeV);
 }
 else
 if ( iaction == 0 )
 {
 iIleNC++;
 NCall->Fill(enGeV);
 }
 if (t.evthdr_ph_sigcor>0) //dla niezerowej energii
 {
 if ( iaction == 1 || iaction == 0 )
 {
 if (iaction == 1) //CC
 {
 if (czastka == 1)
 {
 itCCjakoCC++;
 CCjakoCC->Fill(enGeV);
 }
 else
 if (czastka == 0) //NC
 {
 itCCjakoNC++;
 CCjakoNC->Fill(enGeV);
 }
 }
 else //czyli iaction==0
 {//tutaj NC
 if (czastka == 0)
 {
 itNCjakoNC++;
 NCjakoNC->Fill(enGeV);
 }
 else
 if (czastka == 1)
 {
 itNCjakoCC++;
 NCjakoCC->Fill(enGeV);
 }
 }
 }
 }
 ...
 } //koniec głównej pętli po przypadkach
 //liczenie Pur i Eff na podstawie wyliczonych uprzednio wartosci
 float PurCC = (float) itCCjakoCC / (itCCjakoCC + itNCjakoCC);
 float EffCC = (float) itCCjakoCC / iIleCC;
 float PurNC = (float) itNCjakoNC / (itNCjakoNC + itCCjakoNC);
 float EffNC = (float) itNCjakoNC / iIleNC;
 //wypisanie na ekranie danych liczbowych wybranych
 //interesujących zmiennych
 ...
 //rysowanie histogramów jako wywolania wczesniej
 //zadeklarowanych funkcji, np.:
 hist_TEST2(); //rozklad dowolnej zmiennej
 Pur_Eff(CCjakoCC,CCjakoNC,NCjakoNC,NCjakoCC,CCall,NCall,itCCjakoCC,
 itCCjakoNC,itNCjakoNC,itNCjakoCC,iIleCC,iIleNC); //Pur i Eff
 ...
 }
 

Do prawidłowego działania programu niezbędne są dodatkowe dwa pliki: MyClass.h oraz MyClass.C. Generuje się je automatycznie w celu obsługi plików typu SNTP (zobacz: www.physics.ox.ac.uk/minos/software/oo/html/reco_8h-source.html). Dzięki temu poszczególne zmienne dostępne są w programie głównym selekcja.C bezpośrednio poprzez zastosowanie konwencji: klasa.galaz_zmienna, np.:

 float liczba_fotoelektronow = t.evthdr_ph_pe;
 

W pliku MyClass.h dokonuje się dodawania kolejnych plików zawierających przypadki oddziaływań instrukcją:

 MyClass::MyClass(TTree *tree) //konstruktor
 {
 if (tree == 0) //łaczenie się z plikami, jeśli
 //jeszcze nie były zdefiniowane
 {
 (...)
 TChain * chain = new TChain("NtpSt",""); \\tworzy łańcuch
 cout<<"Wczytuje pliki..."<<endl;
 //Następuje wczytywanie po kolei plików do łańcucha 'chain'
 chain->Add("/sciezka/dostepu/plik1.root");
 chain->Add("/sciezka/dostepu/plik2.root");
 chain->Add("/sciezka/dostepu/plik3.root");
 (...)
 tree = chain;
 }
 Init(tree); //inicjuje drzewo
 }
 

W metodzie cięć użyto następujących plików Monte Carlo:

 //Daleki detektor, właczone pole B, wiazka niskoenergetyczna LE-10-185
 f21011002_0000_L010185N_D00.sntp.cedar.root
 f21011004_0000_L010185N_D00.sntp.cedar.root
 f21011005_0000_L010185N_D00.sntp.cedar.root
 f21011006_0000_L010185N_D00.sntp.cedar.root
 f21011008_0000_L010185N_D00.sntp.cedar.root
 f21011009_0000_L010185N_D00.sntp.cedar.root
 f21011011_0000_L010185N_D00.sntp.cedar.root
 f21011012_0000_L010185N_D00.sntp.cedar.root
 f21011013_0000_L010185N_D00.sntp.cedar.root
 f21011014_0000_L010185N_D00.sntp.cedar.root
 f21011015_0000_L010185N_D00.sntp.cedar.root
 f21011017_0000_L010185N_D00.sntp.cedar.root
 f21011019_0000_L010185N_D00.sntp.cedar.root
 f21011020_0000_L010185N_D00.sntp.cedar.root
 f21011021_0000_L010185N_D00.sntp.cedar.root
 f21011022_0000_L010185N_D00.sntp.cedar.root
 f21011025_0000_L010185N_D00.sntp.cedar.root
 

komentarze

Przeszukaj ePrace

Copyright © 2008-2010 EPrace oraz autorzy prac.