Sposób użycia algorytmu - metoda RS
Algorytm do selekcji oddziaływań CC i NC metodą wieloparametryczną przy użyciu Range Searching składa się z dwóch oddzielnych programów. Pierwszy z nich, prostszy, służy do utworzenia wzorca - dwóch plików binarnych zawierających trzy współrzędne wszystkich wczytywanych po kolei przypadków (patrz: Metoda RS). Poniżej znajduje się pełny kod wspomnianego programu. Kod programu właściwego znajduje się dalej.
#include <stdio.h>
void tworzy_drzewa() //tworzy_drzewa.C
{
double maks2=0, maks3=0;
gROOT->Reset();
gSystem->Load(łibCandNtupleSR.so");
gSystem->Load(łibTruthHelperNtuple.so");
gSystem->Load(łibMCNtuple.so");
gSystem->Load(łibStandardNtuple.so");
gROOT->LoadMacro("MyClass.C"); //pliki MyClass.* niezbędne!
MyClass t; //podstawowa klasa do obsługi zmiennych
Long64_t nentries = t.fChain->GetEntries(); //liczba przypadków
t.fChain.SetBranchStatus("*",0); //wylacza wszystkie zmienne
t.fChain.SetBranchStatus("evthdr*",1); //wlacza sama galaz 'evthdr'
t.fChain.SetBranchStatus("mc.iaction",1); //wlacza zmienna iaction
\\plik zawierajacy tlo (NC):
Float_t a, b, c;
TFile *foutB = new TFile("background.root","RECREATE");
TTree *outTreeB = new TTree("background","Background Sample");
outTreeB->Branch(ą",&a,"F");
outTreeB->Branch("b",&b,"F");
outTreeB->Branch(ć",&c,"F");
\\plik zawierajacy sygnal (CC):
Float_t aa, bb, cc;
TFile *foutP = new TFile("physics.root","RECREATE");
TTree *outTreeP = new TTree("physics","Physics Sample");
outTreeP->Branch(ą",&aa,"F");
outTreeP->Branch("b",&bb,"F");
outTreeP->Branch(ć",&cc,"F");
cout<<endl<<Źapisuje "<<nentries<<" przypadkow do plikow..."<<endl;
for (Int_t jentry = 0; jentry < nentries; jentry++)
{
maks2=0; maks3=0;
int iaction = t.mc_iaction[0]; //CC czy NC
//obliczanie niezbędnych wielkosci:
for (int i=0; i<t.evthdr_ntrack; i++)
{
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
}
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 wejsciowych:
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;
if (t.evthdr_ph_sigcor>0) //warunek na niezerowa energię
{
if (iaction==0) //NC = tło
{
a=zmienna1; //a,b,c - wspolrzedne przypadku NC
b=zmienna2;
c=zmienna3;
outTreeB->Fill(); //zapis 3 zmiennych do pamięci
}
if (iaction==1) //CC = sygnał
{
aa=zmienna1; //aa,bb,cc - wspolrzedne przypadku CC
bb=zmienna2;
cc=zmienna3;
outTreeP->Fill(); //zapis 3 zmiennych do pamięci
}
}
} //koniec pętli po przypadkach
//zapisywanie i zamykanie plików
foutB->Write();
foutB->Close();
foutP->Write();
foutP->Close();
cout<<"Pliki wzorcowe NC i CC zapisane!"<<endl;
}
Powyższy program tworzy dwa pliki wzorcowe na podstawie plików zawierających przypadki oddziaływań Monte Carlo (ścieżki w MyClass.h - patrz Dodatek B):
//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
f21011026_0000_L010185N_D00.sntp.cedar.root
f21011027_0000_L010185N_D00.sntp.cedar.root
f21011029_0000_L010185N_D00.sntp.cedar.root
f21011030_0000_L010185N_D00.sntp.cedar.root
f21011031_0000_L010185N_D00.sntp.cedar.root
f21011034_0000_L010185N_D00.sntp.cedar.root
f21011036_0000_L010185N_D00.sntp.cedar.root
f21011037_0000_L010185N_D00.sntp.cedar.root
f21011038_0000_L010185N_D00.sntp.cedar.root
f21011039_0000_L010185N_D00.sntp.cedar.root
f21011040_0000_L010185N_D00.sntp.cedar.root
f21011042_0000_L010185N_D00.sntp.cedar.root
f21011043_0000_L010185N_D00.sntp.cedar.root
f21011044_0000_L010185N_D00.sntp.cedar.root
f21011045_0000_L010185N_D00.sntp.cedar.root
f21011046_0000_L010185N_D00.sntp.cedar.root
f21011049_0000_L010185N_D00.sntp.cedar.root
f21011052_0000_L010185N_D00.sntp.cedar.root
f21011053_0000_L010185N_D00.sntp.cedar.root
f21011055_0000_L010185N_D00.sntp.cedar.root
f21011056_0000_L010185N_D00.sntp.cedar.root
f21011058_0000_L010185N_D00.sntp.cedar.root
f21011060_0000_L010185N_D00.sntp.cedar.root
f21011062_0000_L010185N_D00.sntp.cedar.root
f21011065_0000_L010185N_D00.sntp.cedar.root
f21011066_0000_L010185N_D00.sntp.cedar.root
f21011068_0000_L010185N_D00.sntp.cedar.root
f21011069_0000_L010185N_D00.sntp.cedar.root
f21011070_0000_L010185N_D00.sntp.cedar.root
f21011079_0000_L010185N_D00.sntp.cedar.root
f21011080_0000_L010185N_D00.sntp.cedar.root
f21011083_0000_L010185N_D00.sntp.cedar.root
f21011084_0000_L010185N_D00.sntp.cedar.root
f21011086_0000_L010185N_D00.sntp.cedar.root
f21011087_0000_L010185N_D00.sntp.cedar.root
f21011088_0000_L010185N_D00.sntp.cedar.root
f21011089_0000_L010185N_D00.sntp.cedar.root
Drugą, główną częścią algorytmu do selekcji oddziaływań CC i NC metodą RS jest program autorstwa B. Koblitza i T. Carliego [18]. Jest on dostępny w internecie pod adresem www.desy.de/~koblitz (wersja pełna range_searching.tgz). Składa się z następujących plików źródłowych: makefile, Range.h, Range.cpp, RSDriver.h, RSDriver.cpp, RangeSearchDict.h, RangeSearchDict.cpp, LinkDef.h, rs.conf. Po ich ściągnięciu do tego samego katalogu, gdzie znajdują się pliki MyClass.h, MyClass.C, background.root, physics.root oraz tworzy_drzewa.C, należy je skompilować używając linuxowej komendy make. Zaleca się jednak podmienienie oryginalnego pliku makefile na kod podany poniżej (ze względu na mniej zawodne działanie):
# Makefile
CXX_COMP = gcc
LDOP = -g -shared
CPPFLAGS = -I$(ROOTSYS)/include
#
HDRS = Range.h RSDriver.h
SRCS = Range.cpp
OBJS = RSDriver.o Range.o
DICTOBJS = RangeSearchDict.o
DICTSRCS = RangeSearchDict.cpp
%.o: %.cpp
$(CXX_COMP) -c $(CPPFLAGS) -o $@ $<
all: libRangeSearch.so
libRangeSearch.so: $(OBJS) $(DICTOBJS)
$(CXX_COMP) -o $@ $(shell root-config --libs) $(LDOP) $^
RangeSearchDict.o: $(HDRS)
rootcint -f $(DICTSRCS) -c $^ LinkDef.h
$(CXX_COMP) -c $(CPPFLAGS) -o $(DICTOBJS) $(DICTSRCS)
clean:
rm *.o *.so *Dict.*
Range.o: Range.h Range.cpp
RSDriver.o: Range.h RSDriver.h
Dzięki temu zostanie utworzona biblioteka libRangeSearch.so niezbędna do analizy Range Searching. Dalszym krokiem jest modyfikacja pliku konfiguracyjnego rs.conf. Poszczególne zmienne powinny zawierać następujące wartości:
background=background.root #drzewo wzorcowe z NC
physics=physics.root #drzewo wzorcowe z CC
backgroundobj=background
physicsobj=physics
varlist=a,b,c #nazwy robocze trzech zmiennych
delta = 0.3,15,15000 #wartosci komórki V
W następnej kolejności należy stworzyć nowy plik wykorzystujący bibliotekę libRangeSearch.so i analizujący przypadki CC i NC metodą Range Searching. Nowy plik należy stworzyć w oparciu o metodę cięć (selekcja.C) z pewnymi modyfikacjami pierwowzoru [Dodatek B]:
#include <stdio.h>
void selekcja2() //selekcja2.C
{
//program analogiczny jak selekcja.C
//dla metody cięć
...
//dodawanie bibliotek
gSystem->Load(łibCandNtupleSR.so");
gSystem->Load(łibTruthHelperNtuple.so");
gSystem->Load(łibMCNtuple.so");
gSystem->Load(łibStandardNtuple.so");
gROOT->LoadMacro("MyClass.C");
gROOT->LoadMacro(łibRangeSearch.so"); //ładuje bibliotekę RangeSearching
MyClass t;
RSDriver *rs=new RSDriver(); //uruchamia Range Searching
...
Long64_t nentries = t.fChain->GetEntries();
...
//pętla po przypadkach
for (Int_t jentry = 0; jentry < nentries; jentry++)
{
int iaction = t.mc_iaction[0];
...
//definicje 3 zmiennych wejsciowych:
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;
float D=0; //prawdopodobienstwo - tego szukamy!
//wejsciowa tablica 3 zmiennych dla danego przypadku
float rsvars[3];
rsvars[0]=zmienna1;
rsvars[1]=zmienna2;
rsvars[2]=zmienna3;
//obliczenie D dla danego przypadku przy użyciu Range
//Searching dla tablicy zmiennych wejsciowych 'rsvars':
D=rs->rscalc(rsvars);
if (t.evthdr_ph_sigcor>0) //bez zerowej energii
{
//rysowanie rozkladu zmiennej D:
if (iaction==0) {histogram_NC -> Fill(D); ileD_NC++;}
if (iaction==1) {histogram_CC -> Fill(D); ileD_CC++;}
//selekcja CC/NC w oparciu o D:
if (D>=0.65) {czastka=1;} //to jest CC
else
if (D<=0.4) {czastka=0;} //to jest NC
else {czastka=-1;} //unknown
}
//na podstawie 'czastka' liczymy Pur i Eff
//analogicznie jak w metodzie cięć
...
} //koniec pętli po przypadkach
...
}
Należy pamiętać, iż przypadki wzorcowe zapisane są w plikach background.root oraz physics.root. Z kolei ścieżki do właściwych plików do przeanalizowania metodą RS w oparciu o wzorzec znajdują się w pliku MyClass.h (patrz: Dodatek B). Nazwy analizowanych plików:
//Daleki detektor, właczone pole B, wiazka niskoenergetyczna LE-10-185
f21011101_0000_L010185N_D00.sntp.cedar.root
f21011102_0000_L010185N_D00.sntp.cedar.root
f21011103_0000_L010185N_D00.sntp.cedar.root
f21011104_0000_L010185N_D00.sntp.cedar.root
f21011105_0000_L010185N_D00.sntp.cedar.root
f21011106_0000_L010185N_D00.sntp.cedar.root
f21011107_0000_L010185N_D00.sntp.cedar.root
f21011109_0000_L010185N_D00.sntp.cedar.root
f21011110_0000_L010185N_D00.sntp.cedar.root
f21011112_0000_L010185N_D00.sntp.cedar.root
f21011118_0000_L010185N_D00.sntp.cedar.root
f21011122_0000_L010185N_D00.sntp.cedar.root
f21011123_0000_L010185N_D00.sntp.cedar.root
f21011124_0000_L010185N_D00.sntp.cedar.root
- MODYFIKACJA ALGORYTMU RANGE SEARCHING
W rozdziale 5.9 przedstawiono sposób modyfikacji metody RS. W celu dokonania tej zmiany należy w pliku RSDriver.cpp podmienić istniejący kod funkcji rscalc na kod podany poniżej (autorstwa [11]):
#include <vector> \\na poczatku pliku
...
\\funkcja 'rscalc':
Float_t RSDriver::rscalc(Float_t *v, Float_t *ph_count, Float_t *bg_count)
{
Float_t l;
TVector *lb,*ub;
Rectangle range;
int i;
vector<Float_t>delta1(varnum);
UInt_t ipet=0;
do
{
lb=new TVector(varnum);
for(i=0;i<varnum; i++) //petla po zmiennych - u nas 3
{
(*lb)(i)=v[i]; //3 zmienne
}
ub=new TVector(*lb);
range.lower=lb;
range.upper=ub;
for(i=0;i< varnum;i++)
{
//delta - własciwa komorka V
//delta1 - dobudowane dodatkowe V'
delta1[i]=delta[i]+delta[i]*0.5*ipet;
(*lb)(i) -= delta1[i]*(1.-shift[i]);
(*ub)(i) += delta1[i]*shift[i];
}
pcount=phRangeSorter->search(&range);
bcount=bgRangeSorter->search(&range);
w2Ph=phRangeSorter->getSumW2();
w2Bg=bgRangeSorter->getSumW2();
delete lb;
delete ub;
++ipet;
}while(pcount<25. || bcount<25. );
//zwiekszaj V' dopoki sygnal i tlo = 25 przypadkow (kazde)
if(bg_count) *bg_count=bcount*bgScale;
if(ph_count) *ph_count=pcount*phScale;
//gdyby było za mało przypadków:
if(pcount<1e-10 && bcount<1e-10) return -1.;
if(bcount<1e-10) return 1.;
if(pcount<1e-10) return 0.;
l=pcount/bcount*phScale/bgScale; //parametr 'l'
//tutaj phScale=bgScale=1
return l/(1+l); //zwraca D dla danego przypadku
}
Dodatkowo w całym pliku RSDriver.cpp i plikach Range.cpp, Range.h oraz RSDriver.h należy zamienić typ float na Float_t we wszystkich miejscach, gdzie tylko występuje.
Zmodyfikowany algorytm selekcji wieloparametrycznej opartej o Range Searcing działa od pięciu do ośmiu razy dłużej niż algorytm podstawowy.
komentarze
Copyright © 2008-2010 EPrace oraz autorzy prac.