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