Program w PYTHIA8 analizujący rozkłady produkowanych inkluzywnie K0 i Lambda

[AddahDraggon]

Wysłałem początkowe pliki, oczywiście będzie tam więcej histogramów; mam problem z wygenerowaniem wydarzeń z centralnej dyfrakcji

[AddahDraggon]

Wygląda na to że jednak jest w porządku - miałem błąd w jednej z flag do Pythii.

[ladamczy]

Wygląda rozsądnie. Robiąc histogramy prosze wyznaczyć prawdopodobieństwo że K0 zostanie zmierzone (filtered) w funkcji kinematyki K0 , najlepiej w 2D w przestrzeni. K0_pt - K0_eta. w akresie Pt - 0 do 2 Gev i eta -3 do 3

PS. Raczej mówi po polsku "przypadki" a nie "wydarzenia"

[ladamczy]

Dobrze by było aby równocześnie robił Pan histogramy dla Lambda, czyli dwa zestawy histogramów jeden dla K0 drugi dla Lambda (i anty-Lambda). To będzie szybsze niż puszczenie dwóch niezaleznych programów jeden dla K0 a drugi dla Lambda.

[ladamczy]

to powinno być -2212

https://github.com/ladamczy/STAR-Analysis/blob/a18d866a51974e8ce5441a28df14fef247ba3865/Inclusive%20Analysis/main999.cc#L105

[AddahDraggon]

Poprawiłem; nigdzie tego konkretnie nie wykorzystywałem, także następny comit będzie już z dodanymi kolejnymi wykresami.

[ladamczy]

Prosze generować centralną dyfrakcje dla 510 GeV

Beams:eCM = 510 SoftQCD:centralDiffractive = on SigmaTotal:zeroAXB = off

Bede nas interesowały energie protonów dyfrakcyjnych a w zasadzie tak zwane xi czyli względna strata energii protonu. Przypadek CD jest scharakteryzowany dwiema wartościami: xi_east, xi_west

Najpierw trzeba zidentyfikować protony dyfrakcyjny

std::abs(event[i].id()) == 2212 && pythia.event[event[i].mother1()].id() == 2212

czyli proton który pochodzi od protonu. Takie protony powinny być zawsze dwa (prosze to sprawdzać) Jeden ma dodatniw pZ drugi ujemne (east, west albi naodwrót)

Jak juz Pan ma protony to liczy Pan xi kazdego z nich:

xi = (255.-event[i].e())/255. (gdzie i to i protonu)

Na początek prosze zrobić rozkłady

2D xi_east,xi_west xi_east*xi_west ln(xi_east/xi_west)

Teraz musimy wybrać przypadki które maja przynajmniej dwie naładowane cząstki w stanie końcowym w zakresie TPC (pT>0.2, |eta|<0.9) metody event[i].isCharged() && event[i].isFinal()

Robimy histogramy

2D xi_east,xi_west xi_east*xi_west ln(xi_east/xi_west)

pod warunkiem ze sa dwie cząstki w stanie końcowym

Teraz dodajemy warunek że nic nie ma w BBC Large

Czyli z przypadków które mają cząstki naładowane w TPC odrzucamy przypadki które mają cząstki naładowane w BBC Large

Tutaj prosze dać wiekszy próg na pT>0.5 zakres rapidity BBC Large 2.1 < |η| < 3.3

Znowu robimy histogramy

2D xi_east,xi_west xi_east*xi_west ln(xi_east/xi_west)

pod warunkiem że przynajmniej dwie cząstki sa TPC i żadnej w BBC large.

Dopiero teraz szukamy K0 które rozpada sie na dwa piony w zakresie akcpetnacji TPC

Znowu robimy histogramy

2D xi_east,xi_west xi_east*xi_west ln(xi_east/xi_west)

z dodatkowym rządaniem że widzimy K0.

Tutaj prosze jeszcze dodać histogram pokazujacy ile jest innych czastek naładowanych w zakresie TPC oprócz dwóch pionów z rozpadu K0.

Znowu robimy histogramy

2D xi_east,xi_west xi_east*xi_west ln(xi_east/xi_west)

ale z dodatkowym wymogiem że tych cząstek w TPC mamy przynajmniej dwie sztuki.

Myślę, że to na początek wystarczy.

[AddahDraggon]

Jak rozumiem, mam również używać flag

  pythia.readString("310:onMode=0");
  pythia.readString("310:OnIfMatch=211 -211");
  pythia.readString("-310:onMode=0");
  pythia.readString("-310:OnIfMatch=-211 211");

?

[ladamczy]

Tak, na obecnym etapie interesuje nas najlepsze możłiwe przefiltrowania.

[ladamczy]

Widzę że kod został dodany do repozytorium ale nie widze plotów/histogramów. Wolałbym odwrotną kolejność. Najpierw ploty/histogramy potem aktualizacja kodu.

[AddahDraggon]

Rozumiem - ploty mam dodawać w tym issue, czy jakoś inaczej?

[ladamczy]

w tym issue

[AddahDraggon]

Zrobiłem histogramy związane z xi, muszę jeszcze dopisać kod do ilości cząstek niepochodzących od K0, ale póki co wysyłam co mam. result.pdf

[ladamczy]

Czy może je Pan zrobić w formie stosunku histogramu po danym kroku do histogramu z kroku poprzedniego.

Czyli TPC/ALL Potem BBCLarge/TPC Potem K0/BBCL

Dziwi mnie efekt BBCL . Czy tu napewno stosuje Pan veto w BBCL?

[AddahDraggon]

Tak, stosuję - cząstek jest ok. 2x mniej po veto dla BBCL niż przed. Robienie ilorazu histogramów ma się odnosić do wszystkich oprócz tego z cząstkami nie od K0?

[ladamczy]

iloraz dla rozkładów xi1*xi2. oraz ln(xi_1/xi_2)

[AddahDraggon]

Zrobione, obecna symulacja jest dla 100 000 przypadków, także pod koniec jest dosyć mało przypadków - mogę ją swobodnie zrobić dla znacznie większej liczby. Przesyłam również plik z dodanym wykresem liczby cząstek wykrytych niepochodzących od K0. 100000_510GeV_central_diffraction_ratios.pdf 100000_510GeV_central_diffraction.pdf

[ladamczy]

Pewny komentarz pózniej ale juz widac ze zamiast rozkładów xi1xi2 prosze robić log(xi1xi2)

[ladamczy]

Kilka uwag:

1. Musimy wybrać odpowiedni zakres kinematyczny naszego pomiaru. Proces centralnej produkcji p+p->p+X+p opisuje się masą stanu X (która jest związana z xi1xi2. M_X^2=xi1xi2*s) oraz rapidity stanu X (które jest związane z xi1/xi2. y_X=ln(xi1/xi2).

2. Aby stan X był chociać częściowo był widziany przez TPC (w postaci przynajmniej 2 śladów) to stan X musi mieć względnie dużą masę i być produkowanym w miare centralnie. Z tych plotów możemy wstępnie przyjąć że naszym zakresem kinemtycznym będzie : xi1xi2>0.005 (aczkolwiek dokładną wartość poznamy jak zrobi Pan wykres log(xi1xi2) oraz |ln(xi1/xi2)|<4

3. Prosze zrobić ploty stosunku log(xi1xi2) po i przed cieciem na TPC z waruniem |ln(xi1/xi2)|<4 stosunku ln(xi1/xi2) po i przed cieciem na TPC z warunkiem xi1xi2>0.005 stosunku 2D log(xi1) vs. log(xi2) po i przed cieciem na TPC z waruniem |ln(xi1/xi2)|<4 Ten ostatni plot jest poto aby sprawdzić czy zamiast ciecia xi1*xi2>0.005 nie lepiej by było mieć ciecie xi1>0.07 i xi2>0.07. Wtedy moglibyśmy zrezygnować w danych z ciecia antyelastycznego bo ono odrzuca trochę dobrych przypadków a protony elastyczne mają xi1 approx 0.0 i ciecia xi1>0.07 i xi2>0.07 automatycznie wycinają elastykę i np tło przypadkowe typu beam-halo proton + pojedyńcza dyfrakcja.

4. Efekt weta w BBCL jest dość dramatyczny. Być może trzeba będzie poprzestać na analizie tylko trygera bz weta ale zanim o tym zdecydujemy to prosze sprawdzić jeszcze efekt weta BBCL ale z wrunkiem |ln(xi1/xi2)|<4 (dla plotu log(xi1xi2)) oraz warunku xi1xi2>0.005 dla plotu ln(xi1/xi2)

5. Prosze zrobić plot (ratio) sumarycznego efektu TPC+BBCL czyli stosunku log(xi1xi2) po i przed cieciem na TPC&BBCL z waruniem |ln(xi1/xi2)|<4 stosunku ln(xi1/xi2) po i przed cieciem na TPC&BBCL z warunkiem xi1xi2>0.005

6. ostatni plot to plot pokazujący rozkład ilości śladów w TPC z warunkiem |ln(xi1/xi2)|<4 & xi1*xi2>0.005 i warunkim na to ze widzimy w TPC dwa ślady z K0_S (podobny plot dla Lambdy)

7. Dobrze by było aby Pan cos pokazał z tych studiów na zebraniu z Czechami w przyszłym tygodniu.

[AddahDraggon]

Czy mam powyższe wykresy dołączyć do wykresów już zrobionych dotychczas? Jak rozumiem, powyższe wykresy dla celów przentacji mam złożyć do jakiegoś rodzaju krótkiej prezentacji w PDFie? Będę musiał chwilowo poczekać z pisaniem kodu do selekcji inkluzywnej żeby zdążyć z tymi wykresami, także wyniki z selekcji inkluzywnej będę mógł zaprzentować dopiero na następnym spotkaniu.

[ladamczy]

Może Pan dołączyć a może to być dodatkowy plik, jak Pan woli. Wyniki tej analizy PYTHII powinny być w formie krótkiej prezentacji w PDF. To jest chyba jedyne co Pan może pokazać na spotkaniu z Czechami.

[AddahDraggon]

Obecne wyniki wyglądają tak (obrazek w PDFie jest niezbyt dobry, dlatego go też dałem osobno): 100000_510GeV_Central_diffraction_v2.pdf W prezentacji o czym powinienem wspomnieć oprócz pokazania jakie cięcia zakładamy jako TPC i BBCL?

[AddahDraggon]

Obrazki oczywiście poprawię gdzie trzeba przed wrzuceniem ich do prezentacji

[ladamczy]

nie wiem czy do 12 uda sie Panu sklecić jakąś prezentacje, ale niech Pan próbuje zdecydujemy czy to ma sens zaraz przed zebraniem.

Generalna uwaga co do plotów: Tylko xi1/xi2 obkładamy logarytmem naturalnym (wtedy ta wielkość odpowiada rapidity stanu centralnego) pozostałe wielkości czyli xi1, xi2, xi1*xi2 obkładamy logarytmem przy podstawie 10.

Plot xi1xi2 pokazuje ze zakres xi1xi2>0.005 jest zbyt restrykcyjny , wydaje się że możemy przyjąć zakres xi1xi2>0.00005 Z drugiej strony rozkład 2D pokazuje że zamiast xi1xi2>0.00005 mozna by żadać xi1>0.007 and xi1>0.007 Wiec prosze wszędzie zamienić xi1*xi2>0.005 na xi1>0.007 and xi1>0.007

Plot ln(xi1/xi2) niech Pan tutaj nie stosuje ciecia |ln(xi1/xi2)|<4 . Chcemy zobaczyć jak ten stosunek szybko spada dla większych wartości.

Weto w BBCL pokazuje że powinniśmy mieć górne cieciei xi1*xi2<0.006 czyli można by sie ppokuśić na zakres fiducial

0.007 < xi1 < 0.08 and 0.007 < xi1 < 0.08

Coś jest nie tak z tymi wydajnościm weta BBCL i połączenia TPC & BBCL w funkcji ln(xi1/xi2) , //// dla ln(xi1/xi2) =0 powinno być maksimum a nie minimum.

[AddahDraggon]

Zmieniłem wartości i cięcia, dla ln(xi1/xi2) =0 faktycznie widać teraz maksimum. Zrobiłęm symulację dla próbki 1 000 000 symulacji, nie mam pojęcia dlaczego cięcie dla wykresów ln(xi_E/xi_W) jest takie równe, w kodzie nie ma tego typu cięć. 1000000_510GeV_Central_diffraction_v2.pdf

[AddahDraggon]

Obawiam się, że na to zebranie nie będę w stanie zrobić tej prezentacji.

[ladamczy]

Wygląda dobrze, ostre obciecie na ln(xi1/xi2) jest jasne ono wynika z ciecia 0.007 < xi1 < 0.08

ln(xi_max/xi_min) = ln (0.08/0.007) = 2.44

Musimy jeszcze popatrzeć na dane ale z analizy PYTHIA wynika że nasz obszar pomiarowy mniej więcej jest dany przez 0.007 < xi1(2) < 0.1
Górne ciecie jest łatwe (to znaczy łatwo zmierzyć xi=0.1) natomiast z dolnym cieiciem bedzie problem bo rozdzielczość pomiaru xi jest pewnie rzedu 0.005 zatem tam sie wszystko zlewa w jeden pik wokół zera. Tam tez będzie duże tło od beam halo. Trzeba popatrzeć do danych ale prawdopodobnie będziemy musieli przyciąć na trzy sigma czyli 0.01 ?

Górne ciecie xi=0.1 wynika z weta na BBCL . Dla trygera bez BBCL zakres możemy poszerzyć do xi<0.2

[AddahDraggon]

Poprawiłem wykresy wydajności, teraz korzystają z obiektów TEfficiency, oraz wykres ilości wykrytych K0 w zależności od całkowitej liczby cząstek. Chciałbym się też zapytać, jaki jest cel tych wykresów - przy warunku widoczności K0 to by nam dało spadek wydajności wykrywania K0 przy danych cięciach, a tak to chyba nie rozumiem co my tym osiągamy.

1000000_510GeV_Central_diffraction_v2.pdf

[ladamczy]

Jednym zdaniem:

Chcemy przeprowadzić pomiar który wnosi cos nowego do obecnego stanu wiedzy.

Zatem nie interesuje nas sam mezon K0 ktorego własnosci sa bardzo dobrze znane na podstawie pomiarów wielu wczesniejszych eksperymentów. Raczej insteresuje nas prawdopodobieństwo(warunkowe) pojawienia sie mezonu K0 pod warunkiem, ze jego źródłem jest jakis konkretny proces fizyczny. W naszym przypadku tym procesem jest centralna dyfrakcja p+p->p+X+p. Centralna dyfrakcja (z punktu widzenia stanu X) jest opisana dwiema wilekościami masa i rapidity tego stanu M_X i y_X. W eksperymencie nie mierzymy wszytskich czastek w stanie X ale korzystamy z zasady zachowania cztero-pedu która wymusza korelacje między stanem X a rozproszonymi protonami. Korelacja to prowadzi do zalezności

M_X^2 = xi_Exi_Ws y_X = ln(xi_E/xi_W)

Zatem punktem wyjścia (zanim zaczniemy cokolwiek mówić o K0) jest określenie z jakiego procesu pochodzi badane K0. Jak Pan zobaczy prawdopodobieństwo produkcji K0 zalezy od M_X i y_X zatem nasz pomiar musi prowadzić do okreslenia prawdopodobieństwa warunkowego P(K0|M_X,y_X) lub P(K0|xi_W,xi_E) w innym wypadku jest on nic nieznaczacy.

Tutaj wchodi drugi problem. Pomiar musimy poprawić na efekty detektorowe (geometryczne i wydajnosciowe). Nie mozemy opublikowac prawdopodobieństwa produkcji P(K0|xi_W,xi_E) pod warunkiem że STAR to mierzy. Bo to np. uniemożliwia porownanie naszych wyników z teoria. Takie porównanie będzie mozliwe tylko wtedy gdy teorie np. przewidywania PYTHII przepuscimy przez symulacje STAR ale ludzie z poza STAR nie moga tego zrobić. Dlatego musimy to poprawić na efekty detektorowe. Poprawki nie moga być zbyt duze. w praktyce przyjmuje sie ze pomiar powinien byc ograniczony do zakresu (tak zwany zkres fidutial ) gdzie wydjaność jest większa 10-20 % . I to jest właśnie to co teraz Pan robi. Mam nadzieje ze to narazie wystarczy

[AddahDraggon]

Tak, rozumiem. Jakie wykresy powinny się znaleźć w prezentacji do najbliższego spotkania oprócz wykresów wydajności cięć TPC i BBCL i płynących z tego ograniczeń?

[ladamczy]

Dobrze by było jeszcze zrobić:

1. rozkłady pT i eta K0 (w zakresie fiducial) (może również L(Lbar) )
2. Prawdopodobieństwo obserwacji K0(L) w funkcji M_x (lub log(xi_w xi_E))
3. Cokolwiek z danych (nawet bez K0) czyli lozkłady M_X (lub log(xi_w xi_E)) log(xi_e/xi_w) osobno log(xi) i/lub 2D log(xi_E) vs. log(xi_W)

Dobrze by było aby Pan przygotował szkic prezentacji na poniedzialek.

[AddahDraggon]

Czy w prezentacji mam umieścić również wykresy bez żadnych cięć na eta i napisać że wybieramy te obecnie używane?

Analogicznie, czy rozkłady eta i pT z K0 powinny być z uwzględnieniem cięć na TPC i BBCL czy bez? (czy oba, dla porównania)?

[ladamczy]

Rozkłady eta i Pt K0 po cieciach na TPC i xi.
Pewnei trzeba je mieć w dwóch wersjach z i bez BBCL (ba mamy dane z i bez BBCL).

Nie rozumiem o jakie wykresy bez ciec na eta Panu chodzi. Najlepiej zrobic jhakiśc szkic prezentacji i ja przedyskutować, Powinna mięc 10-15 minut więc trzeba wybrać tylko te najwazniejsze ploty. Generalnie chodzi o to aby pokaząc jak ciecia na TPC (i ewentualnie BBCL) ograniczaja nasz zakres fiducial

[AddahDraggon]

Mam szkielet prezentacji; przepraszam że tak późno. https://www.overleaf.com/9641628487nxkzntybthgk#2899c7 Chciałbym się zapytać, co powinienem poprawić w tych wykresach, żeby resztę już zrobić na gotowo w podobnym stylu.

[ladamczy]

wykresy wyglądają dobrze. Nie wiem czy to ciecie ln(xie/xiw)<4 powinno być eksponowane bo po cieciu na 0.007<xie(w)<0.08 jest ono nieefektywne ( nic nie wycina). Może prosze napisać alternative cuts ln(xie/xiw) ; log(xie*xiw)

[AddahDraggon]

Podczas robienia kodu do obrazków zauważyłem że jedyne cięcie jakie stosuję to to |ln(xi_E / xiW)|<4, a cięcia 0.007<xi{E,W}<0.08 stosuję tylko do wykresów eta-pT. To chyba oznacza że będę musiał przerobić slajd z wyszczególnionymi cięciami, oraz wszystkie tytuły i ten z podsumowaniem.

[ladamczy]

To niech Pan zmieni tytuły tak aby odpowiadały temu co Pan zrobił.

Slajd z cieiciami i podsumownaime trzeba tylko troche zmienić . Zrobię to teraz.

[ladamczy]

Zmieniłem sljad 2. Podsumowanie może zostać.
Strone 12 niech Pan da jako 9 (przed rozkladami K0)

Na stronie 12 prosze usunąć to kółko i spróbowac narysować obszar

\item $2\log(\xi_{min})< \log(\xi_E\cdot\xiW)<2\log(\xi{max})$ \item $|y_X|=|\ln(\xi_E/\xiW)|<\ln(\xi{max}/\xi_{min})$

czyli

$2\log(0.005)< \log(\xi_E) + \log(\xi_W)<2\log(0.08)$ $|\ln(\xi_E) - \ln(\xi_W)|<\ln(0.08/0.005) = 2.77$

to pewnie będzie jakiś romb a nie kółko

[ladamczy]

zmieniłem również tytuł prezentacji (na jedno K0)

[ladamczy]

Troche mnie niepokoi fakt ze wydajność TPC jest zawsze <80% . spodziewałbym się 100% dla dużych mass.

Czy to nie jest tak że warunek ln(xi_E/xi_W)<4 ma Pan tylko na licznik w tych stosunkach ? Jesli tak jest to prosze to poprawić i dać warunek ln(xi_E/xi_W)<4 również na mianownik.

[AddahDraggon]

Sprawdziłem w programie, cięcie ln(xi_E/xi_W)<4 jest zarówno w mianownikach, jak i w licznikach. Czy przy z grubsza równomiernym rozkładzie cząstek możliwe jest że po prostu TPC zajmuje 80% pełnego kąta bryłowego i niewykryte cząstki uciekają "denkami"?

[ladamczy]

Wracamy do tego issue,

Odpowiadając na ostatnie pytanie: Czastki uciekaja ale to co wyżej oznacza, że nie widzimy nic w obszarze TPC to się może zdażyc dlamałych mass, ale im większa masa stanu X tym to sie powinno zdażać rzadziej. .

Tutaj musimy jescze troche dodać wybików z poziomu true:

1. Jaka jest krotność K0 ?
2. Jaki jest stosunke N_K0/N_events (gdzie N_events jest ilościa przypadków "mierzonych")
3. Jak często oprócz K0 mamy jeszcze przynajmniej dwa ślady w zakresie TPC (wierzchołek pierwotny)
4. Jak liczby z punktów 1,2,3 zależą od masy stanu X
5. jaki liczby z punktów 1,2,3,4 wyglądaja dla Lambda/Lambdabar

[AddahDraggon]

Znalazłem błąd w symulacji, naprawiłem. Plik z histogramami:

10000000_510GeV_Central_diffraction_K0andLambda.zip

Jest tu 10mln przypadków. Zakładam brak cząstek w BBCL. Są histogramy do podpunktów 1,3,4 i 5, podpunkt 2 to liczba przypadków w histogramach 2D krotności/10mln

[AddahDraggon]

Poprawka co do ilości przypadków i stosunku N_part/n_events (n_events to przypadki które przeszły cięcia): Na 10 000 000 przypadków miałem: N_events: 13950 N_K0: 642 N_Lambda: 57 N_K0/N_events: 0.0460 N_Lmbda0/N_events: 0.0041

[ladamczy]

Nie rozumiem histogramu Additional tracks . Po pierwsze w opisie jest at Least two additional ... . Czy to oznacza że Pana selekcja to faktycznie jest K0(Lamba) + przynajmniej dwa inne ślady ? Jesli tak to prosze usunąc te przynajmniej dwa. Dlatego że właśnie chcemy sprawdzić ile jest śladów oprócz K0(Lambda) bez warunku że sa conajmniej dwa.

Po drugie dlaczego jest kilka wejść w tym histogramie dla wartości 0 i 1 a potem olbrzymi skok do 2 i więcej. Co to oznacza ?

[AddahDraggon]

Additional tracks to histogram liczony według positive_tracks+negative_tracks-2*number_of_Lambda(or K0), zatem faktycznie jest źle nazwany.

Zatem ten duży skok prawdopodobnie oznacza że bardzo często oprócz produktów rozpadu K0/Lambda TPC łapie również dwie inne cząstki.

[ladamczy]

No to juz lepiej to wygląda. Właśnie chciałem sie upewnić czy warto patrzeć na przypadki które maja tylko K0. Z tego wynika że znacznie więcej mamy przypadków z dodatkowych sladem. Oczywiście to jest poziom true przy wydajności rzedu 60% to juz nie bedzie tak wyglądało ale mimo wszystko pewnie wiecej niz 50% przypadków bedzie miało dodatkowe slady.

Może Pan to "zasymulować" generujac liczbe losowa z rokładu jednostajnego [0.1] dla kazdego śladu (oprócz tych dwóch z K0) i akceptowac slad jeśłi liczba losowa jest <0.6 . W PYTHII jest generator liczb losowych.

Czy Pan wymusza aby K0 rozpadało sie tylko na pi+pi- ?

[AddahDraggon]

Tak, wymuszam rozpad na pi+pi-