TRestAnalysisPlot: using a std::map observable for cuts results in strange behaviour

RESTFramework
1

REST version : v2.3.4

Hello!

Me and @CMargalejo were trying to apply veto cuts using TRestAnalysisPlot. In the .rml we define in the TRestAnalysisPlot-element

<globalCut variable="veto_MaxPeakAmplitude.second" condition="<500" value="ON" />

where veto_MaxPeakAmplitude is a std::map<int,Double_t>, containing the channel IDs as keys and the amplitudes as values.

However, this cut increases the number of events in the final plot, compared to turning the cut off by setting value="OFF". How can this be?
On the other hand. using a lower threshold in the condition, for example condition="<100" reduces the number of events compared to no cut.

This is how the draw arguments looks like (including many other cuts).

 AnalysisTree->Draw("hitsAna_energy>>HitsEnergy_110(48 , 0 , 19000)", "tckAna_nTracks_X==1 && tckAna_nTracks_Y==1 && hitsAna_nHitsX>0 && hitsAna_nHitsY>0 && hitsAna_xMean*hitsAna_xMean+hitsAna_yMean*hitsAna_yMean<100 && hitsAna_energy<19000 && hitsAna_balanceXYenergy>-0.2 && hitsAna_balanceXYenergy<0.4 && hitsAna_xy2Sigma>0.05 && hitsAna_xy2Sigma<0.7 && hitsAna_xySigmaBalance>-0.3 && hitsAna_xySigmaBalance<0.4 && hitsAna_z2Sigma<4 && hitsAna_xySkew>-1 && hitsAna_xySkew<1 && hitsAna_zSkew>-2 && hitsAna_zSkew<1.8 && veto_MaxPeakAmplitude.second<500", "colz", 9223372036854775807, 0)

Another problem occurs, when I try to take only a single veto panel into account by selecting e.g. channel 4645:

<globalCut variable="veto_MaxPeakAmplitude[4645]" condition="<500" value="ON" />

In this case, not a single event remains after the cut, even if I set a condition that is fulfilled by every veto signal.

Maybe TTree->Draw() cannot handle maps for the cuts? I read this in the root forum, but I am not sure if it is (still) true… If it is true, we would have to either expand TRestAnalysisPlot or change how the observables are stored in TRestRawVetoAnalysisProcess.
I also read somewhere that only a maximum of 16 cuts can be accepted, but I reproduced the problems with only the veto cut active.

Any ideas how to proceed?

2

Hi!

Not sure if that is a valid construction. I remember trying the following post structure, and it worked.

http://rest-forum.unizar.es/t/trestanalysisplot-capability-to-plot-map-variables/322

could you try a similar thing to see if that works?

Trying to imagine a posibility, perhaps when this cut is applied it will create an entry for each veto signal evaluated that fulfills <500?

I didnt know that, where did you read it? We should then place a warning and limit the number of cuts in TRestAnalysisPlot.

My advice is to create a simple TTree with few entries containing hitsAna_energy and the vetoSignal map. Then, since it is a valid Draw construction place the question at the https://root-forum.cern.ch sharing the TTree and asking if it is a valid construction.

3

If at the root-forum tell you this is not possible to do, we can try to find out how to integrate that at TRestAnalysisPlot.

There are other related features here:

4

It is probably not applicable here, because it is only for graphical cuts in histograms: Chapter: Histograms

5

I’ve tried using a cutString as you suggest, but in this case the cut does nothing. I always get the same number of events… I’ve tried both:

<cutString string="veto_MaxPeakAmplitude.first<500" value="ON" />
and
<cutString string="veto_MaxPeakAmplitude[4645]<500" value="ON" />

But none of them works. I’ve also tried with different conditions, e.g., >1000, < 200… and I always get the same number of events.

6

But, the keyword first corresponds to the ID. What if you try:

<histo name="SignalId_145">
<variable name="hitsAna_energy" range="(0,19000)" nbins="48" />
<cutString string="veto_MaxPeakAmplitude.second<500" />
</histo>

But I am afraid it will add an entry for each veto signal that fulfils the condition.

Then, as I proposed earlier, the best would be to prepare a small ROOT file with a pure TTree that contains the std::map and the hitsAna_energy. Test the Draw method, and ask at the ROOT forum the best way to achieve what you want.

7

I think you are right. After the X-ray cuts I have 75 events, and I expect the veto cut to reduce this by a factor of 2. However, after the veto cut we get 328 entries. We have 9 vetoes, so 328/9 ~ 36 entries, which is more or less what I expected to have.
If this is the case, the histogram is kind of rescaled by the number of vetoes.

8

Of course, a way to achieve what you want is to produce an additional observable at the analysisTree, inside the TRestRawVetoAnalysisProcess. For example, something like vetoActivated_500. If any of the veto signals is below <500, then you set this observable to 1, if not you set the observable to 0.

Later, in TRestAnalyisPlot you may use vetoActivated==1 to select those events.

The value of 500 could be left as a free parameter being read at the InitProcess. We do that for example at TRestTrackBlobAnalysisProcess, where we read the value of the radius from the observable name.

 void TRestTrackBlobAnalysisProcess::InitProcess() {
   52     std::vector<string> fObservables;
   53     fObservables = TRestEventProcess::ReadObservables();
   54 
   55     for (unsigned int i = 0; i < fObservables.size(); i++) {
   56         if (fObservables[i].find("Q1_R") != string::npos) fQ1_Observables.push_back(fObservables[i]);
   57         if (fObservables[i].find("Q2_R") != string::npos) fQ2_Observables.push_back(fObservables[i]);
   58 
   59         if (fObservables[i].find("Q1_X_R") != string::npos) fQ1_X_Observables.push_back(fObservables[i]);
   60         if (fObservables[i].find("Q2_X_R") != string::npos) fQ2_X_Observables.push_back(fObservables[i]);
   61 
   62         if (fObservables[i].find("Q1_Y_R") != string::npos) fQ1_Y_Observables.push_back(fObservables[i]);
   63         if (fObservables[i].find("Q2_Y_R") != string::npos) fQ2_Y_Observables.push_back(fObservables[i]);
   64 
   65         if (fObservables[i].find("Qhigh_R") != string::npos) fQhigh_Observables.push_back(fObservables[i]);
   66         if (fObservables[i].find("Qlow_R") != string::npos) fQlow_Observables.push_back(fObservables[i]);
   67         if (fObservables[i].find("Qbalance_R") != string::npos)
   68             fQbalance_Observables.push_back(fObservables[i]);
   69         if (fObservables[i].find("Qratio_R") != string::npos) fQratio_Observables.push_back(fObservables[i]);
   70 
   71         if (fObservables[i].find("Qhigh_X_R") != string::npos)
   72             fQhigh_X_Observables.push_back(fObservables[i]);
   73         if (fObservables[i].find("Qlow_X_R") != string::npos) fQlow_X_Observables.push_back(fObservables[i]);
   74         if (fObservables[i].find("Qbalance_X_R") != string::npos)
   75             fQbalance_X_Observables.push_back(fObservables[i]);
   76         if (fObservables[i].find("Qratio_X_R") != string::npos)
   77             fQratio_X_Observables.push_back(fObservables[i]);
   78 
   79         if (fObservables[i].find("Qhigh_Y_R") != string::npos)
   80             fQhigh_Y_Observables.push_back(fObservables[i]);
   81         if (fObservables[i].find("Qlow_Y_R") != string::npos) fQlow_Y_Observables.push_back(fObservables[i]);
   82         if (fObservables[i].find("Qbalance_Y_R") != string::npos)
   83             fQbalance_Y_Observables.push_back(fObservables[i]);
   84         if (fObservables[i].find("Qratio_Y_R") != string::npos)
   85             fQratio_Y_Observables.push_back(fObservables[i]);
   86     }
   87 
   88     for (unsigned int i = 0; i < fQ1_Observables.size(); i++) {
   89         Double_t r1 = atof(fQ1_Observables[i].substr(4, fQ1_Observables[i].length() - 4).c_str());
   90         fQ1_Radius.push_back(r1);
   91     }

So, the observable here is for example Q1_R5, where 5 is assigned to the fQ1_Radius.

9

Even a simpler way to do that is to create a Double_t observable for the VETO that registers the maximum veto signal amplitude, or the average amplitude. You could then apply directly vetoMaxSignal<500.

Another interesting observable could be the number of signals that go above a given threshold, a fThreshold given through a metadata parameter.

Of course, this could be done differentiating the veto groups also, top, bottom, …

10

I modified TRestRawVetoAnalysisProcess.
Now you can add in the rml a threshold parameter like

<parameter name="threshold" value="500" />

which will add two observables to the analysis tree:

  • VetoAboveThreshold is either 0, if no veto was above threshold, or 1, if at least one veto signal was above threshold.
  • NVetoAboveThreshold gives the number of veto signals above threshold.

In principle the first observable is not necessary, because one could make a cut on NVetoAboveThreshold==0.

2 Likes
11

I placed myself the question at the root-forum. So, does the answer here solve the std::map cut issue?

12

I finally asked in the root forum how to apply a range cut in the selection string of TTree::Draw() with a std::map observable:

This can be achieved by using logical operations in the cut string, without changing TRestAnalysisPlot. For example the following cut works very well:

<globalCutString string="Sum$(veto_PeakTime.second > 200 &amp;&amp; veto_PeakTime.second < 300) == 0" value="ON" />

In this example, the selection requirement is that none of the vetoes has a signal with a peak time in the (200,300) range.

When using the globalCutString-parameter, so far we have to write &amp;&amp; instead of &&, but I asked here in the forum if this can be changed: Issue when using XML reserved characters in the rml file

We cannot apply a cut on a single veto panel, because TTree::Draw() doesn’t support the std::map:operator[], which would allow to select a signle veto ID (e.g. veto_PeakTime[1234] is the peak Time of veto with ID 1234). To achieve this again the recommendation is to use RDataFrame.

Anyways, I think what we have now is good enough (for now).

1 Like
13

I’ve finally implemented this:

<globalCut variable="veto_MaxPeakAmplitude[4645]" condition="<500" value="ON" />

kind of syntax for the expression evaluation process. Once it’s merged you’ll be able to either select specific channels and store them in their own observable or filter on some condition of them. Relevant PR with examples here:

(might be recreated to merge not from my own fork).

1 Like