/[PAMELA software]/PamelaLevel2/doc/examples-1/Trigger.cpp
ViewVC logotype

Annotation of /PamelaLevel2/doc/examples-1/Trigger.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Mon Nov 5 13:10:12 2007 UTC (17 years ago) by pam-fi
Branch: MAIN
CVS Tags: v10RED, v5r00, v6r00, v9r00, HEAD
implemente gpamela and new examples

1 pam-fi 1.1
2     #if !defined(__CINT__) || defined(__MAKECINT__)
3    
4    
5     // ROOT includes //
6     #include <TString.h>
7     #include <TH1F.h>
8     #include <TH2F.h>
9    
10     // standard libraries includes //
11     #include <stdlib.h>
12     #include <iostream>
13     //#include <fstream>
14     #include <iomanip>
15     using namespace std;
16    
17     // PAMELA includes //
18     #include <PamLevel2.h>
19    
20     #endif
21    
22     //======================
23     // global variables
24     //======================
25     Int_t tot=0;
26     Int_t sel[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
27    
28     TrkLevel2 *trkl2 = 0;
29     CaloLevel2 *calol2 = 0;
30     ToFLevel2 *tofl2 = 0;
31     TrigLevel2 *trigl2 = 0;
32     AcLevel2 *acl2 = 0;
33    
34     //======================
35     // HISTOGRAMS declaration
36     //======================
37    
38     TH2F* hxytof[6];
39    
40     void Dump(ToFTrkVar* tof){
41    
42     cout << endl<< "n. tdc hits :"<<tof->npmttdc;
43     cout << endl<< "PMTs :";
44     for(Int_t ip=0; ip< tof->npmttdc; ip++)cout << " "<< tof->pmttdc[ip];
45     cout << endl<< "n. adc hits :"<<tof->npmtadc;
46     cout << endl<< "PMTs :";
47     for(Int_t ip=0; ip< tof->npmtadc; ip++)cout << " "<< tof->pmtadc[ip];
48     cout << endl<< "dedx :";
49     for(Int_t ib=0; ib< tof->npmtadc; ib++)cout << " "<< tof->dedx[ib];
50     cout << endl<< "beta :";
51     for(Int_t ib=0; ib<13 ; ib++)cout << " "<< tof->beta[ib];
52     cout << endl<< "xtofpos :";
53     for(Int_t ib=0; ib<3 ; ib++)cout << " "<< tof->xtofpos[ib];
54     cout << endl<< "ytofpos :";
55     for(Int_t ib=0; ib<3 ; ib++)cout << " "<<tof->ytofpos[ib];
56     cout <<endl<<"============================================="<<endl;
57    
58     }
59    
60    
61     void Dump(ToFLevel2* l2){
62     cout<<"ToF event: "<<l2->npmt()<<" hit PMTs"<<endl;
63     cout<<" id name plane paddle pmt adc tdc"<<endl;
64     for(int ip=0; ip<l2->npmt(); ip++){
65     Int_t iplane = -1;
66     Int_t ipaddle = -1;
67     Int_t ipmt = -1;
68     // cout << l2->GetToFPMT(ip) << endl;
69     Int_t id = l2->GetToFPMT(ip)->pmt_id;
70     float adc = l2->GetToFPMT(ip)->adc;
71     float tdc = l2->GetToFPMT(ip)->tdc_tw; // temporaneo
72     cout <<setw(3)<<id<< setw(10) << l2->GetPMTName(id,iplane,ipaddle,ipmt) << setw(7) << iplane << setw(7)<< ipaddle << setw(7) << ipmt << setw(7) << adc << setw(7) << tdc <<endl;
73     }
74     cout <<"============================================="<<endl;
75    
76     }
77    
78     void Dump(TrigLevel2* l2){
79     cout << "Trigger Configuration (hex): " << hex;
80     cout << l2->trigconf << endl;
81     cout <<"============================================="<<endl;
82     cout << "Trigger Pattern (hex): " << hex;
83     for(int i=0; i<6; i++)cout << setw(5) << l2->patterntrig[i];
84     cout << dec << endl;
85     cout <<"============================================="<<endl;
86     }
87    
88     void EvaluatePatternTrig(UInt_t *patterntrig, ToFLevel2* l2){
89    
90     // TString photoS[48] = {
91     // "S11_1A", "S11_1B", "S11_2A", "S11_2B", "S11_3A", "S11_3B", "S11_4A",
92     // "S11_4B",
93     // "S11_5A", "S11_5B", "S11_6A", "S11_6B", "S11_7A", "S11_7B", "S11_8A",
94     // "S11_8B",
95     // "S12_1A", "S12_1B", "S12_2A", "S12_2B", "S12_3A", "S12_3B", "S12_4A",
96     // "S12_4B", "S12_5A", "S12_5B", "S12_6A", "S12_6B",
97     // "S21_1A", "S21_1B", "S21_2A", "S21_2B",
98     // "S22_1A", "S22_1B", "S22_2A", "S22_2B",
99     // "S31_1A", "S31_1B", "S31_2A", "S31_2B", "S31_3A", "S31_3B",
100     // "S32_1A", "S32_1B", "S32_2A", "S32_2B", "S32_3A", "S32_3B"
101     // };
102     // patterntrig(3) --> S3 [49 to 38] 12 bits: (S32b,S32a,S31b,S31a)
103     // patterntrig(4) --> S2 [37 to 30] 8 bits: (S22b,S22a,S21b,S21a)
104     // patterntrig(5) --> S12 [29 to 18] 12 bits: (S12b,S12a)
105     // patterntrig(6) --> S11 [17 to 2] 16 bits: (S11b,S11a)
106    
107     // UInt_t patterntrig[6];
108     for(int i=0; i<6; i++)*(patterntrig+i)=0;
109    
110     UInt_t word[]={
111     5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, //S11
112     4,4,4,4,4,4,4,4,4,4,4,4, //S12
113     3,3,3,3,3,3,3,3, //S2
114     2,2,2,2,2,2,2,2,2,2,2,2 //S3
115     };
116     UInt_t bit[] ={
117     0,8,1,9,2,10,3,11,4,12,5,13,6,14,7,15, //S11
118     0,6,1,7,2,8,3,9,4,10,5,11, //S12
119     0,2,1,3,4,6,5,7, //S2
120     0,3,1,4,2,5,6,9,7,10,8,11 //S3
121     };
122    
123     for(int ip=0; ip<l2->npmt(); ip++){
124     Int_t iplane = -1;
125     Int_t ipaddle = -1;
126     Int_t ipmt = -1;
127     // cout << l2->GetToFPMT(ip) << endl;
128     Int_t id = l2->GetToFPMT(ip)->pmt_id;
129     float tdc = l2->GetToFPMT(ip)->tdc_tw; // temporaneo
130    
131     if(tdc<4095) *(patterntrig+word[id]) = *(patterntrig+word[id]) | (UInt_t)(pow(2.,(float)bit[id]));
132    
133     }
134    
135    
136    
137    
138     }
139    
140     //===============================================================================
141     //
142     //
143     //
144     //
145     //===============================================================================
146     bool Select( PamLevel2* event ){
147     // return false;
148     TrkParams::SetVerboseMode();
149    
150     //-------------------------------------------
151     // check if the event contains selection info
152     //-------------------------------------------
153     if( event->GetTrkLevel2()==0 ){
154     cout << "Missing TrkLevel2 object "<<endl;
155     return false;
156     }
157     if( event->GetToFLevel2()==0 ){
158     cout << "Missing ToFLevel2 object "<<endl;
159     return false;
160     }
161     if( event->GetCaloLevel2()==0 ){
162     cout << "Missing CaloLevel2 object "<<endl;
163     return false;
164     }
165     if( event->GetTrigLevel2()==0 ){
166     cout << "Missing TrigLevel2 object "<<endl;
167     return false;
168     }
169     if( event->GetAcLevel2()==0 ){
170     cout << "Missing AcLevel2 object "<<endl;
171     return false;
172     }
173     if( event->GetRunInfo()==0 ){
174     cout << "Missing RunInfo object "<<endl;
175     return false;
176     }
177    
178     trkl2 = event->GetTrkLevel2();
179     calol2 = event->GetCaloLevel2();
180     tofl2 = event->GetToFLevel2();
181     trigl2 = event->GetTrigLevel2();
182     acl2 = event->GetAcLevel2();
183    
184     tot++;
185    
186    
187     //------------------------------------------------------------------
188     // Orbital selection
189     //------------------------------------------------------------------
190     bool ORB__OK = false;
191     if(
192     event->GetOrbitalInfo()->Babs>0.22 && //outside SSA
193     true ) ORB__OK = true;
194     if( !ORB__OK )return false;
195     sel[0]++;
196    
197     //---------------------------------------------------------
198     // trigger configuration: TOF 1 (standard)
199     //---------------------------------------------------------
200     bool trigg_patt_ok = false;
201     if ( trigl2->trigconf & (1<<0) ) trigg_patt_ok = true;
202     if(!trigg_patt_ok)return false;
203     sel[1]++;
204    
205     //---------------------------------------------------------
206     // single track
207     //---------------------------------------------------------
208     if( trkl2->GetNTracks()!=1 ) return false; // singola traccia (fisica)
209     event->SetSortingMethod("+CAL"); // usa solo l'info del calorimetro per scegliere l'immagine
210     PamTrack *track = event->GetTrack(0);// qui viene scelta l'"immagine" fisica
211     // if( track->GetScore() <=0 )cout << " >>> Strange track <<< P-score "<< track->GetPScore() << " I-score "<<track->GetIScore()<<endl;
212    
213     if(!track)return false;
214    
215     //----------------------------------------------------------
216     // (CHI**2 standard)
217     //----------------------------------------------------------
218     float p0 = 1.111588e+00;
219     float p1 = 1.707656e+00;
220     float p2 = 1.489693e-01;
221     float def_0 = 0.07;
222     float chi2m025_0 = p0 + fabs(def_0)*p1 + def_0*def_0*p2;
223    
224     //------------------------------------------------------------------
225     // tracker pre-selection
226     //------------------------------------------------------------------
227     TrkTrack *trk = track->GetTrkTrack();
228    
229     float chi2m025 = p0 + fabs(trk->GetDeflection())*p1 + trk->GetDeflection()*trk->GetDeflection()*p2;
230     float chi2m = pow( chi2m025-chi2m025_0+pow(6.4,0.25), 4.);
231    
232     bool TRACK__OK = false;
233     if(
234     // trk->chi2 >0 &&
235     trk->GetNX()>=4 &&
236     trk->GetNY()>=3 &&
237     trk->GetLeverArmX()>=5 &&
238     trk->GetRigidity()> 4. && // unit: GV
239     trk->GetDEDX()<3. && // unit: dE/dX for Z=1 MIP
240     trk->IsInsideCavity() &&
241     trk->chi2 < chi2m &&
242     true ) TRACK__OK = true;
243    
244     if( !TRACK__OK )return false;
245    
246     sel[2]++;
247    
248    
249     //------------------------------------------------------------------
250     // AC
251     //------------------------------------------------------------------
252     if( acl2->IsHit("CARD1-M") || acl2->IsHit("CARD1-E") )return false;
253     if( acl2->IsHit("CARD2-M") || acl2->IsHit("CARD2-E") )return false;
254     if( acl2->IsHit("CARD3-M") || acl2->IsHit("CARD3-E") )return false;
255     if( acl2->IsHit("CARD4-M") || acl2->IsHit("CARD4-E") )return false;
256    
257     if( acl2->IsHit("CAT1-M") || acl2->IsHit("CAT1-E") )return false;
258     if( acl2->IsHit("CAT2-M") || acl2->IsHit("CAT2-E") )return false;
259     if( acl2->IsHit("CAT3-M") || acl2->IsHit("CAT3-E") )return false;
260     if( acl2->IsHit("CAT4-M") || acl2->IsHit("CAT4-E") )return false;
261    
262     // if( acl2->IsHit("CAS1-M") || acl2->IsHit("CAS1-E") )return false;
263     // if( acl2->IsHit("CAS2-M") || acl2->IsHit("CAS2-E") )return false;
264     // if( acl2->IsHit("CAS3-M") || acl2->IsHit("CAS3-E") )return false;
265     // if( acl2->IsHit("CAS4-M") || acl2->IsHit("CAS4-E") )return false;
266     sel[3]++;
267    
268    
269    
270     //------------------------------------------------------------------
271     // CALORIMETER
272     //------------------------------------------------------------------
273     bool NONINT=false;
274     if(
275     track->GetCaloTrack()->qtrack / event->GetCaloLevel2()->qtot > 0.85 &&
276     true)NONINT=true;
277    
278     if(!NONINT)return false;
279     sel[4]++;
280    
281    
282     return true;
283    
284    
285     }
286     //===============================================================
287     // Create histograms
288     //
289     //
290     //
291     //
292     //
293     //===============================================================
294     void CreateHistos( PamLevel2* event , TFile* outf ){
295     cout << "Creating histos..."<<endl;
296     gROOT->cd();//create histos in memory
297     TString hid;
298     TString title;
299    
300     //--------------------------------
301     // histos
302     //--------------------------------
303     for(Int_t i=0; i<6; i++){
304     hid.Form("hxytof%i",i);
305     title.Form("Track coordinates on ToF plane %i",i);
306     hxytof[i] = new TH2F(hid.Data(),title.Data(),6000,-30.,30.,6000,-30.,30.);
307     }
308    
309     //--------------------------------
310     // initialization
311     //--------------------------------
312     cout << "Opening DB connection..."<<endl;
313     event->SetDBConnection();
314    
315     }
316     //===============================================================
317     // Retrieve histograms
318     //
319     //
320     //
321     //
322     //
323     //===============================================================
324     void RetrieveHistos(){
325    
326     gROOT->cd();
327     TString hid;
328    
329     for(Int_t i=0; i<6; i++){
330     hid.Form("hxytof%i",i);
331     hxytof[i] = dynamic_cast<TH2F*>(gDirectory->FindObject(hid.Data()));
332     }
333    
334     }
335     //===============================================================
336     // Get histograms (... non ho capito la differenza con retrieve)
337     //
338     //
339     //
340     //
341     //
342     //===============================================================
343     void GetHistos(){
344    
345     TString hid;
346     }
347     //===============================================================
348     // Fill histograms (called event-by-event)
349     //
350     //
351     //
352     //
353     //
354     //
355     //===============================================================
356     void FillHistos( PamLevel2* event ){
357     // ----------------------------------
358     // retrieve histos
359     // ----------------------------------
360     RetrieveHistos();
361    
362    
363     // ---------------------------------
364     // Retrieve the track
365     // ---------------------------------
366     PamTrack *track = event->GetTrack(0);
367     TrkTrack *trk = track->GetTrkTrack();
368     ToFTrkVar *tof = track->GetToFTrack();
369    
370     // cout << track->GetTrkTrack()->GetRigidity() << endl;
371     // cout << track->GetToFTrack()->beta[12] << endl;
372     // cout << track->GetCaloTrack()->qcyl << endl;
373    
374     // ---------------------------------------------
375     // track coordinates esxtrapolated on tof planes
376     // ---------------------------------------------
377     float ytof[6],xtof[6];
378     for( int i=0; i<6; i++) ytof[i] = tof->ytr_tof[i];
379     for( int i=0; i<6; i++) xtof[i] = tof->xtr_tof[i];
380    
381     for( int i=0; i<6; i++) hxytof[i]->Fill(xtof[i],ytof[i]);
382    
383     Dump(event->GetToFLevel2());
384     Dump(event->GetTrigLevel2());
385     UInt_t patterntrig[6];
386     EvaluatePatternTrig(patterntrig, event->GetToFLevel2());
387     cout << "TDC-Pattern: "<<hex;
388     for(int i=0; i<6; i++)cout << setw(5) << patterntrig[i];
389     cout << dec << endl;
390     cout <<"============================================="<<endl;
391     Dump(tof);
392    
393    
394     }
395     //===============================================================
396     // Save histograms
397     //
398     //
399     //
400     //
401     //
402     //===============================================================
403     void SaveHistos( TFile * outf ){
404    
405     TString hid;
406    
407     gROOT->cd();
408     // ----------------------------------
409     // retrieve histos
410     // ----------------------------------
411     RetrieveHistos();
412    
413     // ----------------------------------
414     // save on file
415     // ----------------------------------
416     outf->cd();
417    
418     for(Int_t i=0; i<6; i++)hxytof[i]->Write();
419    
420     cout << "####------------------------------------------------------------------" <<endl;
421     cout << "#### pre-selected :"<<tot<<endl;
422     if(!tot)return;
423     cout << "#### orbital :"<<sel[0]<< " (" << ((float)sel[0]/(float)tot) <<")"<<endl;
424     cout << "#### trigger configuration :"<<sel[1]<< " (" << ((float)sel[1]/(float)tot) <<")"<<endl;
425     cout << "#### tracker :"<<sel[2]<< " (" << ((float)sel[2]/(float)tot) <<")"<<endl;
426     cout << "#### AC :"<<sel[3]<< " (" << ((float)sel[3]/(float)tot) <<")"<<endl;
427     cout << "#### calo (non-int) :"<<sel[4]<< " (" << ((float)sel[4]/(float)tot) <<")"<<endl;
428     cout << "####------------------------------------------------------------------" <<endl;
429    
430     }

  ViewVC Help
Powered by ViewVC 1.1.23