ToFNuclei/doc/Charge_Selection_Mean_Sigma.txt

// define tracking algorithm
  const char* alg = "NUC"; // or "EXT" "STD" "EXTF" "NUC" "NUCEXT" "NUCEXTF"

//....snip ....

// define PamLevel2  and ToFNuclei
  PamLevel2*  event  = new PamLevel2(dir,filelist,"+AUTO");       // << create pamela event
  ToFNuclei *tofn = new ToFNuclei(event,alg);


// ...snip....

// the n-sigma cuts used for ToFNuclei for the B/C analysis
Float_t nsigma_s12=1.5;  //
Float_t nsigma_s2=2.0; //
Float_t nsigma_s3=1.4; //

 Float_t *charge_layer_trk_raw = NULL;
 Float_t zs[9];  // charge in the 6 layers plus 3 mean charge

 char htitle[50],info[50];

  //==============================================================
  // Here I  read nuclei calibration file
  //==============================================================
  cout<<"reading  nuclei_meansigma  calibration file "<<endl;

sprintf(info,"nuclei_meansigma_10th_1.dat");

cout<<"reading meansigma file "<<info<<endl;

ifstream fin1(info);

  // 5 species:  Li, Be, B, C, O

  TF1 *mean0[9][5];
  TF1 *mean1[9][5];
  TF1 *sigma[9][5];

  for(Int_t i=0;i<9;i++) {
    for(Int_t j=0;j<5;j++) {
      sprintf(htitle, "mean0_%d_%d",j+5,i);
      mean0[i][j] = new TF1(htitle,"pol6",-0.5,1.3);
      sprintf(htitle, "mean1_%d_%d",j+5,i);
      mean1[i][j] = new TF1(htitle,"pol1",1.2,3.);
      sprintf(htitle, "sigma_%d_%d",j+5,i);
      sigma[i][j] = new TF1(htitle,"pol1",-0.5,3.);
    }
  }

  Float_t p0,p1,p2,p3,p4,p5,p6,p0a,p1a;
  Int_t izin,ilayin;


//  S11...S32,S1,S2,S3
// Z= 3,4,5,6,8

for(int i=0;i<9;i++) {
for(int j=0;j<5;j++) {
fin1>>ilayin>>izin>>p0>>p1>>p2>>p3>>p4>>p5>>p6>>p0a>>p1a;
cout<<ilayin<<" "<<izin<<" "<<p0<<" "<<p1<<" "<<p2<<" "<<p3<<" "<<p4<<" "<<p5<<" "<<p6<<" "<<p0a<<" "<<p1a<<endl;
mean0[i][j]->SetParameters(p0,p1,p2,p3,p4,p5,p6);
mean1[i][j]->SetParameters(p0a,p1a);
fin1>>ilayin>>izin>>p0>>p1;
cout<<ilayin<<" "<<izin<<" "<<p0<<" "<<p1<<endl;
sigma[i][j]->SetParameters(p0,p1);
                     }
                     }

fin1.close();


// ....... snip.....

  Long64_t nevents = event->GetEntries();
  cout<<" # events "<<nevents<<endl;

  for(Long64_t iev=0; iev<nevents; iev++) {
    event->GetEntry(iev);

    int   npt = event->GetNTracks(alg); 

    if( npt>0 ) { //if n.physical tracks>0

    PamTrack *pamtrack = event->GetTrack(0,alg); //the method selects the 0th physical track, discarding the image      

    ExtTrack* trktrack = pamtrack->GetExtTrack();
    ToFTrkVar* toftrack = pamtrack->GetToFTrack(); 

//========================================================================
//=======================  only  single tracks    ========================
//========================================================================

     if(
       npt ==1 && //the tracks is single
         trktrack->chi2 > 0 && //the fit converged
         trktrack->nstep < 100 &&    
    true){


      chi2=trktrack->chi2;
      defl = trktrack->GetDeflection();
      rig = 1./defl;
      beta = toftrack->CalcBeta(10.,10.,20.);       //beta12

            if(
                chi2 > 0    &&
               trktrack->GetNX() >= 4  &&
               trktrack->GetNY() >= 3  &&
                defl > 0 &&
                defl < 10 &&
                beta_mean>0.2 &&
                beta_mean<2.0 &&
                true
                ) {

//..........snip........
//=============================================================
//==========  ToFNuclei ===== =================================
//=============================================================

  charge_layer_trk_raw = tofn->Get_Charge_ToF_trk_layer_raw();
  for(Int_t j=0;j<6;j++) {
    zs[j]=charge_layer_trk_raw[j];  // S11,S12,...,S32
  }

// the three mean charges for <S1>, <S2>, <S3>
  zs[6] = 1000.0;
  if(zs[0]<1000 && zs[1]<1000)
    zs[6]=(zs[0]+zs[1])/2;
  else if(zs[0]<1000 && zs[1]>=1000)
    zs[6]=zs[0];
  else if(zs[0]>=1000 && zs[1]<1000)
    zs[6]=zs[1];

  zs[7] = 1000.0;
  if(zs[2]<1000 && zs[3]<1000)
    zs[7]=(zs[2]+zs[3])/2;
  else if(zs[2]<1000 && zs[3]>=1000)
    zs[7]=zs[2];
  else if(zs[2]>=1000 && zs[3]<1000)
    zs[7]=zs[3];

  zs[8] = 1000.0;
  if(zs[4]<1000 && zs[5]<1000)
    zs[8]=(zs[4]+zs[5])/2;
  else if(zs[4]<1000 && zs[5]>=1000)
    zs[8]=zs[4];
  else if(zs[4]>=1000 && zs[5]<1000)
    zs[8]=zs[5];

//================================================================
//============ Defining the cut for the raw selection=============
//============ Define charge using selection criteria  ===========
//================================================================

//===================  S12 & S2mean & S3mean  ====================

Int_t iztof1=0;
Int_t icount=0;

for (Int_t ij=0; ij<5; ij++) {  // Li, Be, B, C, O
  if(
     (rig<=pow(10,1.25) &&
      fabs(zs[1]-mean0[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
      fabs(zs[7]-mean0[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig))) &&
      fabs(zs[8]-mean0[8][ij]->Eval(log10(rig)))< (nsigma_s3*sigma[8][ij]->Eval(log10(rig)))
      ) ||
     (rig>pow(10,1.25) &&
      fabs(zs[1]-mean1[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
      fabs(zs[7]-mean1[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig))) &&
      fabs(zs[8]-mean1[8][ij]->Eval(log10(rig)))< (nsigma_s3*sigma[8][ij]->Eval(log10(rig)))
      )
     ){
   icount++;
   iztof1=ij+3;
   if (ij==4) iztof1 = 8;
        }
                      }

if (icount>1)  cout<<"overlap "<<endl;

//=====================  S12 & S2mean  =========================

Int_t iztof2=0;
icount=0;

for (Int_t ij=0; ij<5; ij++) {  // Li, Be, B, C, O
  if(
     (rig<=pow(10,1.25) &&
      fabs(zs[1]-mean0[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
      fabs(zs[7]-mean0[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig)))
      ) ||
     (rig>pow(10,1.25) &&
      fabs(zs[1]-mean1[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
      fabs(zs[7]-mean1[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig)))
      )
     ){
   icount++;
   iztof2=ij+3;
   if (ij==4) iztof2 = 8;
        }
                      }

// .............snip.......

     } // tracking cuts...
    } // fit converged
   }  // ntracks>1
  }   // event loop


1	mayorov	1.1	// define tracking algorithm
2			const char* alg = "NUC"; // or "EXT" "STD" "EXTF" "NUC" "NUCEXT" "NUCEXTF"
3
4			//....snip ....
5
6			// define PamLevel2 and ToFNuclei
7			PamLevel2* event = new PamLevel2(dir,filelist,"+AUTO"); // << create pamela event
8			ToFNuclei *tofn = new ToFNuclei(event,alg);
9
10
11			// ...snip....
12
13			// the n-sigma cuts used for ToFNuclei for the B/C analysis
14			Float_t nsigma_s12=1.5; //
15			Float_t nsigma_s2=2.0; //
16			Float_t nsigma_s3=1.4; //
17
18			Float_t *charge_layer_trk_raw = NULL;
19			Float_t zs[9]; // charge in the 6 layers plus 3 mean charge
20
21			char htitle[50],info[50];
22
23			//==============================================================
24			// Here I read nuclei calibration file
25			//==============================================================
26			cout<<"reading nuclei_meansigma calibration file "<<endl;
27
28			sprintf(info,"nuclei_meansigma_10th_1.dat");
29
30			cout<<"reading meansigma file "<<info<<endl;
31
32			ifstream fin1(info);
33
34			// 5 species: Li, Be, B, C, O
35
36			TF1 *mean0[9][5];
37			TF1 *mean1[9][5];
38			TF1 *sigma[9][5];
39
40			for(Int_t i=0;i<9;i++) {
41			for(Int_t j=0;j<5;j++) {
42			sprintf(htitle, "mean0_%d_%d",j+5,i);
43			mean0[i][j] = new TF1(htitle,"pol6",-0.5,1.3);
44			sprintf(htitle, "mean1_%d_%d",j+5,i);
45			mean1[i][j] = new TF1(htitle,"pol1",1.2,3.);
46			sprintf(htitle, "sigma_%d_%d",j+5,i);
47			sigma[i][j] = new TF1(htitle,"pol1",-0.5,3.);
48			}
49			}
50
51			Float_t p0,p1,p2,p3,p4,p5,p6,p0a,p1a;
52			Int_t izin,ilayin;
53
54
55			// S11...S32,S1,S2,S3
56			// Z= 3,4,5,6,8
57
58			for(int i=0;i<9;i++) {
59			for(int j=0;j<5;j++) {
60			fin1>>ilayin>>izin>>p0>>p1>>p2>>p3>>p4>>p5>>p6>>p0a>>p1a;
61			cout<<ilayin<<" "<<izin<<" "<<p0<<" "<<p1<<" "<<p2<<" "<<p3<<" "<<p4<<" "<<p5<<" "<<p6<<" "<<p0a<<" "<<p1a<<endl;
62			mean0[i][j]->SetParameters(p0,p1,p2,p3,p4,p5,p6);
63			mean1[i][j]->SetParameters(p0a,p1a);
64			fin1>>ilayin>>izin>>p0>>p1;
65			cout<<ilayin<<" "<<izin<<" "<<p0<<" "<<p1<<endl;
66			sigma[i][j]->SetParameters(p0,p1);
67			}
68			}
69
70			fin1.close();
71
72
73
74			// ....... snip.....
75
76			Long64_t nevents = event->GetEntries();
77			cout<<" # events "<<nevents<<endl;
78
79			for(Long64_t iev=0; iev<nevents; iev++) {
80			event->GetEntry(iev);
81
82			int npt = event->GetNTracks(alg);
83
84			if( npt>0 ) { //if n.physical tracks>0
85
86			PamTrack *pamtrack = event->GetTrack(0,alg); //the method selects the 0th physical track, discarding the image
87
88			ExtTrack* trktrack = pamtrack->GetExtTrack();
89			ToFTrkVar* toftrack = pamtrack->GetToFTrack();
90
91			//========================================================================
92			//======================= only single tracks ========================
93			//========================================================================
94
95			if(
96			npt ==1 && //the tracks is single
97			trktrack->chi2 > 0 && //the fit converged
98			trktrack->nstep < 100 &&
99			true){
100
101
102			chi2=trktrack->chi2;
103			defl = trktrack->GetDeflection();
104			rig = 1./defl;
105			beta = toftrack->CalcBeta(10.,10.,20.); //beta12
106
107			if(
108			chi2 > 0 &&
109			trktrack->GetNX() >= 4 &&
110			trktrack->GetNY() >= 3 &&
111			defl > 0 &&
112			defl < 10 &&
113			beta_mean>0.2 &&
114			beta_mean<2.0 &&
115			true
116			) {
117
118			//..........snip........
119			//=============================================================
120			//========== ToFNuclei ===== =================================
121			//=============================================================
122
123			charge_layer_trk_raw = tofn->Get_Charge_ToF_trk_layer_raw();
124			for(Int_t j=0;j<6;j++) {
125			zs[j]=charge_layer_trk_raw[j]; // S11,S12,...,S32
126			}
127
128			// the three mean charges for <S1>, <S2>, <S3>
129			zs[6] = 1000.0;
130			if(zs[0]<1000 && zs[1]<1000)
131			zs[6]=(zs[0]+zs[1])/2;
132			else if(zs[0]<1000 && zs[1]>=1000)
133			zs[6]=zs[0];
134			else if(zs[0]>=1000 && zs[1]<1000)
135			zs[6]=zs[1];
136
137			zs[7] = 1000.0;
138			if(zs[2]<1000 && zs[3]<1000)
139			zs[7]=(zs[2]+zs[3])/2;
140			else if(zs[2]<1000 && zs[3]>=1000)
141			zs[7]=zs[2];
142			else if(zs[2]>=1000 && zs[3]<1000)
143			zs[7]=zs[3];
144
145			zs[8] = 1000.0;
146			if(zs[4]<1000 && zs[5]<1000)
147			zs[8]=(zs[4]+zs[5])/2;
148			else if(zs[4]<1000 && zs[5]>=1000)
149			zs[8]=zs[4];
150			else if(zs[4]>=1000 && zs[5]<1000)
151			zs[8]=zs[5];
152
153			//================================================================
154			//============ Defining the cut for the raw selection=============
155			//============ Define charge using selection criteria ===========
156			//================================================================
157
158			//=================== S12 & S2mean & S3mean ====================
159
160			Int_t iztof1=0;
161			Int_t icount=0;
162
163			for (Int_t ij=0; ij<5; ij++) { // Li, Be, B, C, O
164			if(
165			(rig<=pow(10,1.25) &&
166			fabs(zs[1]-mean0[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
167			fabs(zs[7]-mean0[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig))) &&
168			fabs(zs[8]-mean0[8][ij]->Eval(log10(rig)))< (nsigma_s3*sigma[8][ij]->Eval(log10(rig)))
169			) \|\|
170			(rig>pow(10,1.25) &&
171			fabs(zs[1]-mean1[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
172			fabs(zs[7]-mean1[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig))) &&
173			fabs(zs[8]-mean1[8][ij]->Eval(log10(rig)))< (nsigma_s3*sigma[8][ij]->Eval(log10(rig)))
174			)
175			){
176			icount++;
177			iztof1=ij+3;
178			if (ij==4) iztof1 = 8;
179			}
180			}
181
182			if (icount>1) cout<<"overlap "<<endl;
183
184			//===================== S12 & S2mean =========================
185
186			Int_t iztof2=0;
187			icount=0;
188
189			for (Int_t ij=0; ij<5; ij++) { // Li, Be, B, C, O
190			if(
191			(rig<=pow(10,1.25) &&
192			fabs(zs[1]-mean0[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
193			fabs(zs[7]-mean0[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig)))
194			) \|\|
195			(rig>pow(10,1.25) &&
196			fabs(zs[1]-mean1[1][ij]->Eval(log10(rig)))< (nsigma_s12*sigma[1][ij]->Eval(log10(rig))) &&
197			fabs(zs[7]-mean1[7][ij]->Eval(log10(rig)))< (nsigma_s2*sigma[7][ij]->Eval(log10(rig)))
198			)
199			){
200			icount++;
201			iztof2=ij+3;
202			if (ij==4) iztof2 = 8;
203			}
204			}
205
206			// .............snip.......
207
208			} // tracking cuts...
209			} // fit converged
210			} // ntracks>1
211			} // event loop
212
213