padsley / k600analyser

Code for the K600 analyser including plugin codes for silicon, clover and HAGAR data
1 stars 4 forks source link

Implementing developments in the master branch #172

Open padsley opened 5 years ago

padsley commented 5 years ago

Putting my project-manager hat on for a moment.

We have made fairly extensive improvements to e.g. the gamma-ray processing codes. These have been made in sub-branches of the K600 analyser.

I would like to implement these improvements in the master branch.

Please let me know which branches contain the improvements by creating an issue on the Github repo and tagging the issue in a reply to this issue. I shall probably implement the developments manually because there are a lot of changes made independently and that's probably the quickest way.

@KevinCWLi - I would like to take the opportunity to make sure that the new ray-tracing algorithm is either implemented instead of the old one or in addition (giving it as another option). Could you like me have a working version of the new raytracing algorithm so that it can be implemented?

KevinCWLi commented 4 years ago

Greetings K600

I've been assigned the project providing an example of the new developments of the k600analyser that have been used for the PhD thesis of KCW Li. The following repos have been updated to be consistent with the submitted results in the PhD thesis:

PR166_2010 (not used for the PhD work, but has been updated with the latest version) PR166_2011 PR194 PR195_2017 PR240 PR251

A good example to look at with relatively-little hardcoding is PR194: 12C(a,a')12C at 0 deg, Ebeam = 160 MeV. In this email, I've focused on the SilverBulletRaytrace development, however the TotalLineshapeCorrection() and X1mapping() functions could possibly also be of interest to future works. I've attempted to break this email up with comments to try and keep sections clear.

------------------------------------------------

The pre-processor variable shown below is how the SilverBulletRaytrace class is implemented. For each wireplane, a SilverBulletRaytrace object should be defined. Since this type of analysis is typically only applicable for the first horizontal wireplane (X1), only a single SilverBulletRaytrace object has been implemented in the example. The green code is where the code is made to work with functions that are not currently defined in the master branch. They can be implemented however people see fit and each green section has been (hopefully) described.

define _SILVERBULLETRAYTRACE

ifdef _SILVERBULLETRAYTRACE

SilverBulletRaytrace *silverBulletRaytrace_X1 = new SilverBulletRaytrace();

endif

------------------------------------------------

The branch is implemented as:

ifdef _SILVERBULLETRAYTRACE

gROOT->ProcessLine("#include \"SilverBulletRaytrace.h\"");
gROOT->ProcessLine(".L SilverBulletRaytrace.c+");
t1->Branch("SilverBulletRaytrace_X1","SilverBulletRaytrace", &silverBulletRaytrace_X1);

endif

------------------------------------------------

This takes the raw drift distances (i.e. without assumption of which wires are inverted) and executes the SBR algorithm of searching through all permutations (under some conditions).

Clearly, it is explicit that only sets of wires with >=3 and <=12 are considered as valid wireplane events... outside of this range, a focal-plane event within the SilverBulletRaytrace class is not created.

ifdef _SILVERBULLETRAYTRACE

if(X1hits_dt>=3 && X1hits_dt<=12)
{
    silverBulletRaytrace_X1->Raytrace(X1.dist, X1.wire, X1hits_dt);
}

endif

------------------------------------------------

For each "Raytrace" function call (see above), a total of N positions, returned by GetNAlternativeEvents(), is created.

Each of these positions can then undergo lineshape correction and the "SilverBulletRaytrace" object, silverBulletRaytrace_X1, then stores each corresponding lineshape-corrected position.

Or in other words, there is a vector of focal-plane positions determined by the Raytrace() function of SilverBulletRaytrace. This class also stores each corresponding line-shape corrected position.

The "TotalLineshapeCorrection" algorithm is not currently implemented in the master branch.

It is used for lineshape correction and could be implemented at the same time... but the point is that here is where the lineshape correction should be done.

ifdef _SILVERBULLETRAYTRACE

for(int i=0; i<silverBulletRaytrace_X1->GetNAlternativeEvents(); i++)
{
    double xPosition = silverBulletRaytrace_X1->GetXPosition(i);

    TotalLineshapeCorrection(correctionParameters, &xPosition);
    //TotalLineshapeCorrection(X1pos, t_Y1, thetaSCAT, &xPosition);

    silverBulletRaytrace_X1->SetXPosition_TLC(i, xPosition);
}

endif

------------------------------------------------

Similarly to the lineshape correction above, the lineshape-corrected positions can be mapped (i.e. offset, for simpler cases) and the SilverBulletRaytrace class also stores each mapped position.

ifdef _SILVERBULLETRAYTRACE

extern bool X1MappingDefined_SBR;

for(int i=0; i<silverBulletRaytrace_X1->GetNAlternativeEvents(); i++)
{
    if(X1MappingDefined_SBR)
    {
        double correctedXPosition = X1Mapping(silverBulletRaytrace_X1->GetXPosition_TLC(i));
        silverBulletRaytrace_X1->SetXPosition_TLC_mapped(i, correctedXPosition);
    }
}

endif

------------------------------------------------

Similarly for the preceding lineshape-correction and peak-position mapping sections, the SilverBulletRaytrace object also stores the excitation-energy associated with each lineshape-corrected and peak-position mapped focal-plane position.

The CalcEx() function is not implemented in the master branch but is provided in a modified FocalPlane.c file.

It obviously doesn't need to be used (unless you want to), any excitation-energy mapping can be implemented here (unless you want to do it after sorting, directly from the TTree).

ifdef _SILVERBULLETRAYTRACE

for(int i=0; i<silverBulletRaytrace_X1->GetNAlternativeEvents(); i++)
{
    if(momentumCalibrationRead_SBR)
    {
        double excitationEnergy = CalcEx(&m_SBR[0], &T_SBR[0], &E_SBR[0], &p_SBR[0], silverBulletRaytrace_X1->GetXPosition_TLC_mapped(i), polarScatteringAngles_SBR[2], momentumCalPars_SBR);

        silverBulletRaytrace_X1->SetExcitationEnergy(i, excitationEnergy);
    }
}

silverBulletRaytrace_X1->DetermineGoodSbrEvent();

endif

------------------------------------------------

At the end of the "main_event" routine, the following class function is called to clear the objects.

ifdef _SILVERBULLETRAYTRACE

silverBulletRaytrace_X1->ClearEvent();

endif

[http://cdn.sun.ac.za/100/ProductionFooter.jpg]http://www.sun.ac.za/english/about-us/strategic-documents

The integrity and confidentiality of this email are governed by these terms. Disclaimerhttp://www.sun.ac.za/emaildisclaimer Die integriteit en vertroulikheid van hierdie e-pos word deur die volgende bepalings bere?l. Vrywaringsklousulehttp://www.sun.ac.za/emaildisclaimer

KevinCWLi commented 4 years ago

It should be mentioned that the DetermineGoodSbrEvent function is to reorder the various positions (and associated parameters) determined by the SilverBulletRaytrace algorithm. These are ordered by the adjusted R-squared metric so that the position with the most linear permutation of wire positions is located at index 0.

================================================ Apologies, I forgot to mention that in the main_bor routine, the following must be executed once:

ifdef _SILVERBULLETRAYTRACE

GeneratePermutationLibrary(12);

endif

This generates the permutation library.

================================================ Hey Phil

Whilst implementing this is not particularly difficult, I appreciate that it could be easier if I did it instead. It should be relatively cut and paste into the master branch... if not we can take a look. Let me know how you wish for this to proceed. You are more than welcome to and substitute the current lineshape-correction function in... or just leave out the optional parts in green (but still calling the DetermineGoodSbrEvent function) and checking if it works.

Kind regards Kevin

On 02 Dec 2019, at 1:03 PM, Kevin Ching Wei Li kcwli@sun.ac.za<mailto:kcwli@sun.ac.za> wrote:

Greetings K600

I've been assigned the project providing an example of the new developments of the k600analyser that have been used for the PhD thesis of KCW Li. The following repos have been updated to be consistent with the submitted results in the PhD thesis:

PR166_2010 (not used for the PhD work, but has been updated with the latest version) PR166_2011 PR194 PR195_2017 PR240 PR251

A good example to look at with relatively-little hardcoding is PR194: 12C(a,a')12C at 0 deg, Ebeam = 160 MeV. In this email, I've focused on the SilverBulletRaytrace development, however the TotalLineshapeCorrection() and X1mapping() functions could possibly also be of interest to future works. I've attempted to break this email up with comments to try and keep sections clear.

------------------------------------------------

The pre-processor variable shown below is how the SilverBulletRaytrace class is implemented. For each wireplane, a SilverBulletRaytrace object should be defined. Since this type of analysis is typically only applicable for the first horizontal wireplane (X1), only a single SilverBulletRaytrace object has been implemented in the example. The green code is where the code is made to work with functions that are not currently defined in the master branch. They can be implemented however people see fit and each green section has been (hopefully) described.

define _SILVERBULLETRAYTRACE

ifdef _SILVERBULLETRAYTRACE

SilverBulletRaytrace *silverBulletRaytrace_X1 = new SilverBulletRaytrace();

endif

------------------------------------------------

The branch is implemented as:

ifdef _SILVERBULLETRAYTRACE

gROOT->ProcessLine("#include \"SilverBulletRaytrace.h\"");
gROOT->ProcessLine(".L SilverBulletRaytrace.c+");
t1->Branch("SilverBulletRaytrace_X1","SilverBulletRaytrace", &silverBulletRaytrace_X1);

endif

------------------------------------------------

This takes the raw drift distances (i.e. without assumption of which wires are inverted) and executes the SBR algorithm of searching through all permutations (under some conditions).

Clearly, it is explicit that only sets of wires with >=3 and <=12 are considered as valid wireplane events... outside of this range, a focal-plane event within the SilverBulletRaytrace class is not created.

ifdef _SILVERBULLETRAYTRACE

if(X1hits_dt>=3 && X1hits_dt<=12)
{
    silverBulletRaytrace_X1->Raytrace(X1.dist, X1.wire, X1hits_dt);
}

endif

------------------------------------------------

For each "Raytrace" function call (see above), a total of N positions, returned by GetNAlternativeEvents(), is created.

Each of these positions can then undergo lineshape correction and the "SilverBulletRaytrace" object, silverBulletRaytrace_X1, then stores each corresponding lineshape-corrected position.

Or in other words, there is a vector of focal-plane positions determined by the Raytrace() function of SilverBulletRaytrace. This class also stores each corresponding line-shape corrected position.

The "TotalLineshapeCorrection" algorithm is not currently implemented in the master branch.

It is used for lineshape correction and could be implemented at the same time... but the point is that here is where the lineshape correction should be done.

ifdef _SILVERBULLETRAYTRACE

for(int i=0; i<silverBulletRaytrace_X1->GetNAlternativeEvents(); i++)
{
    double xPosition = silverBulletRaytrace_X1->GetXPosition(i);

    TotalLineshapeCorrection(correctionParameters, &xPosition);
    //TotalLineshapeCorrection(X1pos, t_Y1, thetaSCAT, &xPosition);

    silverBulletRaytrace_X1->SetXPosition_TLC(i, xPosition);
}

endif

------------------------------------------------

Similarly to the lineshape correction above, the lineshape-corrected positions can be mapped (i.e. offset, for simpler cases) and the SilverBulletRaytrace class also stores each mapped position.

ifdef _SILVERBULLETRAYTRACE

extern bool X1MappingDefined_SBR;

for(int i=0; i<silverBulletRaytrace_X1->GetNAlternativeEvents(); i++)
{
    if(X1MappingDefined_SBR)
    {
        double correctedXPosition = X1Mapping(silverBulletRaytrace_X1->GetXPosition_TLC(i));
        silverBulletRaytrace_X1->SetXPosition_TLC_mapped(i, correctedXPosition);
    }
}

endif

------------------------------------------------

Similarly for the preceding lineshape-correction and peak-position mapping sections, the SilverBulletRaytrace object also stores the excitation-energy associated with each lineshape-corrected and peak-position mapped focal-plane position.

The CalcEx() function is not implemented in the master branch but is provided in a modified FocalPlane.c file.

It obviously doesn't need to be used (unless you want to), any excitation-energy mapping can be implemented here (unless you want to do it after sorting, directly from the TTree).

ifdef _SILVERBULLETRAYTRACE

for(int i=0; i<silverBulletRaytrace_X1->GetNAlternativeEvents(); i++)
{
    if(momentumCalibrationRead_SBR)
    {
        double excitationEnergy = CalcEx(&m_SBR[0], &T_SBR[0], &E_SBR[0], &p_SBR[0], silverBulletRaytrace_X1->GetXPosition_TLC_mapped(i), polarScatteringAngles_SBR[2], momentumCalPars_SBR);

        silverBulletRaytrace_X1->SetExcitationEnergy(i, excitationEnergy);
    }
}

silverBulletRaytrace_X1->DetermineGoodSbrEvent();

endif

------------------------------------------------

At the end of the "main_event" routine, the following class function is called to clear the objects.

ifdef _SILVERBULLETRAYTRACE

silverBulletRaytrace_X1->ClearEvent();

endif

[http://cdn.sun.ac.za/100/ProductionFooter.jpg]http://www.sun.ac.za/english/about-us/strategic-documents

The integrity and confidentiality of this email are governed by these terms. Disclaimerhttp://www.sun.ac.za/emaildisclaimer Die integriteit en vertroulikheid van hierdie e-pos word deur die volgende bepalings bere?l. Vrywaringsklousulehttp://www.sun.ac.za/emaildisclaimer

padsley commented 4 years ago

Ordinarily, I'd be delighted to shift the job off on someone else but there're other things which maybe should change in the master branch so I think that I'm going to have to take a look at it myself.

Dr Philip Adsley MA MSci (Cantab) Claude Leon Postdoctoral Fellow University of the Witwatersrand/iThemba LABS

On Mon, 2 Dec 2019 at 13:14, Kevin Ching Wei Li notifications@github.com wrote:

It should be mentioned that the DetermineGoodSbrEvent function is to reorder the various positions (and associated parameters) determined by the SilverBulletRaytrace algorithm. These are ordered by the adjusted R-squared metric so that the position with the most linear permutation of wire positions is located at index 0.

================================================ Apologies, I forgot to mention that in the main_bor routine, the following must be executed once:

ifdef _SILVERBULLETRAYTRACE

GeneratePermutationLibrary(12);

endif

This generates the permutation library.

================================================ Hey Phil

Whilst implementing this is not particularly difficult, I appreciate that it could be easier if I did it instead. It should be relatively cut and paste into the master branch... if not we can take a look. Let me know how you wish for this to proceed. You are more than welcome to and substitute the current lineshape-correction function in... or just leave out the optional parts in green (but still calling the DetermineGoodSbrEvent function) and checking if it works.

Kind regards Kevin

On 02 Dec 2019, at 1:03 PM, Kevin Ching Wei Li <kcwli@sun.ac.za<mailto: kcwli@sun.ac.za>> wrote:

Greetings K600

I've been assigned the project providing an example of the new developments of the k600analyser that have been used for the PhD thesis of KCW Li. The following repos have been updated to be consistent with the submitted results in the PhD thesis:

PR166_2010 (not used for the PhD work, but has been updated with the latest version) PR166_2011 PR194 PR195_2017 PR240 PR251

A good example to look at with relatively-little hardcoding is PR194: 12C(a,a')12C at 0 deg, Ebeam = 160 MeV. In this email, I've focused on the SilverBulletRaytrace development, however the TotalLineshapeCorrection() and X1mapping() functions could possibly also be of interest to future works. I've attempted to break this email up with comments to try and keep sections clear.

------------------------------------------------

The pre-processor variable shown below is how the SilverBulletRaytrace class is implemented. For each wireplane, a SilverBulletRaytrace object should be defined. Since this type of analysis is typically only applicable for the first horizontal wireplane (X1), only a single SilverBulletRaytrace object has been implemented in the example. The green code is where the code is made to work with functions that are not currently defined in the master branch. They can be implemented however people see fit and each green section has been (hopefully) described.

define _SILVERBULLETRAYTRACE

ifdef _SILVERBULLETRAYTRACE

SilverBulletRaytrace *silverBulletRaytrace_X1 = new SilverBulletRaytrace();

endif

------------------------------------------------

The branch is implemented as:

ifdef _SILVERBULLETRAYTRACE

gROOT->ProcessLine("#include \"SilverBulletRaytrace.h\""); gROOT->ProcessLine(".L SilverBulletRaytrace.c+"); t1->Branch("SilverBulletRaytrace_X1","SilverBulletRaytrace", &silverBulletRaytrace_X1);

endif

------------------------------------------------

This takes the raw drift distances (i.e. without assumption of which

wires are inverted) and executes the SBR algorithm of searching through all permutations (under some conditions).

Clearly, it is explicit that only sets of wires with >=3 and <=12 are

considered as valid wireplane events... outside of this range, a focal-plane event within the SilverBulletRaytrace class is not created.

ifdef _SILVERBULLETRAYTRACE

if(X1hits_dt>=3 && X1hits_dt<=12) { silverBulletRaytrace_X1->Raytrace(X1.dist, X1.wire, X1hits_dt); }

endif

------------------------------------------------

For each "Raytrace" function call (see above), a total of N positions,

returned by GetNAlternativeEvents(), is created.

Each of these positions can then undergo lineshape correction and the

"SilverBulletRaytrace" object, silverBulletRaytrace_X1, then stores each corresponding lineshape-corrected position.

Or in other words, there is a vector of focal-plane positions determined

by the Raytrace() function of SilverBulletRaytrace. This class also stores each corresponding line-shape corrected position.

The "TotalLineshapeCorrection" algorithm is not currently implemented in

the master branch.

It is used for lineshape correction and could be implemented at the same

time... but the point is that here is where the lineshape correction should be done.

ifdef _SILVERBULLETRAYTRACE

for(int i=0; iGetNAlternativeEvents(); i++) { double xPosition = silverBulletRaytrace_X1->GetXPosition(i);

TotalLineshapeCorrection(correctionParameters, &xPosition); //TotalLineshapeCorrection(X1pos, t_Y1, thetaSCAT, &xPosition);

silverBulletRaytrace_X1->SetXPosition_TLC(i, xPosition); }

endif

------------------------------------------------

Similarly to the lineshape correction above, the lineshape-corrected

positions can be mapped (i.e. offset, for simpler cases) and the SilverBulletRaytrace class also stores each mapped position.

ifdef _SILVERBULLETRAYTRACE

extern bool X1MappingDefined_SBR;

for(int i=0; iGetNAlternativeEvents(); i++) { if(X1MappingDefined_SBR) { double correctedXPosition = X1Mapping(silverBulletRaytrace_X1->GetXPosition_TLC(i)); silverBulletRaytrace_X1->SetXPosition_TLC_mapped(i, correctedXPosition); } }

endif

------------------------------------------------

Similarly for the preceding lineshape-correction and peak-position

mapping sections, the SilverBulletRaytrace object also stores the excitation-energy associated with each lineshape-corrected and peak-position mapped focal-plane position.

The CalcEx() function is not implemented in the master branch but is

provided in a modified FocalPlane.c file.

It obviously doesn't need to be used (unless you want to), any

excitation-energy mapping can be implemented here (unless you want to do it after sorting, directly from the TTree).

ifdef _SILVERBULLETRAYTRACE

for(int i=0; iGetNAlternativeEvents(); i++) { if(momentumCalibrationRead_SBR) { double excitationEnergy = CalcEx(&m_SBR[0], &T_SBR[0], &E_SBR[0], &p_SBR[0], silverBulletRaytrace_X1->GetXPosition_TLC_mapped(i), polarScatteringAngles_SBR[2], momentumCalPars_SBR);

silverBulletRaytrace_X1->SetExcitationEnergy(i, excitationEnergy); } }

silverBulletRaytrace_X1->DetermineGoodSbrEvent();

endif

------------------------------------------------

At the end of the "main_event" routine, the following class function is

called to clear the objects.

ifdef _SILVERBULLETRAYTRACE

silverBulletRaytrace_X1->ClearEvent();

endif

[http://cdn.sun.ac.za/100/ProductionFooter.jpg]< http://www.sun.ac.za/english/about-us/strategic-documents>

The integrity and confidentiality of this email are governed by these terms. Disclaimerhttp://www.sun.ac.za/emaildisclaimer Die integriteit en vertroulikheid van hierdie e-pos word deur die volgende bepalings bere?l. Vrywaringsklousulehttp://www.sun.ac.za/emaildisclaimer

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/padsley/k600analyser/issues/172?email_source=notifications&email_token=AB54VOHORNAXQDRYT4XHCY3QWTUZ3A5CNFSM4JLZVON2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEFTERAA#issuecomment-560351360, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB54VOCFH3IFR2EMWBO6LHDQWTUZ3ANCNFSM4JLZVONQ .

padsley commented 4 years ago

An update - the gamma-ray stuff from @HarshnaJ and @pellegri has now been included in the PR298 branch in the main repo. I think that I will merge that branch into the master at some point after DevOps2020 gets merged back into the master because that should then include all of the required modifications (I hope!).