doc/examples/My-Histos.cpp

////////////////////////////////////////////////////////////////
//
//  Histogram set for tracker-efficiency study.
//
//  It contains the functions:
//
//  - CreateHistos()
//  - RetrieveHistos()
//  - FillHistos(PamLevel2* event)
//  - Save Histos()
//
//  Usage (example):
//  [] .L My-histo-set-macro.C+
//  [] .L My-selection.C+
//  [] .L Loop.C+
//  [] Loop("data-dir/","file-list.txt",100000,"+ALL fillTree fillHistos","output-file.root")
//  
////////////////////////////////////////////////////////////////

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMath.h>
#include <TLine.h>
#include <TPolyMarker.h>
#include <TSelector.h>
#include <TFile.h>

#include <stdlib.h>
#include <iostream>
using namespace std;

#include <PamLevel2.h>
#include <CaloAxis.h>

#endif
//======================
// HITOGRAMS declaration
//======================

TH2F *hchi2def;
TH2F *hchi2rig;
TH2F *hdedxtofvsdedxtrk;
TH2F *hdedxtrkvsrig ;
TH2F *hdedxtofvsrig ;
TH2F *hbetavsrig ;
TH2F *hdedxtrkvsdef ;
TH2F *hbetavsdef ;
TH2F *hbetavsdef2 ;

TH1F *herrdef;
TH1F *herrx[6];
TH1F *herry[6];

TH1F *hdef ;
TH1F *hrig1 ;
TH1F *hrig2 ;
TH1F *hen1 ;
TH1F *hen2 ;


Int_t tot=0;
Int_t sel[]={0,0,0,0,0,0,0,0};

//===============================================================
// Get ToF dE/dx
//===============================================================
// Float_t GetdEdx( PamTrack *trk ){
//     Float_t dedx = 0.;
//     Int_t nval=0;
//     Bool_t ISNULL = false;
//     for (Int_t i=0; i<trk->GetToFTrack()->npmtadc; i++){
//      if(trk->GetToFTrack()->dedx[i] == 0)ISNULL=true;        
//      dedx += (trk->GetToFTrack()->dedx).At(i)*(Int_t)(!ISNULL); 
//      nval += (Int_t)(!ISNULL);
//     };
//     //
//     if(nval)dedx=dedx/nval;
// //    cout <<" -- dedx  "<<dedx<<endl;
//     return(dedx);    
// }

// float GetZ(float beta, float dedx){

//     int nspec=2;
//     float p0[]={0.904,3.373};
//     float p1[]={0.295,2.756};
//     float zz[]={1.,4.};

//     if(beta>1)beta=1;
//     float f0 = p0[0]/beta+p1[0];
//     float f1 = p0[1]/beta+p1[1];
//     float z0 = zz[0];
//     float z1 = zz[1];
//     float aa = (z0-z1)/(f0-f1);
//     float bb = -aa*f0+z0;
//     if( (bb+aa*dedx)==0 )return 0;
//     return sqrt(bb+aa*dedx);

// }


//===============================================================
// Create histograms
//===============================================================
void CreateHistos( TFile* outf ){

    gROOT->cd();//create histos in memory

    TString hid;
    TString title;
    //===========================================
    // rigidity binning
    Int_t   nbins = 200;
    Float_t rmin  = .1;
    Float_t rmax  = 1000.;
    Float_t dr = (log10(rmax)-log10(rmin))/nbins;
    TArrayD* arr = new TArrayD(nbins+1);
    arr->AddAt((Double_t)rmin,0);
    for(Int_t i=1; i<=nbins; i++){
        Double_t previous = arr->At(i-1);
        Double_t value=0;
        value = (Double_t) pow(10,log10((Float_t)previous)+dr);
        arr->AddAt(value,i);
    };
    Double_t *rbinpt = arr->GetArray();
    //===========================================

    hchi2def = new TH2F("hchi2def","(chi**2)**0.25 vs deflection",100,0.,2.5,200,0.,10.);
    hchi2rig = new TH2F("hchi2rig","#chi^{2} vs R",nbins,rbinpt,1000,0.,500.);
    hdedxtofvsdedxtrk  = new TH2F("hdedxtofvsdedxtrk","dE/dx (ToF-average) vs dE/dx (trk-average)",600,0.,40.,600,0.,40.);
    hdedxtrkvsrig      = new TH2F("hdedxtrkvsrig","dE/dx (Trk-average) vs R",nbins,rbinpt,600,0.,40.);
    hdedxtofvsrig      = new TH2F("hdedxtofvsrig","dE/dx (ToF-average) vs R",nbins,rbinpt,600,0.,40.);
    hbetavsrig         = new TH2F("hbetavsrig","beta vs R",nbins,rbinpt,200,0.4,1.5);
    hrig1              = new TH1F("hrig1","Rigidity distribution Z=1",nbins,rbinpt);
    hrig2              = new TH1F("hrig2","Rigidity distribution Z=2",nbins,rbinpt);
    hen1               = new TH1F("hen1","Kin.Energy per nucleon Z=1",nbins,rbinpt);
    hen2               = new TH1F("hen2","Kin.Energy per nucleon Z=2",nbins,rbinpt);

    hdedxtrkvsdef      = new TH2F("hdedxtrkvsdef","dE/dx (Trk-average) vs deflection",500,-10.,10.,600,0.,40.);
    hbetavsdef         = new TH2F("hbetavsdef","beta vs deflection",500,-10.,10.,200,0.4,1.5);
    hbetavsdef2        = new TH2F("hbetavsdef2","beta vs deflection (Z=1)",500,-10.,10.,200,0.4,1.5);
    hdef               = new TH1F("hdef","Deflection distribution",5000,-2.5,2.5);

    herrdef            = new TH1F("herrdef","Deflection error",1000,0.,0.01);
    
    for(Int_t ip=0; ip<6; ip++){
        hid.Form("herrx%i",ip);
        title = ""; title.Form("Spatial resolution - plane %i (x)",ip);
        herrx[ip] = new TH1F(hid.Data(),title.Data(),500,0.,0.005);
        hid.Form("herry%i",ip);
        title = ""; title.Form("Spatial resolution - plane %i (y)",ip);
        herry[ip] = new TH1F(hid.Data(),title.Data(),500,0.,0.005);
    }

}
//===============================================================
// Retrieve histograms
//===============================================================
void RetrieveHistos(){

    TString hid;

    hchi2def           = dynamic_cast<TH2F*>(gDirectory->FindObject("hchi2def"));
    hchi2rig           = dynamic_cast<TH2F*>(gDirectory->FindObject("hchi2rig"));
    hdedxtofvsdedxtrk  = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtofvsdedxtrk"));
    hdedxtrkvsrig      = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtrkvsrig"));
    hdedxtofvsrig      = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtofvsrig"));
    hbetavsrig         = dynamic_cast<TH2F*>(gDirectory->FindObject("hbetavsrig"));
    hdedxtrkvsdef      = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtrkvsdef"));
    hbetavsdef         = dynamic_cast<TH2F*>(gDirectory->FindObject("hbetavsdef"));
    hbetavsdef2        = dynamic_cast<TH2F*>(gDirectory->FindObject("hbetavsdef2"));
    hdef               = dynamic_cast<TH1F*>(gDirectory->FindObject("hdef"));
    hrig1              = dynamic_cast<TH1F*>(gDirectory->FindObject("hrig1"));
    hrig2              = dynamic_cast<TH1F*>(gDirectory->FindObject("hrig2"));
    hen1               = dynamic_cast<TH1F*>(gDirectory->FindObject("hen1"));
    hen2               = dynamic_cast<TH1F*>(gDirectory->FindObject("hen2"));
    herrdef            = dynamic_cast<TH1F*>(gDirectory->FindObject("herrdef"));
    for(Int_t ip=0; ip<6; ip++){
        hid.Form("herrx%i",ip);
        herrx[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
        hid.Form("herry%i",ip);
        herry[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
    }
}
//===============================================================
// Fill histograms (called event-by-event)
//===============================================================
 
void FillHistos( PamLevel2* event ){
//  ----------------------------------
//  retrieve histos
//  ----------------------------------
    RetrieveHistos();

//  ----------------------------------
//  single-track events
//  ----------------------------------
    float beta = 0;
    float dedxtof = 0;
    float dedxtrk = 0;
    float rigidity = 0;
    float deflection = 0;
    float chi2 = 0;
    float errdef = 0.;
    float errx[] = {10.,10.,10.,10.,10.,10.};
    float erry[] = {10.,10.,10.,10.,10.,10.};
    
    tot++;
        
    if( event->GetTrkLevel2()->GetNTracks()==0 ) return; 
    PamTrack *track = event->GetTrack(0);
    //------------------------------
    // track selection
    //------------------------------
    if( 
        true
        ){
        
        sel[0]++;
        //-------------------
        // Calo variables
        //-------------------
        
        
        //-------------------
        // ToF variables
        //-------------------
        beta    = track->GetToFTrack()->beta[12];
//      dedxtof = GetdEdx(track);
        
        //-------------------
        // tracker variables
        //-------------------
        dedxtrk    = track->GetTrkTrack()->GetDEDX()/1.3; //<<<ATTENZIONE TEMPORANEO!!!!
        rigidity   = track->GetTrkTrack()->GetRigidity();
        deflection = track->GetTrkTrack()->GetDeflection();
        chi2       = track->GetTrkTrack()->chi2;
        errdef     = sqrt(track->GetTrkTrack()->coval[4][4]); 
        for(Int_t ip=0; ip<6; ip++){
            if(track->GetTrkTrack()->XGood(ip))errx[ip]=track->GetTrkTrack()->resx[ip];
            if(track->GetTrkTrack()->YGood(ip))erry[ip]=track->GetTrkTrack()->resy[ip];
        }

        //-------------------
        // chi2 selection
        //-------------------
        hchi2def->Fill( fabs(deflection), pow((double)chi2,0.25) );
        hchi2rig->Fill( rigidity, chi2 );
        if(pow((double)chi2,0.25) > (2.5+1.85*fabs(deflection)) )return;
        sel[1]++;
        
//      hdedxtofvsdedxtrk->Fill(dedxtrk,dedxtof);
        hdedxtrkvsrig->Fill(rigidity,dedxtrk);
//      hdedxtofvsrig->Fill(rigidity,dedxtof);
        hbetavsrig->Fill(rigidity,beta);

        hdedxtrkvsdef->Fill(deflection,dedxtrk);
        hbetavsdef->Fill(deflection,beta);

        //------------------------
        // rough charge selection
        //------------------------
        if(dedxtrk > (6.+16./(rigidity*rigidity)))return;
        //--------
        // helium
        //--------
        if( dedxtrk > (2.8+3.3/(rigidity*rigidity)) && deflection>0 ){
            sel[5]++;
            hrig2->Fill(rigidity);
            float m=3.725;
            float a=4.;
            float z=2.;
            float en=(sqrt( (z*rigidity)*(z*rigidity)+m*m )-m)/a;
            hen2->Fill(en);
        }

        //--------
        // Z=1
        //--------
        if(dedxtrk > (2.8+3.3/(rigidity*rigidity)))return;
        if(dedxtrk < (0.4+0.35/(rigidity*rigidity)))return;
        sel[2]++;

        if(deflection>0){
            hrig1->Fill(rigidity);
            float m=0.938;
            float a=1.;
            float z=1.;
            float en=(sqrt( (z*rigidity)*(z*rigidity)+m*m )-m)/a;
            hen1->Fill(en);
        }
        //------------------------
        // 6-planes lever-arm
        //------------------------
        if( !(track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(5)) )return;
        sel[3]++;

        for(Int_t ip=0; ip<6; ip++)herrx[ip]->Fill(errx[ip]);
        for(Int_t ip=0; ip<6; ip++)herry[ip]->Fill(erry[ip]);

        //---------------------------
        // cut on spatial resolution
        //---------------------------
        if( errx[0] > 0.0004 )return;
        if( errx[5] > 0.0004 )return;
        sel[4]++;

        herrdef->Fill(errdef);
        hbetavsdef2->Fill(deflection,beta);
        hdef->Fill(deflection);
    }
    
//    delete track;
   
    
}
//===============================================================
// Save histograms
//===============================================================
void SaveHistos( TFile * outf ){

    TString hid; 

    gROOT->cd();
//  ----------------------------------
//  retrieve histos
//  ----------------------------------
    RetrieveHistos();

//  ----------------------------------
//  save on file
//  ----------------------------------
    outf->cd();
    
    hchi2def->Write();
    hchi2rig->Write();
    hdedxtofvsdedxtrk->Write(); 
    hdedxtrkvsrig->Write();   
    hdedxtofvsrig->Write(); 
    hbetavsrig->Write();  
    hdedxtrkvsdef->Write(); 
    hbetavsdef->Write(); 
    hbetavsdef2->Write(); 
    hdef->Write(); 
    hrig1->Write(); 
    hrig2->Write(); 
    hen1->Write(); 
    hen2->Write(); 
    herrdef->Write(); 
    for(Int_t ip=0; ip<6; ip++){
        herrx[ip]->Write(); 
        herry[ip]->Write(); 
    }

    cout << "------------------------------------------------------------------" <<endl;
    cout << "Pre-selected events :"<<tot<<endl;
    if(!tot)return;
    cout << "chi2 selection                            :"<<sel[1]<< "  (" << ((float)sel[1]/(float)tot) <<")"<<endl;
    cout << "Z=2 selection                             :"<<sel[5]<< "  (" << ((float)sel[5]/(float)tot) <<")"<<endl;
    cout << "Z=1 selection                             :"<<sel[2]<< "  (" << ((float)sel[2]/(float)tot) <<")"<<endl;
    cout << "(6-plane lever-arm)                       :"<<sel[3]<< "  (" << ((float)sel[3]/(float)tot) <<")"<<endl;
    cout << "ext.spatial res. on planes 1-6 < 4 micron :"<<sel[4]<< "  (" << ((float)sel[4]/(float)tot) <<")"<<endl;
    cout << "------------------------------------------------------------------" <<endl;

}

1	pam-fi	1.1	////////////////////////////////////////////////////////////////
2			//
3			// Histogram set for tracker-efficiency study.
4			//
5			// It contains the functions:
6			//
7			// - CreateHistos()
8	pam-fi	1.2	// - RetrieveHistos()
9	pam-fi	1.1	// - FillHistos(PamLevel2* event)
10			// - Save Histos()
11			//
12			// Usage (example):
13			// [] .L My-histo-set-macro.C+
14			// [] .L My-selection.C+
15			// [] .L Loop.C+
16			// [] Loop("data-dir/","file-list.txt",100000,"+ALL fillTree fillHistos","output-file.root")
17			//
18			////////////////////////////////////////////////////////////////
19
20			#if !defined(__CINT__) \|\| defined(__MAKECINT__)
21
22			#include <TString.h>
23			#include <TH1F.h>
24			#include <TH2F.h>
25			#include <TMath.h>
26			#include <TLine.h>
27			#include <TPolyMarker.h>
28			#include <TSelector.h>
29			#include <TFile.h>
30
31			#include <stdlib.h>
32			#include <iostream>
33			using namespace std;
34
35			#include <PamLevel2.h>
36			#include <CaloAxis.h>
37
38			#endif
39			//======================
40			// HITOGRAMS declaration
41			//======================
42	pam-fi	1.2
43			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	pam-fi	1.3	// 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	pam-fi	1.2
102	pam-fi	1.3	// }
103	pam-fi	1.2
104
105	pam-fi	1.1	//===============================================================
106			// Create histograms
107			//===============================================================
108			void CreateHistos( TFile* outf ){
109
110			gROOT->cd();//create histos in memory
111
112	pam-fi	1.2	TString hid;
113			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","(chi2)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	pam-fi	1.1
158			}
159	pam-fi	1.2	//===============================================================
160			// Retrieve histograms
161			//===============================================================
162			void RetrieveHistos(){
163
164			TString hid;
165	pam-fi	1.1
166	pam-fi	1.2	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	pam-fi	1.1	//===============================================================
189			// Fill histograms (called event-by-event)
190			//===============================================================
191
192			void FillHistos( PamLevel2* event ){
193			// ----------------------------------
194			// retrieve histos
195			// ----------------------------------
196	pam-fi	1.2	RetrieveHistos();
197	pam-fi	1.1
198			// ----------------------------------
199	pam-fi	1.2	// single-track events
200	pam-fi	1.1	// ----------------------------------
201	pam-fi	1.2	float beta = 0;
202			float dedxtof = 0;
203			float dedxtrk = 0;
204			float rigidity = 0;
205			float deflection = 0;
206			float chi2 = 0;
207			float errdef = 0.;
208			float errx[] = {10.,10.,10.,10.,10.,10.};
209			float erry[] = {10.,10.,10.,10.,10.,10.};
210	pam-fi	1.1
211	pam-fi	1.2	tot++;
212
213	pam-fi	1.3	if( event->GetTrkLevel2()->GetNTracks()==0 ) return;
214	pam-fi	1.2	PamTrack *track = event->GetTrack(0);
215			//------------------------------
216			// track selection
217			//------------------------------
218			if(
219			true
220			){
221
222			sel[0]++;
223			//-------------------
224			// Calo variables
225			//-------------------
226
227
228			//-------------------
229			// ToF variables
230			//-------------------
231	pam-fi	1.3	beta = track->GetToFTrack()->beta[12];
232	pam-fi	1.2	// dedxtof = GetdEdx(track);
233
234			//-------------------
235			// tracker variables
236			//-------------------
237	pam-fi	1.3	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	pam-fi	1.2	for(Int_t ip=0; ip<6; ip++){
243	pam-fi	1.3	if(track->GetTrkTrack()->XGood(ip))errx[ip]=track->GetTrkTrack()->resx[ip];
244			if(track->GetTrkTrack()->YGood(ip))erry[ip]=track->GetTrkTrack()->resy[ip];
245	pam-fi	1.2	}
246	pam-fi	1.1
247	pam-fi	1.2	//-------------------
248			// chi2 selection
249			//-------------------
250			hchi2def->Fill( fabs(deflection), pow((double)chi2,0.25) );
251			hchi2rig->Fill( rigidity, chi2 );
252			if(pow((double)chi2,0.25) > (2.5+1.85*fabs(deflection)) )return;
253			sel[1]++;
254
255			// hdedxtofvsdedxtrk->Fill(dedxtrk,dedxtof);
256			hdedxtrkvsrig->Fill(rigidity,dedxtrk);
257			// hdedxtofvsrig->Fill(rigidity,dedxtof);
258			hbetavsrig->Fill(rigidity,beta);
259
260			hdedxtrkvsdef->Fill(deflection,dedxtrk);
261			hbetavsdef->Fill(deflection,beta);
262
263			//------------------------
264			// rough charge selection
265			//------------------------
266			if(dedxtrk > (6.+16./(rigidity*rigidity)))return;
267			//--------
268			// helium
269			//--------
270			if( dedxtrk > (2.8+3.3/(rigidity*rigidity)) && deflection>0 ){
271			sel[5]++;
272			hrig2->Fill(rigidity);
273			float m=3.725;
274			float a=4.;
275			float z=2.;
276			float en=(sqrt( (zrigidity)(zrigidity)+mm )-m)/a;
277			hen2->Fill(en);
278			}
279	pam-fi	1.1
280	pam-fi	1.2	//--------
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( (zrigidity)(zrigidity)+mm )-m)/a;
293			hen1->Fill(en);
294	pam-fi	1.1	}
295	pam-fi	1.2	//------------------------
296			// 6-planes lever-arm
297			//------------------------
298	pam-fi	1.3	if( !(track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(5)) )return;
299	pam-fi	1.2	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	pam-fi	1.1	}
315
316	pam-fi	1.3	// delete track;
317	pam-fi	1.2
318
319	pam-fi	1.1	}
320			//===============================================================
321			// Save histograms
322			//===============================================================
323			void SaveHistos( TFile * outf ){
324
325	pam-fi	1.2	TString hid;
326
327	pam-fi	1.1	gROOT->cd();
328			// ----------------------------------
329			// retrieve histos
330			// ----------------------------------
331	pam-fi	1.2	RetrieveHistos();
332	pam-fi	1.1
333			// ----------------------------------
334			// save on file
335			// ----------------------------------
336			outf->cd();
337
338	pam-fi	1.2	hchi2def->Write();
339			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	pam-fi	1.1	}
369