Kolekcje zderzeń i cięcia
W poprzedniej części analizowaliśmy zderzenia w funkcji parametru zderzenia oraz liczby cząstek. Niestety parametr zderzenia nie jest dostępny w pomiarach. Dlatego zwykle się go estymuje na podstawie liczby cząstek - krotności. Przy czym nie zawsze krotność oznacza całkowitą liczbę cząstek - czasem np. oznacza ona liczbę cząstek w pewnym obszarze kinematycznym.
W naszej wersji skupimy się na prostej estymacji parametru zderzenia i zobaczymy jak skuteczna jest ta metoda. Zaczniemy po pierwsze od oszacowania jak powinniśmy podzielić nasze zderzenia aby mieć centralność 0-10%, 10-40% i 40-100%. Wystarczy więc policzyć tzw. kwantyle dla histogramu z rozkładem liczby cząstek:
1 2 3 4 5 6 7 8 9 10 11 | void glauber_estim() { auto file = new Hal::AnaFile("/opt/temp/data.root"); auto h = file->GetHistogramPassed(Hal::ECutUpdate::kEvent, 0, 2); Double_t probSum[] = {0.6, 0.9}; Double_t q[3]; h->ResetStats(); h->GetQuantiles(2, q, probSum); for (int i = 0; i < 2; i++) cout << q[i] << endl; h->Draw();} |
Ponieważ musimy oddzielić 3 przedziały centralności wystarczy znaleźć 2 liczby. Pierwsza liczba odpowiada 60% najbardziej centralnych zderzeń i odpowiada centralności 40-100%, druga liczba odpowiada 90% najmniej centralnych zderzeń i pozwala wyznaczyć granicę między przedziałami 40-10 i 10-0%. Makro mi wyrzuciło wartości: 237 i 838.
Nie będziemy teraz generować danych w locie jak było to w poprzedniej lekcji ale użyjemy danych wygenerowanych przez makro glauber_gen z lekcji 1. Ponieważ w obu przypadkach używaliśmy tych samych ustawień będzie to działało dobrze. Nowe makro będzie wyglądać tak:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 | void glauber() { auto run = new Hal::AnalysisManager(); auto source = new Hal::RootSource("/opt/temp/glauber.root"); run->SetOutput("/opt/temp/dataA.root"); run->SetSource(source); HalOTF::Reader* reader = new HalOTF::Reader(); run->AddTask(reader); Hal::EventAna* ana = new Hal::EventAna(); ana->SetFormatOption(Hal::EventAna::EFormatOption::kReaderAccess); Hal::EventFieldMonitorXY prop(Hal::DataFieldID::Event::EMc::kB + Hal::DataFieldID::ImStep, Hal::DataFieldID::Event::EBasic::kTracksNo + Hal::DataFieldID::ImStep); prop.SetXaxis(100, 0, 20); prop.SetYaxis(200, 0, 1500); Hal::EventFieldMonitorX prop1(Hal::DataFieldID::Event::EMc::kB + Hal::DataFieldID::ImStep); Hal::EventFieldMonitorX prop2(Hal::DataFieldID::Event::EBasic::kTracksNo + Hal::DataFieldID::ImStep); prop1.SetXaxis(100, 0, 20); prop2.SetXaxis(200, 0, 1500); ana->AddCutMonitor(prop); ana->AddCutMonitor(prop1); ana->AddCutMonitor(prop2); Hal::EventTotalTrackNoCut nTrack1; nTrack1.SetMinMax(0, 237); ana->AddCut(nTrack1, "{0}+im"); Hal::EventTotalTrackNoCut nTrack2; nTrack2.SetMinMax(238, 838); ana->AddCut(nTrack2, "{1}+im"); Hal::EventTotalTrackNoCut nTrack3; nTrack3.SetMinMax(839, 2000); ana->AddCut(nTrack3, "{2}+im"); run->AddTask(ana); run->Init(); run->Run(0, 100000);} |
Początek makra wygląda podobnie jak w poprzednim wypadku. Source jest jednak inny, ponieważ nie analizujemy danych generowanych w locie ale z pliku ROOT-owego używamy klasy Hal::RootSource. Następnie używamy readera do konwersji danych do formatu HAL-owego. Analizę konfigurujemy w identyczny sposób jak poprzedni z wyjątkiem linii 23-33.
Te nowe linie to linie dodające tzw. cięcia. Cięcia to kryteria selekcji danych. W naszym wypadku chcemy podzielić zderzenia na 3 centralności, dlatego potrzebujemy 3 cięć. Generalnie to analizę można by przeprowadzić na 2 sposoby:
- stworzyć 3 taski typu EventAna i dla każdego ustawić cięciem inną centralność
- stworzyć 1 task typu EventAna ale dodać 3 cięcia tworząc 3 różne "kolekcje danych"
W tym makrze wybraliśmy tę drugą drogę. Jak widzimy do selekcji danych używamy klasy Hal::EventTotalTrackNoCut. Wybiera ona (tnie) zderzenia wg. kryterium całkowitej liczby cząstek. Wartościami SetMinMax ustawiamy zakres akceptowalny. Magia dzieje się w liniach 26,29 i 32. W linii 25 nasze cięcie jest konfigurowane tak, że wybiera zderzenia mające od 0 do 237 cząstek (włącznie), dodajemy je z flagą "{0}+im". Flaga "im" znana jest nam z poprzedniej lekcji i oznacza tyle że mamy "format zespolony" i patrzymy na część "urojoną" tj. tą z symulacji. Flaga "{0}" oznacza, że zderzenie zostanie dodane do kolekcji 0. Kolejne zderzenia dodawane są do kolekcji 1, 2. Czym są kolekcje?
Analizy w HAL-u działaja trochę jak sieci neuronowe, w zależności od sieci istnieje różna liczba warstw. W analizach typu EventAna istnieje 1 warstwa tych "neuronów" czy też "kolekcji" - warstwa "eventowa". Na początku analiza nie ma żadnej kolekcji, ale za każdym razem jak dodajemy cięcie, sprawdzane jest czy kolekcja do której dodajemy cięcie istnieje, jeśli nie - to kolekcja tworzona (niektóre analizy mogą mieć jednak maksymalna liczbę danego typu kolekcji - zwykle jak użytkownik "przegnie" to dostaje ostrzeżenie). W naszej analizie dodając cięcia do 3 nieistniejących kolekcji wymuszamy ich stworzenie. Nasza sieć ma więc takie "trzy neurony". Gdy AnalysisManager wywołuje naszą analizę/task to do każdego z naszych trzech "neuronów/kolekcji" płyną zderzenia, każdy neuron ma jednak swoje monitory. Dlatego w wyniku analizy otrzymamy nie 3 monitory ale 9, po 3 dla każdego z różnych cięć.
Jak więc widać manipulując liczbą kolekcji można wykonywać kilka analiz jednocześnie - choć jak pokażę później potęga kolekcji nie ujawnia się w analizach zderzeń, ale cząstek bądź par cząstek.
Ciekawostka: alternatywną opcją dodawania cięć jest opcja poniżej, zadziała ona tak samo bo HAL tworzy sobie kopię cięć przy dodawaniu.
1 2 3 4 5 6 7 | Hal::EventTotalTrackNoCut nTrack; nTrack.SetMinMax(0, 237); ana->AddCut(nTrack, "{0}+im"); nTrack.SetMinMax(238, 838); ana->AddCut(nTrack, "{1}+im"); nTrack.SetMinMax(839, 2000); ana->AddCut(nTrack, "{2}+im"); |
No dobra, a teraz pytanie jak wyglądają nasze rozkłady i jak się do nich dobrać? Użyjmy makra poniżej:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 | void glauber_show2() { auto file = new Hal::AnaFile("/opt/temp/dataA.root"); Hal::Std::MakeBeautiful(); gStyle->SetOptStat(0); TCanvas* c = new TCanvas(); c->Divide(2, 3); Int_t cent[3]; for (int i = 0; i < 3; i++) { cent[i] = file->GetPassedNo(Hal::ECutUpdate::kEvent, i); cout << cent[i] << endl; } for (int i = 0; i < 3; i++) { auto h_imp_good = file->GetHistogramPassed(Hal::ECutUpdate::kEvent, i, 1); auto h_tr_good = file->GetHistogramPassed(Hal::ECutUpdate::kEvent, i, 2); auto h_imp_bad = file->GetHistogramFailed(Hal::ECutUpdate::kEvent, i, 1); auto h_tr_bad = file->GetHistogramFailed(Hal::ECutUpdate::kEvent, i, 2); h_imp_good->SetFillColorAlpha(kGreen, 0.5); h_tr_good->SetFillColorAlpha(kGreen, 0.5); c->cd(1 + i * 2); h_imp_bad->Draw(); h_imp_good->Draw("SAME"); c->cd(2 + i * 2); h_tr_bad->Draw(); h_tr_good->Draw("SAME"); gPad->SetLogy(); }} |
Otrzymamy coś takiego:
Patrząc na terminal powinniśmy otrzymać liczbę zderzeń, które przeszły kolejne cięcia, jeśli nasze cięcia dobrze ustawiliśmy to pierwsza liczba powinna stanowić około 60% sumy wszystkich trzech wyświetlonych liczb, druga 30% a trzecia 10%.
Skupmy się teraz na interpretacji wyników. Po lewej widzimy rozkład liczby cząstek wg. której wybieraliśmy dane, zielony histogram to zderzenia, które przeszły cięcia, a czerwony to odrzucone cięcia. Widzimy, że nasze cięcia tak jak powinny wybierają określony zakres liczby cząstek. Po prawej widzimy analogiczny rozkład ale dla parametru zderzenia, ale tu widzimy, że nie mamy już tak ładnego "ostrego" rozdzielenia różnych wartości parametrów zderzenia. Podobnie w realnym eksperymencie, wybierając zderzenia wg mierzalnych kryteriów takich jak krotność nie jesteśmy w stanie idealnie odseparować zderzeń wg. określonego parametru zderzenia.