/[PAMELA software]/tracker/flight/TrkNuclei/doc/example/test.C
ViewVC logotype

Contents of /tracker/flight/TrkNuclei/doc/example/test.C

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Tue Jan 20 10:29:57 2009 UTC (15 years, 10 months ago) by pam-fi
Branch: MAIN, tracker
CVS Tags: V00, HEAD
Changes since 1.1: +0 -0 lines
File MIME type: text/plain
tracker utilities

1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TString.h>
3 #include <TH1F.h>
4 #include <TH2F.h>
5 #include <TF1.h>
6 #include <TFile.h>
7 #include <TEventList.h>
8 #include <TMultiGraph.h>
9 #include <TCanvas.h>
10 #include <TCut.h>
11
12 #include <stdlib.h>
13 #include <iostream>
14 #include <iomanip>
15 using namespace std;
16
17 #include <PamLevel2.h>
18 #include <TrkNuclei.h>
19
20 #endif
21
22 //-------------------------------------------------
23 // functions prototypes
24 //-------------------------------------------------
25 float GetDEDX_ToFPlane(int plane, ToFTrkVar *tof);
26 int GetHitPMTsOutsideTrack(ToFLevel2* l2,ToFTrkVar* tof,int plane);
27 bool ToFOK_New(ToFLevel2* tofl2, ToFTrkVar *tof);
28 int GetClustersOusideTrack(int iv, int ip, TrkLevel2 *trkl2, TrkTrack *trk, float lr);
29 int GetPlanesWithClustersOutsideTrack(TrkLevel2 *trkl2, TrkTrack *trk, float lr);
30 bool MultipleTracksInsideTracker(TrkLevel2* trkl2, TrkTrack *trk, ToFLevel2 *tofl2, ToFTrkVar *tof);
31 //--------------------------------------------------------------------------------
32 //
33 //
34 //
35 //
36 //--------------------------------------------------------------------------------
37 TH2F* test(TString ddir, TString list){
38
39 PamLevel2 *pam = new PamLevel2(ddir,list,"+AUTO");
40 TChain *chain = pam->GetPamTree();
41
42 TrkNuclei *trknuc = new TrkNuclei();
43 TH2F *h = new TH2F("h","Z (trk) vs beta",100,0.3,1.1,100,0.,9.);
44 TH2F *hh = new TH2F("hh","Z (trk) vs rigidity",100,0.1,1000,100,0.,9.);
45 TH2F *hhh = new TH2F("hhh","Z vs Z",100,0.,9.,100,0.,9.);
46
47
48 ULong64_t sel=0;
49 ULong64_t tot=0;
50 // --------------------
51 // loop over the events
52 // --------------------
53 for(int iev=0; iev<chain->GetEntries(); iev++){
54 // for(int iev=0; iev<50000; iev++){
55
56 // read the event
57 pam->GetEntry(iev);
58 tot++;
59
60 // retrieve pamela objects
61 ToFLevel2 *tofl2 = pam->GetToFLevel2();
62 TrkLevel2 *trkl2 = pam->GetTrkLevel2();
63
64 ToFTrkVar *tof = pam->GetTrack(0)->GetToFTrack();
65 TrkTrack *trk = pam->GetTrack(0)->GetTrkTrack();
66
67 // -------------------------------------------------------
68 // these are part of the cuts to reject interacting events
69 // developed for antiprotons (slightly modified... search XXX)
70 // -------------------------------------------------------
71 bool INTERACTION = true;
72 if(
73 ToFOK_New(tofl2,tof) &&
74 !MultipleTracksInsideTracker(trkl2,trk,tofl2,tof) &&
75 true)INTERACTION=false;
76 if(INTERACTION)continue;
77
78
79 // -------------------------------------------------------
80 // retrieve the beta (Wolfgang)
81 // -------------------------------------------------------
82 // NB: this loop over the stored tracks IS NOT in the intention of
83 // the tracker classes implementation.
84 // CalcBeta should be a method of the ToFTrkVar class and should
85 // not need "notrack" as input.
86 int notrack;
87 for(notrack=0; notrack<tofl2->ntrk(); notrack++){
88 if( ((ToFTrkVar*)(tofl2->ToFTrk->At(notrack)))->trkseqno == trk->seqno )break;
89 };
90 float b2 = pam->GetToFLevel2()->CalcBeta(notrack,5.,15.,4.); // medium quality
91 float b3 = pam->GetToFLevel2()->CalcBeta(notrack,3.,20.,3.); // high quality
92
93 if( b3<=0. || b3>2. )continue;
94
95 // -------------------------------------------------------
96 // now fill the TrkNuclei class
97 // -------------------------------------------------------
98 trknuc->Reset();
99 if( !trknuc->Set(trk) )continue;
100
101 sel++;
102
103 // -------------------------------------------------------
104 // the class contains a copy of the TrkTrack object
105 // so that you can retrieve any TrkTrack variable/methods
106 // as trknuc->GetTrkTrack()->....
107 // -------------------------------------------------------
108 float dedx = trknuc->GetTrkTrack()->GetDEDX();
109 // -------------------------------------------------------
110 // the methods to get the charge (as a function of beta) are
111 //
112 // TrkNuclei::GetZ_Beta(ip,iv,beta) // single view (x or y)
113 // TrkNuclei::GetZ_Beta(ip,beta) // single plane (average between x and y)
114 // TrkNuclei::GetZ_Beta(beta) // average over all planes
115 //
116 // -------------------------------------------------------
117 float z1 = trknuc->GetZ_Beta( b3 );
118 // -------------------------------------------------------
119 // the calibration mip-to-charge is implemented as a
120 // TF2 static object, function of beta and dedx.
121 // You can access directly this function as follows:
122 // -------------------------------------------------------
123 float z2 = TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(b3,dedx);
124 // -------------------------------------------------------
125 // or you can change the calibration function as you like,
126 // with the following methods (must be called once):
127 //
128 // TrkNuclei_parameters::Get()->SetZ_Beta(TF2 *myfunc);
129 //
130 // -------------------------------------------------------
131
132
133 // -------------------------------------------------------
134 // Also implemented is the mip-to-charge as a function
135 // of rigidity.
136 //
137 // TrkNuclei::GetZ_Rigidity(ip,iv,rig) // single view (x or y)
138 // ecc...
139 //
140 // At present however the TF2 object is NULL.
141 // You can set your own function as explained above
142 //
143 // TrkNuclei_parameters::Get()->SetZ_Rigidity(TF2 *myfunc);
144 //
145 // -------------------------------------------------------
146 float rig = trknuc->GetTrkTrack()->GetRigidity();
147
148 float z3 = trknuc->GetZ_Rigidity( rig );
149 float z4 = TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,dedx);
150
151 h->Fill(b3,z1);
152 hh->Fill(rig,z3);
153 hhh->Fill(z1,z3);
154
155 cout << setw(10) << b3 <<setw(10) << z1 <<setw(10)<<1./rig<<setw(10)<< z3 << endl;
156
157
158 }
159
160 cout << endl ;
161 cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-"<<endl;
162 cout << sel<<" selected events over "<<tot<<endl;
163 cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-"<<endl;
164
165 // return h;
166 // return hh;
167 return hhh;
168
169 }
170 /////////////////////////////////////////////////////////////////////////////
171 //
172 //
173 // Selection function from antiproton analysis (to be reviewed!!!!)
174 //
175 //
176 /////////////////////////////////////////////////////////////////////////////
177 //----------------------------------------------------------------------------
178 //
179 //----------------------------------------------------------------------------
180 float GetDEDX_ToFPlane(int plane, ToFTrkVar *tof){
181
182 if(!tof)return 0.;
183 float dedx = 0.;
184 float ndedx = 0;
185 ToFLevel2 l2 = ToFLevel2();
186 for(int ip=0; ip<tof->npmtadc; ip++){
187
188 Int_t iplane = -1;
189 Int_t ipaddle = -1;
190 Int_t ipmt = -1;
191 Int_t id = tof->pmtadc[ip];
192
193 l2.GetPMTName(id,iplane,ipaddle,ipmt);
194
195 if(iplane!=plane)continue;
196 dedx += ( tof->adcflag[ip]==0 )*tof->dedx[ip];
197 ndedx += ( tof->adcflag[ip]==0 );
198 }
199 if(ndedx>0)return dedx/ndedx;
200 return 0.;
201 }
202 //--------------------------------------------------------------------------------
203 //
204 // conta il numero di PMT con segnale (tdc) fuori dalla traccia
205 //
206 //--------------------------------------------------------------------------------
207 int GetHitPMTsOutsideTrack(ToFLevel2* l2,ToFTrkVar* tof,int plane){
208 int nn=0;
209 // cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"<<endl;
210 for(int ip=0; ip<l2->npmt(); ip++){
211 Int_t iplane = -1;
212 Int_t ipaddle = -1;
213 Int_t ipmt = -1;
214 Int_t id = l2->GetToFPMT(ip)->pmt_id;
215 l2->GetPMTName(id,iplane,ipaddle,ipmt);
216
217 if(iplane!=plane)continue;
218 // cout << endl<<plane << " -- "<<id;
219 float tdc = l2->GetToFPMT(ip)->tdc; if(tdc >= 4095.)continue;///to avoid adc pile-up
220
221 int iht = -1;
222 int iha = -1;
223 for(iht=0; iht<tof->npmttdc; iht++)if(id==tof->pmttdc[iht])break;
224 if( iht>=0 && iht <tof->npmttdc )continue;
225 for(iha=0; iha<tof->npmtadc; iha++)if(id==tof->pmtadc[iha])break;
226 if( iha>=0 && iha <tof->npmtadc )continue;
227 nn++;
228 // cout << " OUTSIDE ";
229 }
230 return nn;
231 }
232 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
233 //
234 //
235 //
236 //
237 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
238 //
239 // it requires:
240 // - not more than 1 hit paddle on S11/S12/S21/S22
241 // - at least 1 hit paddle on S1/S2
242 // - if there is a hit paddle (S1/S2), this must be associated to the track
243 //
244 //--------------------------------------------------------------------------------
245 bool ToFOK_New(ToFLevel2* tofl2, ToFTrkVar *tof){
246
247 if( !tofl2)return false;
248
249 // cout << "((( TOFOK )))"<<endl;
250
251 bool TOF__OK=false;
252
253 // -----------------------------------------
254 // at least one hit paddle on both S1 and S2
255 // -----------------------------------------
256 if(
257 //-----------------------------------
258 // not more than 1 paddle per layer
259 //-----------------------------------
260 tofl2 &&
261 tofl2->GetNHitPaddles(0) <= 1 && //S11
262 tofl2->GetNHitPaddles(1) <= 1 && //S12
263 tofl2->GetNHitPaddles(2) <= 1 && //S21
264 tofl2->GetNHitPaddles(3) <= 1 && //S22
265 // tofl2->GetNHitPaddles(4) <= 1 && //S31
266 // tofl2->GetNHitPaddles(5) <= 1 && //S32
267 //-----------------------------------
268 // 1 paddle per layer
269 //-----------------------------------
270 // tofl2->GetNHitPaddles(0) == 1 && //S11
271 // tofl2->GetNHitPaddles(1) == 1 && //S12
272 // tofl2->GetNHitPaddles(2) == 1 && //S21
273 // tofl2->GetNHitPaddles(3) == 1 && //S22
274 // tofl2->GetNHitPaddles(4) == 1 && //S31
275 // tofl2->GetNHitPaddles(5) == 1 && //S32
276 //-----------------------------------
277 // at least 1 paddle per plane
278 //-----------------------------------
279 (tofl2->GetNHitPaddles(0)==1 || tofl2->GetNHitPaddles(1)==1) && //S1
280 (tofl2->GetNHitPaddles(2)==1 || tofl2->GetNHitPaddles(3)==1) && //S2
281 // (tofl2->GetNHitPaddles(4)==1 || tofl2->GetNHitPaddles(5)==1) && //S3
282 //-----------------------------------
283 // n.hit pmts outside track
284 //-----------------------------------
285 // GetHitPMTsOutsideTrack(tofl2,tof,0)<2 && //S11 //troppo (puo` succedere )
286 // GetHitPMTsOutsideTrack(tofl2,tof,1)<2 && //S12
287 // GetHitPMTsOutsideTrack(tofl2,tof,2)<2 && //S21
288 // GetHitPMTsOutsideTrack(tofl2,tof,3)<2 && //S22
289 GetHitPMTsOutsideTrack(tofl2,tof,0)<3 && //S11
290 GetHitPMTsOutsideTrack(tofl2,tof,1)<3 && //S12
291 // GetHitPMTsOutsideTrack(tofl2,tof,2)<3 && //S21
292 // GetHitPMTsOutsideTrack(tofl2,tof,3)<3 && //S22
293 true ) TOF__OK = true;
294 // cout << "////"<<tofl2->GetNHitPaddles(0)<<endl;
295 // cout << "////"<<tofl2->GetNHitPaddles(1)<<endl;
296 // cout << "////"<<tofl2->GetNHitPaddles(2)<<endl;
297 // cout << "////"<<tofl2->GetNHitPaddles(3)<<endl;
298 // cout << "////"<<tofl2->GetNHitPaddles(4)<<endl;
299 // cout << "////"<<tofl2->GetNHitPaddles(5)<<endl;
300
301 if( !TOF__OK )return false;
302
303 // ----------------------------------------
304 // dopo che wolfgang mi ha illuminato....
305 // ----------------------------------------
306 int hitplane[]={0,0,0,0,0,0};
307 int hitpaddle[]={-1,-1,-1,-1,-1,-1};
308 for(Int_t ip=0; ip< tof->npmttdc; ip++){
309 Int_t iplane = -1;
310 Int_t ipaddle = -1;
311 Int_t ipmt = -1;
312 Int_t id = tof->pmttdc[ip];
313 tofl2->GetPMTName(id,iplane,ipaddle,ipmt);
314 if(tof->tdcflag[ip]==0){
315 hitplane[iplane]=1; //there is a true tdc signal associated to the track
316 if( hitpaddle[iplane]>=0 && hitpaddle[iplane]!=ipaddle )cout<< "ORRORE!!!!"<<endl;
317 hitpaddle[iplane]=ipaddle;//store the id of the paddle associated to the track
318 }
319 };
320 //
321 for(int iplane=0; iplane<4; iplane++){//loop over S11/S12/S21/S22
322 //retrieve the id of the paddle traversed by the track
323 // (...there is some tolerance...)
324 int ipaddle = tofl2->GetPaddleIdOfTrack(tof->xtr_tof[iplane],tof->ytr_tof[iplane],iplane);
325 //check if the traversed paddle is hit
326 bool OK = true;
327 if(
328 tofl2->GetNHitPaddles(iplane)>0 &&//if there is a hit paddle in this plane...
329 !tofl2->HitPaddle(iplane,ipaddle) && //...and the paddle traversed by the track is not hit
330 true)OK = false;//..discard the event
331 /// hence try to recover some events...
332 if(
333 tofl2->GetNHitPaddles(iplane)>0 &&//if there is a hit paddle...
334 ipaddle==-1 &&//...and the track does not traverse any paddle...
335 tofl2->HitPaddle(iplane,hitpaddle[iplane]) &&//... BUT there are tdc signals belonging to a hit paddle
336 true)OK=true;//recover
337 TOF__OK = TOF__OK && OK;
338 }
339
340 if( !TOF__OK )return false;
341
342
343 return true;
344
345 }
346 ///////////////////////////////////////////////////////////
347 // conta i cluster fuori traccia sulla singola vista
348 // se lr=+1(-1) conta solo quelli a destra(sinistra) della traccia
349 // vengono esclusi dal conteggio:
350 // - i cluster bad
351 // - i cluster saturati (...misembrava ragionevole)
352 // - i cluster entro toll dalla traccia
353 // - i cluster con signal < dedxcut
354 //
355 int GetClustersOusideTrack(int iv, int ip, TrkLevel2 *trkl2, TrkTrack *trk, float lr){
356
357 if(!trkl2)return -1;
358 if(ip<-1||ip>5)return -1;
359 if(lr!=0 && trk==0)return -1;
360
361
362 float dedxcut = 0.7;//mip
363 float toll = 0.01;//cm (100 micron)
364
365 int nhits=0;
366 int nsing=0;
367
368 bool LEFT = false;
369 bool RIGHT = false;
370 if(lr>0)RIGHT=true;
371 if(lr<0)LEFT=true;
372
373 //
374 if(iv==0)nsing=trkl2->nclsx();
375 if(iv==1)nsing=trkl2->nclsy();
376
377 float trkcoord = 0.;
378 if(trk && iv==0)trkcoord=trk->xv[ip];
379 if(trk && iv==1)trkcoord=trk->yv[ip];
380
381 //
382 for(int ii=0; ii<nsing; ii++){
383 TrkSinglet *s = 0;
384 if(iv==0)s = trkl2->GetSingletX(ii);
385 if(iv==1)s = trkl2->GetSingletY(ii);
386 if(!s)continue;
387 if(s->plane != ip+1)continue;
388 // cout << ii << s->plane << s->IsBad() << endl;
389 // if(!s->IsBad())cout << s->sgnl;
390 // cout << s->plane << " - "<< (s->coord[0] - trkcoord) << " "<<(s->coord[1] - trkcoord)<<" "<<s->sgnl<<endl;
391 bool SIDEOK=true;
392 if(
393 LEFT &&
394 ((s->coord[0] - trkcoord)> -1*toll
395 ||
396 (s->coord[1] - trkcoord)> -1*toll) &&
397 true)SIDEOK=false; ;
398 if(
399 RIGHT &&
400 ((s->coord[0] - trkcoord)< +1*toll
401 ||
402 (s->coord[1] - trkcoord)< +1*toll) &&
403 true)SIDEOK=false; ;
404
405 if(
406 // (ip==-1 || s->plane == ip+1) &&
407 !s->IsBad() &&
408 s->GetSignal() > dedxcut &&
409 !s->IsSaturated() &&
410 SIDEOK &&
411 true){
412 // cout << "*";
413 nhits++;
414 }
415 /// cout << endl;
416 }
417 // cout << "-----"<<endl;
418 return nhits;
419
420 }
421 //----------------------------------------------------------------------------
422 //
423 //----------------------------------------------------------------------------
424 int GetPlanesWithClustersOutsideTrack(TrkLevel2 *trkl2, TrkTrack *trk, float lr){
425 int np=0;
426 for(int ip=0; ip<6; ip++)
427 np+=(GetClustersOusideTrack(0,ip,trkl2,trk,lr)>0||GetClustersOusideTrack(1,ip,trkl2,trk,lr)>0);
428 return np;
429 }
430 ///////////////////////////////////////////////////////////
431 bool MultipleTracksInsideTracker(TrkLevel2* trkl2, TrkTrack *trk, ToFLevel2 *tofl2, ToFTrkVar *tof){
432
433
434 if(!trkl2)return true;
435 if(!tofl2)return true;
436
437 if(!trk)return true;
438 if(!tof)return true;
439
440 // cout << "MultipleTracksInsideTracker ?"<<endl;
441 float rig = trk->GetRigidity();
442
443 // n.PMTs outside the track
444 int npmt = 0;
445 if(tofl2&&tof)
446 for(int ip=0; ip<4; ip++)npmt+=GetHitPMTsOutsideTrack(tofl2,tof,ip);
447
448 // n.hits on the left/right of the track (x-view)
449 int npxl=0;
450 int npxr=0;
451 for(int ip=0; ip<5; ip++)
452 npxl+=(GetClustersOusideTrack(0,ip,trkl2,trk,-1)>0) ;
453 for(int ip=0; ip<5; ip++)
454 npxr+=(GetClustersOusideTrack(0,ip,trkl2,trk,+1)>0) ;
455
456 // n.hits outside the track (y-view)
457 int npy=0;
458 for(int ip=0; ip<5; ip++)
459 npy+=(GetClustersOusideTrack(1,ip,trkl2,trk,0)>0) ;
460
461 // upper limit on tof dedx
462 bool troppo1=false;
463 bool troppo2=false;
464 //
465 //
466 // ------------------------------------------------------
467 // XXX
468 // commented to adapt the routine to nuclei analysis !!!
469 // ------------------------------------------------------
470 //
471 // float cuts11up = F_ftofdedxcut(0,rig);
472 // float cuts12up = F_ftofdedxcut(1,rig);
473 // if( tofl2 && GetDEDX_ToFPlane(0,tof) > cuts11up)troppo1=true;
474 // if( tofl2 && GetDEDX_ToFPlane(1,tof) > cuts12up)troppo1=true;
475 // //
476 // float cuts21up = F_ftofdedxcut(2,rig);
477 // float cuts22up = F_ftofdedxcut(3,rig);
478 // if( tofl2 && GetDEDX_ToFPlane(2,tof) > cuts21up)troppo2=true;
479 // if( tofl2 && GetDEDX_ToFPlane(3,tof) > cuts22up)troppo2=true;
480
481 // multiple tracks
482 bool INTERACTING = false;
483 if( (npxl>2 || npxr>2) && npy>2)INTERACTING=true;
484 if( (npxl>1 || npxr>1) && npy>1 && (npmt!=0||troppo2) )INTERACTING=true;
485
486 // cout << npxl << " "<<npxr<< " "<<npy<<endl;
487 // cout << npmt << endl;
488 // if( (npxl>2 || npxr>2) && npy>2)cout << " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 33"<<endl;
489 // if( (npxl>1 || npxr>1) && npy>1 && npmt!=0 )cout << " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% npmt !!"<<endl;
490 // if( (npxl>1 || npxr>1) && npy>1 && troppo2 )cout << " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% troppo2 !!"<<endl;
491
492 return INTERACTING;
493
494 }

  ViewVC Help
Powered by ViewVC 1.1.23