--- PamelaLevel2/doc/examples/My-Histos.cpp 2006/12/11 18:26:55 1.1 +++ PamelaLevel2/doc/examples/My-Histos.cpp 2007/01/03 13:28:49 1.2 @@ -5,6 +5,7 @@ // It contains the functions: // // - CreateHistos() +// - RetrieveHistos() // - FillHistos(PamLevel2* event) // - Save Histos() // @@ -38,8 +39,68 @@ //====================== // HITOGRAMS declaration //====================== -TH1F *hresangx; -TH1F *hresangy; + +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; inpmtadc; i++){ + if(trk->dedx[i] == 0)ISNULL=true; + dedx += (trk->dedx).At(i)*(Int_t)(!ISNULL); + nval += (Int_t)(!ISNULL); + }; + // + if(nval)dedx=dedx/nval; +// cout <<" -- dedx "<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 //=============================================================== @@ -47,11 +108,82 @@ gROOT->cd();//create histos in memory - hresangx = new TH1F("hresangx","Angular residual on the top calo plane (x)",400,-40.,40.); - hresangy = new TH1F("hresangy","Angular residual on the top calo plane (y)",400,-40.,40.); + 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(gDirectory->FindObject("hchi2def")); + hchi2rig = dynamic_cast(gDirectory->FindObject("hchi2rig")); + hdedxtofvsdedxtrk = dynamic_cast(gDirectory->FindObject("hdedxtofvsdedxtrk")); + hdedxtrkvsrig = dynamic_cast(gDirectory->FindObject("hdedxtrkvsrig")); + hdedxtofvsrig = dynamic_cast(gDirectory->FindObject("hdedxtofvsrig")); + hbetavsrig = dynamic_cast(gDirectory->FindObject("hbetavsrig")); + hdedxtrkvsdef = dynamic_cast(gDirectory->FindObject("hdedxtrkvsdef")); + hbetavsdef = dynamic_cast(gDirectory->FindObject("hbetavsdef")); + hbetavsdef2 = dynamic_cast(gDirectory->FindObject("hbetavsdef2")); + hdef = dynamic_cast(gDirectory->FindObject("hdef")); + hrig1 = dynamic_cast(gDirectory->FindObject("hrig1")); + hrig2 = dynamic_cast(gDirectory->FindObject("hrig2")); + hen1 = dynamic_cast(gDirectory->FindObject("hen1")); + hen2 = dynamic_cast(gDirectory->FindObject("hen2")); + herrdef = dynamic_cast(gDirectory->FindObject("herrdef")); + for(Int_t ip=0; ip<6; ip++){ + hid.Form("herrx%i",ip); + herrx[ip] = dynamic_cast(gDirectory->FindObject(hid.Data())); + hid.Form("herry%i",ip); + herry[ip] = dynamic_cast(gDirectory->FindObject(hid.Data())); + } +} //=============================================================== // Fill histograms (called event-by-event) //=============================================================== @@ -60,103 +192,177 @@ // ---------------------------------- // retrieve histos // ---------------------------------- -// hresangx = dynamic_cast(gDirectory->FindObject("hresangx")); -// hresangy = dynamic_cast(gDirectory->FindObject("hresangy")); + RetrieveHistos(); // ---------------------------------- -// calo variables +// single-track events // ---------------------------------- - CaloAxis *x_axis = new CaloAxis(); - CaloAxis *y_axis = new CaloAxis(); - - float rcil = 1.;// tolerance (cm) - x_axis->FitAxis(event->GetCaloLevel1(),0,rcil); - y_axis->FitAxis(event->GetCaloLevel1(),1,rcil); - - float qtrack = x_axis->GetQaxis()+y_axis->GetQaxis(); - int ntrack = x_axis->GetN()+y_axis->GetN(); + 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.}; - float caloangx = (float)atan((double)x_axis->par[1])*180./3.1415026; - float caloangy = (float)atan((double)y_axis->par[1])*180./3.1415026; - - CaloStrip top = CaloStrip(); - top.Set(0,0,48); - float zcalotop = top.GetZ(); - -// ---------------------------------- -// tof variables -// ---------------------------------- -// get track stored by the tof - ToFTrkVar *tof = event->GetToFStoredTrack(-1); - if(!tof){ - cout << " no ToF-track stored "<GetNTracks()==0 ) return; + PamTrack *track = event->GetTrack(0); + //------------------------------ + // track selection + //------------------------------ + if( + true + ){ + + sel[0]++; + //------------------- + // Calo variables + //------------------- + + + //------------------- + // ToF variables + //------------------- + beta = track->beta[12]; +// dedxtof = GetdEdx(track); + + //------------------- + // tracker variables + //------------------- + dedxtrk = track->GetDEDX()/1.3; //<<GetRigidity(); + deflection = track->GetDeflection(); + chi2 = track->chi2; + errdef = sqrt(track->coval[4][4]); + for(Int_t ip=0; ip<6; ip++){ + if(track->XGood(ip))errx[ip]=track->resx[ip]; + if(track->YGood(ip))erry[ip]=track->resy[ip]; + } - float resangx = -999; - float resangy = -999; - if( event->GetNTracks()==1 ){ - - PamTrack *track = event->GetTrack(0); - if( track->chi2 >0 ){ - - Float_t ztraj[13]; - - Int_t i=0; - ztraj[i++] = event->GetZTOF(11); - ztraj[i++] = event->GetZTOF(12); - ztraj[i++] = event->GetZTOF(21); - ztraj[i++] = event->GetZTOF(22); - ztraj[i++] = event->GetZTrk(1); - ztraj[i++] = event->GetZTrk(2); - ztraj[i++] = event->GetZTrk(3); - ztraj[i++] = event->GetZTrk(4); - ztraj[i++] = event->GetZTrk(5); - ztraj[i++] = event->GetZTrk(6); - ztraj[i++] = event->GetZTOF(31); - ztraj[i++] = event->GetZTOF(32); - ztraj[i++] = zcalotop; - Trajectory traj = Trajectory(13,ztraj); - traj.DoTrack2( track->al ); - - resangx = caloangx - traj.thx[12]; - resangy = caloangy - traj.thy[12]; + //------------------- + // 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); } - if(track) delete track; + //------------------------ + // 6-planes lever-arm + //------------------------ + if( !(track->XGood(0) && track->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); } -// ---------------------------------- -// fill histograms -// ---------------------------------- - hresangx->Fill(resangx); - hresangy->Fill(resangy); - - + delete track; + + } //=============================================================== // Save histograms //=============================================================== void SaveHistos( TFile * outf ){ + TString hid; + gROOT->cd(); // ---------------------------------- // retrieve histos // ---------------------------------- -// hresangx = dynamic_cast(gDirectory->FindObject("hresangx")); -// hresangy = dynamic_cast(gDirectory->FindObject("hresangy")); + RetrieveHistos(); // ---------------------------------- // save on file // ---------------------------------- - outf->cd(); - if(hresangx) hresangx->Write(); - if(hresangy) hresangy->Write(); + 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 << "------------------------------------------------------------------" <