/[PAMELA software]/PamelaLevel2/doc/examples/My-Histos.cpp
ViewVC logotype

Diff of /PamelaLevel2/doc/examples/My-Histos.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by pam-fi, Mon Dec 11 18:26:55 2006 UTC revision 1.3 by pam-fi, Mon Jan 15 11:51:39 2007 UTC
# Line 5  Line 5 
5  //  It contains the functions:  //  It contains the functions:
6  //  //
7  //  - CreateHistos()  //  - CreateHistos()
8    //  - RetrieveHistos()
9  //  - FillHistos(PamLevel2* event)  //  - FillHistos(PamLevel2* event)
10  //  - Save Histos()  //  - Save Histos()
11  //  //
# Line 38  using namespace std; Line 39  using namespace std;
39  //======================  //======================
40  // HITOGRAMS declaration  // HITOGRAMS declaration
41  //======================  //======================
42  TH1F *hresangx;  
43  TH1F *hresangy;  TH2F *hchi2def;
44    TH2F *hchi2rig;
45    TH2F *hdedxtofvsdedxtrk;
46    TH2F *hdedxtrkvsrig ;
47    TH2F *hdedxtofvsrig ;
48    TH2F *hbetavsrig ;
49    TH2F *hdedxtrkvsdef ;
50    TH2F *hbetavsdef ;
51    TH2F *hbetavsdef2 ;
52    
53    TH1F *herrdef;
54    TH1F *herrx[6];
55    TH1F *herry[6];
56    
57    TH1F *hdef ;
58    TH1F *hrig1 ;
59    TH1F *hrig2 ;
60    TH1F *hen1 ;
61    TH1F *hen2 ;
62    
63    
64    Int_t tot=0;
65    Int_t sel[]={0,0,0,0,0,0,0,0};
66    
67    //===============================================================
68    // Get ToF dE/dx
69    //===============================================================
70    // Float_t GetdEdx( PamTrack *trk ){
71    //     Float_t dedx = 0.;
72    //     Int_t nval=0;
73    //     Bool_t ISNULL = false;
74    //     for (Int_t i=0; i<trk->GetToFTrack()->npmtadc; i++){
75    //      if(trk->GetToFTrack()->dedx[i] == 0)ISNULL=true;        
76    //      dedx += (trk->GetToFTrack()->dedx).At(i)*(Int_t)(!ISNULL);
77    //      nval += (Int_t)(!ISNULL);
78    //     };
79    //     //
80    //     if(nval)dedx=dedx/nval;
81    // //    cout <<" -- dedx  "<<dedx<<endl;
82    //     return(dedx);    
83    // }
84    
85    // float GetZ(float beta, float dedx){
86    
87    //     int nspec=2;
88    //     float p0[]={0.904,3.373};
89    //     float p1[]={0.295,2.756};
90    //     float zz[]={1.,4.};
91    
92    //     if(beta>1)beta=1;
93    //     float f0 = p0[0]/beta+p1[0];
94    //     float f1 = p0[1]/beta+p1[1];
95    //     float z0 = zz[0];
96    //     float z1 = zz[1];
97    //     float aa = (z0-z1)/(f0-f1);
98    //     float bb = -aa*f0+z0;
99    //     if( (bb+aa*dedx)==0 )return 0;
100    //     return sqrt(bb+aa*dedx);
101    
102    // }
103    
104    
105  //===============================================================  //===============================================================
106  // Create histograms  // Create histograms
107  //===============================================================  //===============================================================
# Line 47  void CreateHistos( TFile* outf ){ Line 109  void CreateHistos( TFile* outf ){
109    
110      gROOT->cd();//create histos in memory      gROOT->cd();//create histos in memory
111    
112      hresangx = new TH1F("hresangx","Angular residual on the top calo plane (x)",400,-40.,40.);      TString hid;
113      hresangy = new TH1F("hresangy","Angular residual on the top calo plane (y)",400,-40.,40.);      TString title;
114        //===========================================
115        // rigidity binning
116        Int_t   nbins = 200;
117        Float_t rmin  = .1;
118        Float_t rmax  = 1000.;
119        Float_t dr = (log10(rmax)-log10(rmin))/nbins;
120        TArrayD* arr = new TArrayD(nbins+1);
121        arr->AddAt((Double_t)rmin,0);
122        for(Int_t i=1; i<=nbins; i++){
123            Double_t previous = arr->At(i-1);
124            Double_t value=0;
125            value = (Double_t) pow(10,log10((Float_t)previous)+dr);
126            arr->AddAt(value,i);
127        };
128        Double_t *rbinpt = arr->GetArray();
129        //===========================================
130    
131        hchi2def = new TH2F("hchi2def","(chi**2)**0.25 vs deflection",100,0.,2.5,200,0.,10.);
132        hchi2rig = new TH2F("hchi2rig","#chi^{2} vs R",nbins,rbinpt,1000,0.,500.);
133        hdedxtofvsdedxtrk  = new TH2F("hdedxtofvsdedxtrk","dE/dx (ToF-average) vs dE/dx (trk-average)",600,0.,40.,600,0.,40.);
134        hdedxtrkvsrig      = new TH2F("hdedxtrkvsrig","dE/dx (Trk-average) vs R",nbins,rbinpt,600,0.,40.);
135        hdedxtofvsrig      = new TH2F("hdedxtofvsrig","dE/dx (ToF-average) vs R",nbins,rbinpt,600,0.,40.);
136        hbetavsrig         = new TH2F("hbetavsrig","beta vs R",nbins,rbinpt,200,0.4,1.5);
137        hrig1              = new TH1F("hrig1","Rigidity distribution Z=1",nbins,rbinpt);
138        hrig2              = new TH1F("hrig2","Rigidity distribution Z=2",nbins,rbinpt);
139        hen1               = new TH1F("hen1","Kin.Energy per nucleon Z=1",nbins,rbinpt);
140        hen2               = new TH1F("hen2","Kin.Energy per nucleon Z=2",nbins,rbinpt);
141    
142        hdedxtrkvsdef      = new TH2F("hdedxtrkvsdef","dE/dx (Trk-average) vs deflection",500,-10.,10.,600,0.,40.);
143        hbetavsdef         = new TH2F("hbetavsdef","beta vs deflection",500,-10.,10.,200,0.4,1.5);
144        hbetavsdef2        = new TH2F("hbetavsdef2","beta vs deflection (Z=1)",500,-10.,10.,200,0.4,1.5);
145        hdef               = new TH1F("hdef","Deflection distribution",5000,-2.5,2.5);
146    
147        herrdef            = new TH1F("herrdef","Deflection error",1000,0.,0.01);
148        
149        for(Int_t ip=0; ip<6; ip++){
150            hid.Form("herrx%i",ip);
151            title = ""; title.Form("Spatial resolution - plane %i (x)",ip);
152            herrx[ip] = new TH1F(hid.Data(),title.Data(),500,0.,0.005);
153            hid.Form("herry%i",ip);
154            title = ""; title.Form("Spatial resolution - plane %i (y)",ip);
155            herry[ip] = new TH1F(hid.Data(),title.Data(),500,0.,0.005);
156        }
157    
158  }  }
159    //===============================================================
160    // Retrieve histograms
161    //===============================================================
162    void RetrieveHistos(){
163    
164        TString hid;
165    
166        hchi2def           = dynamic_cast<TH2F*>(gDirectory->FindObject("hchi2def"));
167        hchi2rig           = dynamic_cast<TH2F*>(gDirectory->FindObject("hchi2rig"));
168        hdedxtofvsdedxtrk  = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtofvsdedxtrk"));
169        hdedxtrkvsrig      = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtrkvsrig"));
170        hdedxtofvsrig      = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtofvsrig"));
171        hbetavsrig         = dynamic_cast<TH2F*>(gDirectory->FindObject("hbetavsrig"));
172        hdedxtrkvsdef      = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtrkvsdef"));
173        hbetavsdef         = dynamic_cast<TH2F*>(gDirectory->FindObject("hbetavsdef"));
174        hbetavsdef2        = dynamic_cast<TH2F*>(gDirectory->FindObject("hbetavsdef2"));
175        hdef               = dynamic_cast<TH1F*>(gDirectory->FindObject("hdef"));
176        hrig1              = dynamic_cast<TH1F*>(gDirectory->FindObject("hrig1"));
177        hrig2              = dynamic_cast<TH1F*>(gDirectory->FindObject("hrig2"));
178        hen1               = dynamic_cast<TH1F*>(gDirectory->FindObject("hen1"));
179        hen2               = dynamic_cast<TH1F*>(gDirectory->FindObject("hen2"));
180        herrdef            = dynamic_cast<TH1F*>(gDirectory->FindObject("herrdef"));
181        for(Int_t ip=0; ip<6; ip++){
182            hid.Form("herrx%i",ip);
183            herrx[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
184            hid.Form("herry%i",ip);
185            herry[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
186        }
187    }
188  //===============================================================  //===============================================================
189  // Fill histograms (called event-by-event)  // Fill histograms (called event-by-event)
190  //===============================================================  //===============================================================
# Line 60  void FillHistos( PamLevel2* event ){ Line 193  void FillHistos( PamLevel2* event ){
193  //  ----------------------------------  //  ----------------------------------
194  //  retrieve histos  //  retrieve histos
195  //  ----------------------------------  //  ----------------------------------
196  //     hresangx = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangx"));      RetrieveHistos();
 //     hresangy = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangy"));  
197    
198  //  ----------------------------------  //  ----------------------------------
199  //  calo variables  //  single-track events
200  //  ----------------------------------  //  ----------------------------------
201      CaloAxis *x_axis = new CaloAxis();      float beta = 0;
202      CaloAxis *y_axis = new CaloAxis();      float dedxtof = 0;
203        float dedxtrk = 0;
204      float rcil = 1.;// tolerance (cm)      float rigidity = 0;
205      x_axis->FitAxis(event->GetCaloLevel1(),0,rcil);      float deflection = 0;
206      y_axis->FitAxis(event->GetCaloLevel1(),1,rcil);      float chi2 = 0;
207        float errdef = 0.;
208      float qtrack = x_axis->GetQaxis()+y_axis->GetQaxis();      float errx[] = {10.,10.,10.,10.,10.,10.};
209      int   ntrack = x_axis->GetN()+y_axis->GetN();          float erry[] = {10.,10.,10.,10.,10.,10.};
210            
211      float caloangx = (float)atan((double)x_axis->par[1])*180./3.1415026;      tot++;
212      float caloangy = (float)atan((double)y_axis->par[1])*180./3.1415026;          
213        if( event->GetTrkLevel2()->GetNTracks()==0 ) return;
214      CaloStrip top = CaloStrip();      PamTrack *track = event->GetTrack(0);
215      top.Set(0,0,48);      //------------------------------
216      float zcalotop = top.GetZ();      // track selection
217        //------------------------------
218  //  ----------------------------------      if(
219  //  tof variables          true
220  //  ----------------------------------          ){
221  //  get track stored by the tof              
222      ToFTrkVar *tof = event->GetToFStoredTrack(-1);          sel[0]++;
223      if(!tof){          //-------------------
224          cout << " no ToF-track stored "<<endl;          // Calo variables
225          return;          //-------------------
226      }          
227            
228  //  ----------------------------------          //-------------------
229  //  tracker variables          // ToF variables
230  //  ----------------------------------          //-------------------
231            beta    = track->GetToFTrack()->beta[12];
232    //      dedxtof = GetdEdx(track);
233            
234            //-------------------
235            // tracker variables
236            //-------------------
237            dedxtrk    = track->GetTrkTrack()->GetDEDX()/1.3; //<<<ATTENZIONE TEMPORANEO!!!!
238            rigidity   = track->GetTrkTrack()->GetRigidity();
239            deflection = track->GetTrkTrack()->GetDeflection();
240            chi2       = track->GetTrkTrack()->chi2;
241            errdef     = sqrt(track->GetTrkTrack()->coval[4][4]);
242            for(Int_t ip=0; ip<6; ip++){
243                if(track->GetTrkTrack()->XGood(ip))errx[ip]=track->GetTrkTrack()->resx[ip];
244                if(track->GetTrkTrack()->YGood(ip))erry[ip]=track->GetTrkTrack()->resy[ip];
245            }
246    
247      float resangx = -999;          //-------------------
248      float resangy = -999;          // chi2 selection
249      if( event->GetNTracks()==1 ){          //-------------------
250                    hchi2def->Fill( fabs(deflection), pow((double)chi2,0.25) );
251          PamTrack *track = event->GetTrack(0);          hchi2rig->Fill( rigidity, chi2 );
252          if( track->chi2 >0  ){          if(pow((double)chi2,0.25) > (2.5+1.85*fabs(deflection)) )return;
253            sel[1]++;
254              Float_t ztraj[13];          
255                //      hdedxtofvsdedxtrk->Fill(dedxtrk,dedxtof);
256              Int_t i=0;          hdedxtrkvsrig->Fill(rigidity,dedxtrk);
257              ztraj[i++] = event->GetZTOF(11);  //      hdedxtofvsrig->Fill(rigidity,dedxtof);
258              ztraj[i++] = event->GetZTOF(12);          hbetavsrig->Fill(rigidity,beta);
259              ztraj[i++] = event->GetZTOF(21);  
260              ztraj[i++] = event->GetZTOF(22);          hdedxtrkvsdef->Fill(deflection,dedxtrk);
261              ztraj[i++] = event->GetZTrk(1);          hbetavsdef->Fill(deflection,beta);
262              ztraj[i++] = event->GetZTrk(2);  
263              ztraj[i++] = event->GetZTrk(3);          //------------------------
264              ztraj[i++] = event->GetZTrk(4);          // rough charge selection
265              ztraj[i++] = event->GetZTrk(5);          //------------------------
266              ztraj[i++] = event->GetZTrk(6);          if(dedxtrk > (6.+16./(rigidity*rigidity)))return;
267              ztraj[i++] = event->GetZTOF(31);          //--------
268              ztraj[i++] = event->GetZTOF(32);          // helium
269              ztraj[i++] = zcalotop;          //--------
270              Trajectory traj = Trajectory(13,ztraj);          if( dedxtrk > (2.8+3.3/(rigidity*rigidity)) && deflection>0 ){
271              traj.DoTrack2( track->al );              sel[5]++;
272                            hrig2->Fill(rigidity);
273              resangx = caloangx - traj.thx[12];              float m=3.725;
274              resangy = caloangy - traj.thy[12];              float a=4.;
275                float z=2.;
276                float en=(sqrt( (z*rigidity)*(z*rigidity)+m*m )-m)/a;
277                hen2->Fill(en);
278            }
279    
280            //--------
281            // Z=1
282            //--------
283            if(dedxtrk > (2.8+3.3/(rigidity*rigidity)))return;
284            if(dedxtrk < (0.4+0.35/(rigidity*rigidity)))return;
285            sel[2]++;
286    
287            if(deflection>0){
288                hrig1->Fill(rigidity);
289                float m=0.938;
290                float a=1.;
291                float z=1.;
292                float en=(sqrt( (z*rigidity)*(z*rigidity)+m*m )-m)/a;
293                hen1->Fill(en);
294          }          }
295          if(track) delete track;          //------------------------
296            // 6-planes lever-arm
297            //------------------------
298            if( !(track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(5)) )return;
299            sel[3]++;
300    
301            for(Int_t ip=0; ip<6; ip++)herrx[ip]->Fill(errx[ip]);
302            for(Int_t ip=0; ip<6; ip++)herry[ip]->Fill(erry[ip]);
303    
304            //---------------------------
305            // cut on spatial resolution
306            //---------------------------
307            if( errx[0] > 0.0004 )return;
308            if( errx[5] > 0.0004 )return;
309            sel[4]++;
310    
311            herrdef->Fill(errdef);
312            hbetavsdef2->Fill(deflection,beta);
313            hdef->Fill(deflection);
314      }      }
315            
316  //  ----------------------------------  //    delete track;
317  //  fill histograms    
318  //  ----------------------------------      
     hresangx->Fill(resangx);  
     hresangy->Fill(resangy);  
   
   
319  }  }
320  //===============================================================  //===============================================================
321  // Save histograms  // Save histograms
322  //===============================================================  //===============================================================
323  void SaveHistos( TFile * outf ){  void SaveHistos( TFile * outf ){
324    
325        TString hid;
326    
327      gROOT->cd();      gROOT->cd();
328  //  ----------------------------------  //  ----------------------------------
329  //  retrieve histos  //  retrieve histos
330  //  ----------------------------------  //  ----------------------------------
331  //     hresangx = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangx"));      RetrieveHistos();
 //     hresangy = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangy"));  
332    
333  //  ----------------------------------  //  ----------------------------------
334  //  save on file  //  save on file
335  //  ----------------------------------  //  ----------------------------------
   
336      outf->cd();      outf->cd();
337            
338      if(hresangx) hresangx->Write();      hchi2def->Write();
339      if(hresangy) hresangy->Write();      hchi2rig->Write();
340        hdedxtofvsdedxtrk->Write();
341        hdedxtrkvsrig->Write();  
342        hdedxtofvsrig->Write();
343        hbetavsrig->Write();  
344        hdedxtrkvsdef->Write();
345        hbetavsdef->Write();
346        hbetavsdef2->Write();
347        hdef->Write();
348        hrig1->Write();
349        hrig2->Write();
350        hen1->Write();
351        hen2->Write();
352        herrdef->Write();
353        for(Int_t ip=0; ip<6; ip++){
354            herrx[ip]->Write();
355            herry[ip]->Write();
356        }
357    
358        cout << "------------------------------------------------------------------" <<endl;
359        cout << "Pre-selected events :"<<tot<<endl;
360        if(!tot)return;
361        cout << "chi2 selection                            :"<<sel[1]<< "  (" << ((float)sel[1]/(float)tot) <<")"<<endl;
362        cout << "Z=2 selection                             :"<<sel[5]<< "  (" << ((float)sel[5]/(float)tot) <<")"<<endl;
363        cout << "Z=1 selection                             :"<<sel[2]<< "  (" << ((float)sel[2]/(float)tot) <<")"<<endl;
364        cout << "(6-plane lever-arm)                       :"<<sel[3]<< "  (" << ((float)sel[3]/(float)tot) <<")"<<endl;
365        cout << "ext.spatial res. on planes 1-6 < 4 micron :"<<sel[4]<< "  (" << ((float)sel[4]/(float)tot) <<")"<<endl;
366        cout << "------------------------------------------------------------------" <<endl;
367    
368  }  }
369    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23