Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-03 08:13:10

0001 #include <TH2D.h>
0002 #include <TFile.h>
0003 #include <TF1.h>
0004 #include <TGraph.h>
0005 #include  <TCanvas.h>
0006 #include  <TColor.h>
0007 #include <TROOT.h>
0008 #include <TStyle.h>
0009 #include <TLegend.h>
0010 #include <TLatex.h>
0011 #include <TLine.h>
0012 #include <TEfficiency.h>
0013 #include <TGraphAsymmErrors.h>
0014 
0015 #include "sPhenixStyle.C"
0016 
0017 void look_purity()
0018 {
0019   SetsPhenixStyle();
0020   /*
0021   gROOT->SetStyle("Plain");
0022   gStyle->SetOptStat(0);
0023   gStyle->SetOptFit(0);
0024   gStyle->SetOptTitle(0);
0025   */
0026 
0027   //gStyle->SetOptStat(1);
0028 
0029   // determines if plots are labeled "HIJING"
0030   bool hijing = false;
0031 
0032   //int n_maps_layer = 3;
0033   int n_maps_layer = 0;
0034   //int n_intt_layer = 4;
0035   int n_intt_layer = 0;
0036   int n_tpc_layer = 48;
0037     
0038   double ptmax = 40.0;
0039   double hptmax = 40.0;
0040   //double ptstep = 0.25;
0041   //double ptstep = 1.0;
0042   //double ptstep = 2.0;
0043   //double ptstep = 0.5;
0044   double ptstep = 0.4;
0045 
0046   // set to false only to generate pT resolution plots without fits
0047   // BEWARE: false means that the 4 sigma cuts are meaningless - they are not done with fitted parameters
0048   bool pt_resolution_fit = true;
0049   //bool pt_resolution_fit = false;
0050   
0051   TFile *fin = new TFile("root_files/purity_out.root"); 
0052   //TFile *fin = new TFile("root_files/purity_out_100pions_noINTT_primvert.root"); 
0053 
0054 
0055   if(!fin)
0056     {
0057       cout << "Failed to find input file" << endl;
0058       exit(1);
0059     }
0060 
0061   TCanvas *c1 = new TCanvas("c1","c1",5,5,900,500);
0062   c1->Divide(3,1);
0063 
0064   //==========================
0065   // 2D histo for embedded pions
0066   TH2D *hpt_dca2d = 0;
0067   fin->GetObject("hpt_dca2d",hpt_dca2d);
0068   c1->cd(1);
0069   gPad->SetLeftMargin(0.12);
0070   hpt_dca2d->GetYaxis()->SetTitleOffset(1.5);
0071   hpt_dca2d->SetMarkerStyle(20);
0072   hpt_dca2d->SetMarkerSize(0.3);
0073   hpt_dca2d->Draw();
0074 
0075   // 2D histo for embedded pions
0076   TH2D *hpt_dcaZ = 0;
0077   fin->GetObject("hpt_dcaZ",hpt_dcaZ);
0078   c1->cd(2);
0079   gPad->SetLeftMargin(0.12);
0080   hpt_dcaZ->SetMarkerStyle(20);
0081   hpt_dcaZ->SetMarkerSize(0.3);
0082   hpt_dcaZ->GetYaxis()->SetTitleOffset(1.5);
0083   hpt_dcaZ->Draw();
0084 
0085   TH2D *hpt_compare = 0;
0086   fin->GetObject("hpt_compare",hpt_compare);
0087   c1->cd(3);
0088   gPad->SetLeftMargin(0.12);
0089   hpt_compare->GetYaxis()->SetTitleOffset(1.5);
0090   hpt_compare->SetMarkerStyle(20);
0091   hpt_compare->SetMarkerSize(0.3);
0092   hpt_compare->Draw();
0093   
0094   //========================================
0095   // extract DCA resolution vs pT from 2D histo hpt_dca2d
0096   // by fitting NPT pT slices
0097 
0098   TCanvas *c2 = new TCanvas("c2","c2",20,20,800,600);
0099 
0100   int NPT = (int) (ptmax/ptstep);
0101 
0102   double pT[200];
0103   double dca2d[200];
0104   for(int i = 0;i<NPT;i++)
0105     {
0106       // divide pT range into bins of width ptstep GeV/c at 0.25, 0.75, 1.25, .....
0107       // bins will be (e.g.)from 0-0.5, 0.5-1.0, 1.0-1.5, .....
0108       double ptlo = (double) i * ptstep + 0.0;
0109       double pthi = (double) i * ptstep + ptstep;
0110 
0111       int binlo = hpt_dca2d->GetXaxis()->FindBin(ptlo);
0112       int binhi = hpt_dca2d->GetXaxis()->FindBin(pthi);
0113 
0114       std::cout << "ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
0115 
0116       TH1D *h = new TH1D("h","dca2d resolution",2000, -0.1, 0.1);
0117       hpt_dca2d->ProjectionY("h",binlo,binhi);
0118       h->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0119       h->GetXaxis()->SetTitle("#Delta dca2d (cm)");
0120       h->GetXaxis()->SetTitleOffset(1.0);
0121       //if(i<8) h->Rebin(4);
0122       h->DrawCopy();
0123 
0124       double mean = h->GetMean();
0125       double sigma = h->GetRMS();
0126       double yield = h->Integral();
0127   
0128      pT[i] = (ptlo + pthi) / 2.0;
0129     
0130       double low = -0.01, high=0.01;
0131       if(n_maps_layer == 0)
0132     {
0133       low = 3.0 * low;
0134       high = 3.0 * high;
0135     } 
0136       if(pT[i] < 6.0)
0137     {
0138       low = 3.0*low;
0139       high = 3.0*high;
0140     }
0141 
0142       TF1 *f = new TF1("f","gaus",low, high);
0143       f->SetParameter(1, yield/100.0);
0144       f->SetParameter(2, 0.0);
0145       f->SetParameter(3,0.002);
0146       h->Fit(f,"R");
0147 
0148  
0149 
0150       dca2d[i] = f->GetParameter(2);
0151       cout << " pT " << pT[i] << " dca2d " << dca2d[i] << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
0152     }
0153 
0154   //============================================
0155   // plot the dca2d resolution extracted above for embedded pions
0156 
0157   TCanvas *c3 = new TCanvas("c3","c3",100,100,800,600);
0158   TGraph *grdca2d = new TGraph(NPT,pT,dca2d);
0159   grdca2d->SetMarkerStyle(20);
0160   grdca2d->SetMarkerSize(1.2);
0161   grdca2d->SetMarkerColor(kRed);
0162   grdca2d->SetName("dca2d_resolution");
0163   grdca2d->SetTitle("dca2d resolution");
0164 
0165   TH1D *hdummy = new TH1D("hdummy","#Delta dca2d vs p_{T}",100,0.0,hptmax);
0166   hdummy->SetMinimum(0);
0167   if(n_maps_layer == 0)
0168     hdummy->SetMaximum(0.015);
0169   else
0170     hdummy->SetMaximum(0.0051);
0171   hdummy->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0172   hdummy->GetYaxis()->SetTitle("DCA(r#phi) resolution (cm)");
0173   hdummy->GetYaxis()->SetTitleOffset(1.5);
0174   gPad->SetLeftMargin(0.15);
0175   hdummy->Draw();
0176   grdca2d->Draw("p");
0177 
0178   char dcalab[500];
0179   if(hijing)
0180     //sprintf(dcalab, "sPHENIX   HIJING   b=0-4 fm");
0181     sprintf(dcalab, "sPHENIX   HIJING   b=0-12 fm, 100 kHz pileup");
0182   else
0183   sprintf(dcalab, "sPHENIX  100 pions");
0184   TLegend *ldca = new TLegend(0.3, 0.75, 1.05, 0.90,dcalab,"NDC");
0185   ldca->SetBorderSize(0);
0186   ldca->SetFillColor(0);
0187   ldca->SetFillStyle(0);
0188   char lstr1[500];
0189   //sprintf(lstr1,"MVTX+INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0190   sprintf(lstr1,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0191   ldca->AddEntry(grdca2d, lstr1, "p");
0192   ldca->Draw();
0193 
0194   //===================================
0195 
0196   //========================================
0197   // extract DCA Z resolution vs pT from 2D histo hpt_dcaZ
0198   // by fitting NPT pT slices
0199 
0200   TCanvas *c2a = new TCanvas("c2a","c2a",20,20,800,600);
0201 
0202   double dcaZ[200];
0203   for(int i = 0;i<NPT;i++)
0204     {
0205       // divide pT range into bins of width 0.5 GeV/c at 0.25, 0.75, 1.25, .....
0206       // bins will be from 0-0.5, 0.5-1.0, 1.0-1.5, .....
0207       double ptlo = (double) i * ptstep + 0.0;
0208       double pthi = (double) i * ptstep + ptstep;
0209 
0210       int binlo = hpt_dcaZ->GetXaxis()->FindBin(ptlo);
0211       int binhi = hpt_dcaZ->GetXaxis()->FindBin(pthi);
0212 
0213       std::cout << "ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
0214       TH1D *h = new TH1D("h","DCA (Z) resolution (cm)",2000, -0.1, 0.1);
0215       hpt_dcaZ->ProjectionY("h",binlo,binhi);
0216       h->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0217       h->GetYaxis()->SetTitle("DCA(Z) resolution (cm)");
0218       h->GetYaxis()->SetTitleOffset(1.0);
0219       //if(i<8) h->Rebin(4);
0220       h->DrawCopy();
0221 
0222       double mean = h->GetMean();
0223       double sigma = h->GetRMS();
0224       double yield = h->Integral();
0225   
0226      pT[i] = (ptlo + pthi) / 2.0;
0227     
0228       double low = -0.01, high=0.01;
0229       if(n_maps_layer == 0)
0230     {
0231       low = 3.0 * low;
0232       high = 3.0 * high;
0233     }
0234   
0235       if(pT[i] < 6.0)
0236     {
0237       low = 3.0*low;
0238       high = 3.0*high;
0239     }
0240 
0241       TF1 *f = new TF1("f","gaus",low, high);
0242       f->SetParameter(1, yield/100.0);
0243       f->SetParameter(2, 0.0);
0244       f->SetParameter(3,0.002);
0245       h->Fit(f,"R"); 
0246 
0247       dcaZ[i] = f->GetParameter(2);
0248       cout << " pT " << pT[i] << " dcaZ " << dcaZ[i] << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
0249     }
0250 
0251   //============================================
0252   // plot the dcaZ resolution extracted above for embedded pions
0253 
0254   TCanvas *c3a = new TCanvas("c3a","c3a",100,100,800,600);
0255   TGraph *grdcaZ = new TGraph(NPT,pT,dcaZ);
0256   grdcaZ->SetMarkerStyle(20);
0257   grdcaZ->SetMarkerSize(1.2);
0258   grdcaZ->SetMarkerColor(kRed);
0259   grdcaZ->SetName("dcaZ_resolution");
0260   grdcaZ->SetTitle("dcaZ resolution");
0261 
0262   TH1D *hdummya = new TH1D("hdummya","#Delta dcaZ vs p_{T}",100,0.0,hptmax);
0263   hdummya->SetMinimum(0);
0264 
0265   if(n_maps_layer == 0)  
0266     hdummya->SetMaximum(0.051);
0267   else
0268     hdummya->SetMaximum(0.0051);
0269   hdummya->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0270   hdummya->GetYaxis()->SetTitle("DCA(Z) resolution  (cm)");
0271   hdummya->GetYaxis()->SetTitleOffset(1.5);
0272   gPad->SetLeftMargin(0.15);
0273   hdummya->Draw();
0274   grdcaZ->Draw("p");
0275 
0276   TLegend *ldcaZ = new TLegend(0.3, 0.75, 1.05, 0.90, dcalab,"NDC");
0277   ldcaZ->SetBorderSize(0);
0278   ldcaZ->SetFillColor(0);
0279   ldcaZ->SetFillStyle(0);
0280   //char lstr1[500];
0281   //sprintf(lstr1,"MVTX+INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0282   sprintf(lstr1,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0283   ldcaZ->AddEntry(grdcaZ, lstr1, "p");
0284   ldcaZ->Draw();
0285 
0286   //===================================
0287 
0288 
0289   // extract pT resolution vs pT from hpt_compare
0290 
0291   TCanvas *c4 = new TCanvas("c4","c4",60,60,800,600);
0292 
0293   double dpT[200];
0294 
0295   dpT[0] = 1.0; // throw it away, stats are too poor
0296   for(int i = 1;i<NPT;i++)
0297     //for(int i = 0;i<1;i++)
0298     {
0299       // divide pT range into bins of width 0.5 GeV/c at 0.25, 0.75, 1.25, .....
0300       // bins will be from 0-0.5, 0.5-1.0, 1.0-1.5, .....
0301       double ptlo = (double) i * ptstep + 0.0;
0302       double pthi = (double) i * ptstep + ptstep;
0303 
0304       int binlo = hpt_compare->GetXaxis()->FindBin(ptlo);
0305       int binhi = hpt_compare->GetXaxis()->FindBin(pthi);
0306 
0307       TH1D *hpt = new TH1D("hpt","pT resolution ",200, 0, 2.0);    
0308       hpt_compare->ProjectionY("hpt",binlo,binhi);
0309       hpt->GetXaxis()->SetTitle("#Delta p_{T}/p_{T}");
0310       hpt->GetXaxis()->SetTitleOffset(1.0);
0311       if(i>30) hpt->Rebin(4);
0312       hpt->DrawCopy();
0313 
0314       std::cout << "ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << " integral " << hpt->Integral() << std::endl;
0315 
0316       TF1 *f = new TF1("f","gaus",0.8,1.2);
0317       hpt->Fit(f,"R");
0318 
0319       pT[i] = (ptlo + pthi) / 2.0;
0320 
0321       dpT[i] = f->GetParameter(2);
0322       cout << " pT " << pT[i] << " dpT " << dpT[i] << " integral " << hpt->Integral() << std::endl;
0323     }
0324 
0325   //==========================================
0326   // Plot pT resolution extracted above for embedded pions
0327 
0328   TCanvas *c5 = new TCanvas("c5","c5",100,100,800,800);
0329   c5->SetLeftMargin(0.14);
0330   TGraph *grdpt = new TGraph(NPT,pT,dpT);
0331   grdpt->SetMarkerStyle(20);
0332   grdpt->SetMarkerSize(1.1);
0333   grdpt->SetName("pt_resolution");
0334   grdpt->SetTitle("pT resolution");
0335 
0336   TH1D *hdummy2 = new TH1D("hdummy2","#Delta p_{T} vs p_{T}",100,0.0,hptmax);
0337   hdummy2->SetMinimum(0);
0338   //hdummy2->SetMaximum(0.05);
0339   //hdummy2->SetMaximum(0.12);
0340   hdummy2->SetMaximum(0.3);
0341   hdummy2->GetXaxis()->SetTitle("p_{T}");
0342   hdummy2->GetYaxis()->SetTitle("#Delta p_{T}/p_{T}");
0343   hdummy2->GetYaxis()->SetTitleOffset(1.5);
0344   hdummy2->Draw();
0345   grdpt->Draw("p");
0346 
0347  TLegend *ldpt= new TLegend(0.25, 0.75, 0.95, 0.90, dcalab, "NDC");
0348   ldpt->SetBorderSize(0);
0349   ldpt->SetFillColor(0);
0350   ldpt->SetFillStyle(0);
0351   //char lstr1[500];
0352   //sprintf(lstr1,"MVTX+INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0353   sprintf(lstr1,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0354   ldpt->AddEntry(grdpt, lstr1, "p");
0355   ldpt->Draw();
0356 
0357   // Parameterize pT resolution
0358   
0359   //TF1 *fpt = new TF1("fpt","sqrt([0]*[0] + [1]*[1]*x*x)", 0, 35.0);
0360   TF1 *fpt = new TF1("fpt","[0] + [1]*x", 2.0, hptmax);
0361   fpt->SetParameter(0,0.005);
0362   //fpt->FixParameter(0,0.005);
0363   fpt->SetParameter(1,0.0015);
0364   if(pt_resolution_fit)  
0365     grdpt->Fit(fpt,"R");
0366 
0367   char lab[1000];
0368   //sprintf(lab,"#frac{#Deltap_{T}}{p_{T}} = #sqrt{%.4f^{2} + (%.6f #times p_{T})^{2}}", fpt->GetParameter(0), fpt->GetParameter(1))
0369   sprintf(lab,"#frac{#Deltap_{T}}{p_{T}} = %.4f + %.6f #times p_{T}", fpt->GetParameter(0), fpt->GetParameter(1));
0370   TLatex *mres = new TLatex(0.45,0.25,lab);
0371   //mres->SetTextSize(0.1);
0372   mres->SetNDC();
0373   if(pt_resolution_fit)  
0374     mres->Draw();
0375 
0376   // For making a cut on the momentum difference between rgpT and rpT 
0377   double pT_sigmas = 6.0;
0378   double const_term =  fpt->GetParameter(0);
0379   double linear_term = fpt->GetParameter(1);
0380 
0381   TH1D *hpt_truth;
0382   fin->GetObject("hpt_truth",hpt_truth);
0383   if(!hpt_truth)
0384     {
0385       cout << "Failed to get hpt_truth, quit!" << endl;
0386       exit(1);
0387     }
0388   
0389   //TCanvas *ctruth = new TCanvas("ctruth","ctruth", 5,5,800,600);
0390 
0391   TH1D *hpt_matched = new TH1D("hpt_matched","hpt_matched", 500, 0.0, hptmax);
0392   double eff_pt[200];
0393 
0394   for(int i = 0;i<NPT;i++)
0395     {
0396       // divide pT range into bins of width 0.5 GeV/c at 0.25, 0.75, 1.25, .....
0397       // bins will be from 0-0.5, 0.5-1.0, 1.0-1.5, .....
0398       double ptlo = (double) i * ptstep + 0.0;
0399       double pthi = (double) i * ptstep + ptstep;
0400       double ptval = (ptlo+pthi)/2.0;
0401 
0402       int binlo = hpt_dca2d->GetXaxis()->FindBin(ptlo);
0403       int binhi = hpt_dca2d->GetXaxis()->FindBin(pthi);
0404 
0405       TH1D *hpt1 = new TH1D("hpt1","pT resolution ",2000, 0.0, 2.0);    
0406       hpt_compare->ProjectionY("hpt1",binlo,binhi);
0407 
0408       double ptres = sqrt(pow(const_term,2) + pow(linear_term * ptval,2)); 
0409       double ptreslo = 1.0 - ptres * pT_sigmas;
0410       double ptreshi = 1.0 + ptres * pT_sigmas;
0411 
0412       int momlo = hpt1->GetXaxis()->FindBin(ptreslo);
0413       int momhi = hpt1->GetXaxis()->FindBin(ptreshi);
0414 
0415       hpt_matched->Fill(ptval, hpt1->Integral(momlo, momhi));
0416 
0417       int tlo = hpt_truth->FindBin(ptlo);
0418       int thi = hpt_truth->FindBin(pthi);
0419       double truth_yield = hpt_truth->Integral(tlo, thi);;
0420       eff_pt[i] = hpt1->Integral(momlo, momhi) / truth_yield;
0421       //eff_pt[i] = hpt1->Integral() / truth_yield;
0422       cout << " pT " << ptval << " ptres " << ptres << " ptreslo " << ptreslo << " ptreshi " << ptreshi << " momlo " << momlo << " momhi " << momhi << endl;
0423       cout << "      truth_yield " << truth_yield << " yield " << hpt1->Integral(momlo,momhi) << " eff_pt " << eff_pt[i] << endl;
0424 
0425       /*
0426       if(i == 20)
0427     {
0428       ctruth->cd(0);
0429       hpt1->DrawCopy();
0430       TF1 *fres = new TF1("fres","gaus");
0431       fres->SetParameter(0, 20.0);
0432       fres->FixParameter(1, 1.0);
0433       fres->FixParameter(2, ptres);
0434       fres->SetLineColor(kRed);
0435       hpt1->Fit(fres);
0436       fres->DrawCopy("same");
0437       ctruth->Update();
0438       cout << hpt1->Integral() << endl;
0439       cout << ptres << endl;
0440       int k = 0;
0441       int y=0;
0442       cin >> k >> y;
0443     }
0444       */
0445 
0446       cout << " ptval " << ptval 
0447        << " ptreslo " << ptreslo
0448        << " ptreshi " << ptreshi
0449        << " total " << hpt1->Integral()
0450        << " momlo " << momlo
0451        << " momhi " << momhi
0452        << "  cut " << hpt1->Integral(momlo,momhi)
0453        << " tlo  " << tlo
0454        << " thi " << thi
0455        << " truth " << truth_yield
0456        << " eff_pt " << eff_pt[i]
0457        << endl;
0458 
0459       delete hpt1;
0460 
0461     }  
0462   
0463   cout << " create canvas c7" << endl;
0464   TCanvas *c7 = new TCanvas("c7","c7",60,60,800,800);
0465 
0466   TH1F *hd = new TH1F("hd","hd",100, 0.0, hptmax);
0467   hd->SetMinimum(0.0);
0468   hd->SetMaximum(1.0);
0469   hd->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0470   hd->GetYaxis()->SetTitle("Single track efficiency");
0471   hd->GetYaxis()->SetTitleOffset(1.0);
0472   hd->Draw();
0473 
0474   TGraph *gr_eff = new TGraph(NPT,pT,eff_pt);
0475   gr_eff->SetName("single_track_efficiency");
0476   gr_eff->SetMarkerStyle(20);
0477   gr_eff->SetMarkerSize(1);
0478   gr_eff->SetMarkerColor(kRed);
0479 
0480   gr_eff->Draw("p");
0481 
0482  TLegend *leff = new TLegend(0.40, 0.35, 0.99, 0.50, dcalab, "NDC");
0483   leff->SetBorderSize(0);
0484   leff->SetFillColor(0);
0485   leff->SetFillStyle(0);
0486   char lstr[500];
0487   //sprintf(lstr,"MVTX+INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0488   sprintf(lstr,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
0489   leff->AddEntry(gr_eff, lstr, "p");
0490   //leff->AddEntry(lmax, "max possible efficiency", "l");
0491   leff->Draw();
0492 
0493   /*
0494   //=======================
0495   // Get track purity for Hijing events by 
0496   // finding tracks within pt_sigmas of the 
0497   // truth pT
0498 
0499   TH2D *hpt_hijing_compare = 0;
0500   fin->GetObject("hpt_hijing_compare",hpt_hijing_compare);
0501   if(!hpt_hijing_compare)
0502     {
0503       cout << "Did not get hpt_hijing_compare - quit!" << endl;
0504       exit(1);
0505     }
0506 
0507   static const int NVARBINS = 36;
0508   double xbins[NVARBINS+1] = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,
0509                   1.1, 1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
0510                   2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
0511                   4.5,5.0,5.5,6.0,7.0,8.0,10.0};
0512   
0513   TH1D *hpt_hijing_allreco = new TH1D("hpt_hijing_allreco","hpt_hijing_allreco",NVARBINS,xbins);
0514   TH1D *hpt_hijing_matched = new TH1D("hpt_hijing_matched","hpt_hijing_matched",NVARBINS,xbins);
0515   hpt_hijing_matched->GetYaxis()->SetTitle("Purity");
0516 
0517   //double purity_pt[NVARBINS];
0518   
0519   for (int i=0;i<NVARBINS;i++)  
0520     {
0521       double ptlo = xbins[i];
0522       double pthi = xbins[i+1];
0523       double ptval = (ptlo+pthi)/2.0;
0524 
0525       int binlo = hpt_hijing_compare->GetXaxis()->FindBin(ptlo);
0526       int binhi = hpt_hijing_compare->GetXaxis()->FindBin(pthi);
0527 
0528       TH1D *hpt1 = new TH1D("hpt1","pT resolution ",2000, 0.0, 2.0);    
0529       hpt_hijing_compare->ProjectionY("hpt1",binlo,binhi);
0530 
0531       double ptres = sqrt(pow(const_term,2) + pow(linear_term * ptval,2)); 
0532       double ptreslo = 1.0 - ptres * pT_sigmas;
0533       double ptreshi = 1.0 + ptres * pT_sigmas;
0534 
0535       int momlo = hpt1->GetXaxis()->FindBin(ptreslo);
0536       int momhi = hpt1->GetXaxis()->FindBin(ptreshi);
0537 
0538       // This to avoid TEfficiency's stupid insistence that the histograms cannot be filled with weights
0539       for (int j=0;j<(int)hpt1->Integral(momlo,momhi);j++)
0540     hpt_hijing_matched->Fill(ptval);
0541       for (int j=0;j<(int)hpt1->Integral();j++)
0542     hpt_hijing_allreco->Fill(ptval);
0543 
0544       //purity_pt[i] = hpt1->Integral(momlo, momhi) / hpt1->Integral();
0545 
0546       delete hpt1;
0547     }
0548 
0549   cout << " create canvas cpur " << endl;
0550   TCanvas *cpur = new TCanvas("cpur","cpur",60,60,1000,600);
0551 
0552   TEfficiency* pEff = 0;
0553   //if(TEfficiency::CheckConsistency(*hpt_hijing_matched,*hpt_hijing_allreco))
0554     {
0555       pEff = new TEfficiency(*hpt_hijing_matched,*hpt_hijing_allreco);
0556       char tname[500];
0557       sprintf(tname,"Reconstruction efficiency (%.1f#sigma p_{T}) ; p_{T} (GeV/c) ; reco'd tracks within %.0f #sigma p_{T}",pT_sigmas,pT_sigmas);
0558       //pEff->SetTitle("Reconstruction efficiency (3.0#sigma p_{T}) ; p_{T} (GeV/c) ; reco'd tracks within 3 #sigma p_{T}");
0559       pEff->SetTitle(tname);
0560       pEff->SetMarkerStyle(20);
0561       pEff->SetMarkerColor(kBlack);
0562       pEff->SetMarkerSize(1);
0563       pEff->Draw();
0564       gPad->Update();
0565       pEff->GetPaintedGraph()->GetYaxis()->SetTitleOffset(1.0);
0566       pEff->GetPaintedGraph()->GetYaxis()->SetRangeUser(0.0,1.1);
0567       pEff->GetPaintedGraph()->GetXaxis()->SetTitleOffset(1.1);
0568     }
0569   */
0570 
0571   //=========================
0572   // plot quality for Hijing tracks
0573 
0574   TCanvas *cpq = new TCanvas("cpq","cpq",40,40,1200,600);
0575   TH1D * hquality;
0576   fin->GetObject("hquality",hquality);
0577   hquality->Draw();
0578 
0579   //======================
0580   // Plot number of hits per track
0581   TCanvas *c8 = new TCanvas("c8","c8",40,40,1200,600);
0582   TH1D * hnhits;
0583   fin->GetObject("hnhits",hnhits);
0584   hnhits->GetXaxis()->SetTitle("hits per track");
0585   hnhits->Draw();
0586   cout << "hnhits integral = " << hnhits->Integral() << endl;
0587 
0588 
0589   //================================
0590   // extract three DCA 2d distributions inside pT bins
0591 
0592   TCanvas *cdcadist = new TCanvas("cdcadist","cdcadist",5,20,1200,800);
0593   cdcadist->Divide(3,1);
0594 
0595   double dcaptbinlo[3] = {0.5, 1.0, 2.0};
0596   double dcaptbinhi[3] = {1.0, 2.0, 50.0};
0597   for(int i = 0;i<3;i++)
0598     {
0599       double ptlo = dcaptbinlo[i];
0600       double pthi = dcaptbinhi[i];
0601 
0602       int binlo = hpt_dca2d->GetXaxis()->FindBin(ptlo);
0603       int binhi = hpt_dca2d->GetXaxis()->FindBin(pthi);
0604 
0605       std::cout << "DCA distr: ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
0606 
0607       TH1D *h = new TH1D("h","dca2d resolution",2000, -0.1, 0.1);
0608       hpt_dca2d->ProjectionY("h",binlo,binhi);
0609       h->GetXaxis()->SetTitle("DCA (r#phi) (cm)");
0610       h->GetXaxis()->SetTitleOffset(1.0);
0611       h->SetMinimum(0.5);
0612 
0613       cdcadist->cd(i+1);
0614       gPad->SetLogy(1);
0615       h->Sumw2();
0616       h->Rebin(4);
0617       //h->SetMarkerStyle(20);
0618       h->SetMarkerSize(0.8);
0619       h->GetXaxis()->SetNdivisions(505);
0620       h->DrawCopy("E");
0621 
0622       double mean = h->GetMean();
0623       double sigma = h->GetRMS();
0624       double yield = h->Integral();
0625   
0626       double low = -0.01, high=0.01;
0627       if(n_maps_layer == 0)
0628     {
0629       low = 3.0 * low;
0630       high = 3.0 * high;
0631     }
0632 
0633       if(pthi < 6.0)
0634     {
0635       low = 3.0*low;
0636       high = 3.0*high;
0637     }
0638 
0639       TF1 *f = new TF1("f","gaus",low, high);
0640       f->SetParameter(1, yield/100.0);
0641       f->SetParameter(2, 0.0);
0642       f->SetParameter(3,0.002);
0643       f->SetLineColor(kRed);
0644       h->Fit(f,"R");
0645 
0646       char fitr[500];
0647       sprintf(fitr,"p_{T} = %.1f-%.1f GeV/c",ptlo,pthi);
0648       TLatex *l1 = new TLatex(0.2,0.971,fitr);
0649       l1->SetNDC();
0650       l1->SetTextSize(0.07);
0651       l1->Draw();
0652 
0653       sprintf(fitr,"#sigma = %.1f #mum", f->GetParameter(2)*10000);
0654       TLatex *l1a = new TLatex(0.2,0.892,fitr);
0655       l1a->SetNDC();
0656       l1a->SetTextSize(0.07);
0657       l1a->Draw();
0658 
0659        cout << " pT range " << ptlo << " " << pthi << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
0660     }
0661 
0662   //================================
0663   // extract three DCA Z distributions inside pT bins
0664 
0665   TCanvas *cdcazdist = new TCanvas("cdcazdist","cdcazdist",5,20,1200,800);
0666   cdcazdist->Divide(3,1);
0667 
0668   for(int i = 0;i<3;i++)
0669     {
0670       double ptlo = dcaptbinlo[i];
0671       double pthi = dcaptbinhi[i];
0672 
0673       int binlo = hpt_dcaZ->GetXaxis()->FindBin(ptlo);
0674       int binhi = hpt_dcaZ->GetXaxis()->FindBin(pthi);
0675 
0676       std::cout << "DCA Z distr: ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
0677 
0678       TH1D *h = new TH1D("h","dcaZ resolution",2000, -0.1, 0.1);
0679       hpt_dcaZ->ProjectionY("h",binlo,binhi);
0680       h->GetXaxis()->SetTitle("DCA (Z) (cm)");
0681       h->GetXaxis()->SetTitleOffset(1.0);
0682 
0683       cdcazdist->cd(i+1);
0684       gPad->SetLogy(1);
0685       h->Sumw2();
0686       h->Rebin(4);
0687       h->SetMarkerSize(0.8);
0688       h->GetXaxis()->SetNdivisions(505);
0689       h->SetMinimum(0.5);
0690       h->DrawCopy("E");
0691 
0692       double mean = h->GetMean();
0693       double sigma = h->GetRMS();
0694       double yield = h->Integral();
0695   
0696       double low = -0.01, high=0.01;
0697       if(n_maps_layer == 0)
0698     {
0699       low = 3.0 * low;
0700       high = 3.0 * high;
0701     }
0702 
0703       if(pthi < 6.0)
0704     {
0705       low = 3.0*low;
0706       high = 3.0*high;
0707     }
0708 
0709       TF1 *f = new TF1("f","gaus",low, high);
0710       f->SetParameter(1, yield/100.0); 
0711      f->SetParameter(2, 0.0);
0712       f->SetParameter(3,0.002);
0713       f->SetLineColor(kRed);
0714       h->Fit(f,"R");
0715 
0716       char fitr[500];
0717       sprintf(fitr,"p_{T} = %.1f-%.1f GeV/c",ptlo,pthi);
0718       TLatex *l1 = new TLatex(0.2,0.971,fitr);
0719       l1->SetNDC();
0720       l1->SetTextSize(0.07);
0721       l1->Draw();
0722 
0723       sprintf(fitr,"#sigma = %.1f #mum", f->GetParameter(2)*10000);
0724       TLatex *l1a = new TLatex(0.2,0.892,fitr);
0725       l1a->SetNDC();
0726       l1a->SetTextSize(0.07);
0727       l1a->Draw();
0728 
0729        cout << " pT range " << ptlo << " " << pthi << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
0730     }
0731 
0732   /*
0733 
0734   //=========================
0735 
0736   TCanvas *cdca = new TCanvas("cdca","cdca",5,20,1200,800);
0737   //cdca->SetTopMargin(0.2);
0738   cdca->Divide(3,1);
0739 
0740   TH1D *hdca2d[3];
0741   fin->GetObject("hdca2d0",hdca2d[0]);
0742   fin->GetObject("hdca2d1",hdca2d[1]);
0743   fin->GetObject("hdca2d2",hdca2d[2]);
0744 
0745   cdca->cd(1);
0746   gPad->SetLogy(1);
0747   //gPad->SetRightMargin(0.1);
0748   hdca2d[0]->GetXaxis()->SetNdivisions(505);
0749 
0750   hdca2d[0]->Draw();
0751 
0752   cdca->cd(2);
0753   gPad->SetLogy(1);
0754   hdca2d[1]->GetXaxis()->SetNdivisions(505);
0755   hdca2d[1]->Draw();
0756 
0757   cdca->cd(3);
0758   gPad->SetLogy(1);
0759   hdca2d[2]->GetXaxis()->SetNdivisions(505);
0760   hdca2d[2]->Draw();
0761 
0762   TF1 *fdca = new TF1("fdca","gaus",-0.1,0.1);
0763   fdca->SetLineColor(kRed);
0764   cdca->cd(1);
0765   hdca2d[0]->Fit(fdca);
0766 
0767   TLatex *l1a = new TLatex(0.2,0.971,"p_{T} = 0.5-1.0 GeV/c");
0768   l1a->SetNDC();
0769   l1a->SetTextSize(0.07);
0770   l1a->Draw();
0771   char fitr[500];
0772   sprintf(fitr,"#sigma = %.1f #mum",fdca->GetParameter(2)*10000);
0773   TLatex *l1 = new TLatex(0.2,0.892,fitr);
0774   l1->SetNDC();
0775   l1->SetTextSize(0.07);
0776   l1->Draw();
0777 
0778   cdca->cd(2);
0779   hdca2d[1]->Fit(fdca);
0780 
0781   TLatex *l2a = new TLatex(0.2,0.971,"p_{T} = 1.0-2.0 GeV/c");
0782   l2a->SetNDC();
0783   l2a->SetTextSize(0.07);
0784   l2a->Draw();
0785 
0786   sprintf(fitr,"#sigma = %.1f #mum",fdca->GetParameter(2)*10000);
0787   TLatex *l2 = new TLatex(0.2,0.892,fitr);
0788   l2->SetNDC();
0789   l2->SetTextSize(0.07);
0790   l2->Draw();
0791 
0792 
0793   cdca->cd(3);
0794   hdca2d[2]->Fit(fdca);
0795   sprintf(fitr,"#splitline{p_{T} > 2.0 GeV/c}{#sigma = %.1f #mum}",fdca->GetParameter(2)*10000);
0796 
0797   TLatex *l3a = new TLatex(0.2,0.971,"p_{T} = 1.0-2.0 GeV/c");
0798   l3a->SetNDC();
0799   l3a->SetTextSize(0.07);
0800   l3a->Draw();
0801 
0802   sprintf(fitr,"#sigma = %.1f #mum",fdca->GetParameter(2)*10000);
0803   TLatex *l3 = new TLatex(0.2,0.892,fitr);
0804   l3->SetNDC();
0805   l3->SetTextSize(0.07);
0806   l3->Draw();
0807   */
0808 
0809   TCanvas *ceta = new TCanvas("ceta","ceta",10,10,600,600);
0810 
0811   TH1D *hgeta = 0;
0812   fin->GetObject("hgeta",hgeta);
0813   if(!hgeta)
0814     {
0815       cout << "Did not get hgeta" << endl;
0816       exit(1);
0817     }
0818   TH1D *hreta = 0;
0819   fin->GetObject("hreta",hreta);
0820 
0821   hreta->GetXaxis()->SetTitle("#eta");
0822   hreta->GetYaxis()->SetNdivisions(505);
0823   hreta->SetLineColor(kRed);
0824   hreta->SetMinimum(0.0);
0825   hreta->Draw();
0826 
0827   hgeta->Draw("same");
0828 
0829 
0830  TLegend *leta = new TLegend(0.3, 0.35, 1.0, 0.60, dcalab, "NDC");
0831   leta->SetBorderSize(0);
0832   leta->SetFillColor(0);
0833   leta->SetFillStyle(0);
0834   leta->AddEntry(hgeta,"Truth","l");
0835   leta->AddEntry(hreta,"Reconstructed","l");
0836   leta->Draw();
0837   ceta->Update();
0838 
0839   TCanvas *cfake = new TCanvas("cfake","cfake",4,4,800,600);
0840   cfake->SetLeftMargin(0.15);
0841   TH2D *hpt_nfake = 0;
0842   fin->GetObject("hpt_nfake",hpt_nfake);
0843   if(!hpt_nfake)
0844     {
0845       cout << "Did not get hpt_nfake" << endl;
0846       exit(1);
0847     }
0848 
0849   int n_noise = 4;
0850   double ynoise[10][200];
0851   double ynoise_ref[200];
0852   //double dptfake = 1.0;
0853   double dptfake = 0.2;
0854   double fake_ptmax = 20.0;
0855   int nptfake = (int) (fake_ptmax/dptfake);
0856   double ptfake[200];
0857   for(int ipt=0;ipt<nptfake;ipt++)
0858     {
0859       double ptlo = dptfake * ipt;
0860       double  pthi = dptfake*ipt + dptfake;
0861       ptfake[ipt] = (ptlo+pthi) / 2.0;
0862       int binlo = hpt_nfake->GetXaxis()->FindBin(ptlo);
0863       int binhi = hpt_nfake->GetXaxis()->FindBin(pthi);
0864       
0865       TH1D *hnoise = new TH1D("hnoise","",73,-5,68);      
0866       hpt_nfake->ProjectionY("hnoise",binlo,binhi);
0867       double ynoise_ref = hnoise->Integral();
0868       for(int j=0;j<n_noise;j++)
0869     {
0870       int bin = hnoise->FindBin(j);
0871       ynoise[j][ipt]= hnoise->GetBinContent(bin);
0872       ynoise[j][ipt] /= ynoise_ref;
0873     }
0874       hnoise->Delete();
0875     }
0876   
0877 
0878   TH1D *hfdummy = new TH1D("hfdummy","",100,0.0,fake_ptmax);
0879   hfdummy->SetMaximum(1.05);
0880   hfdummy->SetMinimum(0.001);
0881   hfdummy->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0882   hfdummy->GetYaxis()->SetTitle("Fraction of tracks");
0883   gPad->SetLogy(1);
0884   hfdummy->Draw();
0885 
0886   TLegend *lfake = new TLegend(0.35,0.55,0.87,0.89,"","NDC");
0887   lfake->SetBorderSize(1);
0888   TGraph *hnoise[10];
0889   int fake_col[5] = {kBlack,kRed,kBlue,kMagenta,kViolet};
0890   for(int j=0;j<n_noise;j++)
0891     {
0892       hnoise[j] = new TGraph(nptfake,ptfake,ynoise[j]);      
0893       hnoise[j]->SetMarkerStyle(20);
0894       hnoise[j]->SetMarkerSize(1.0);
0895       hnoise[j]->SetMarkerColor(fake_col[j]);
0896       if(j == 0)
0897     hnoise[j]->Draw("p");
0898       else
0899     hnoise[j]->Draw("p same");  
0900       char lab[500];
0901       sprintf(lab,"%i mismatched hits",j);
0902       lfake->AddEntry(hnoise[j],lab,"p");
0903     }
0904   lfake->Draw();
0905 
0906   TH2D *hcorr_nfake_nmaps = 0;
0907   fin->GetObject("hcorr_nfake_nmaps",hcorr_nfake_nmaps);
0908   if(hcorr_nfake_nmaps)
0909     {
0910       TCanvas *cmm = new TCanvas("cmm","cmm",4,4,800,600);
0911       cmm->SetLeftMargin(0.15);
0912       
0913       gPad->SetLogy(1);
0914       
0915       TLegend *lmm = new TLegend(0.35,0.62,0.9,0.90,"","NDC");
0916       lmm->SetBorderSize(1);
0917       for(int imaps =0;imaps<n_maps_layer+1;imaps++)
0918     {
0919       double nmapslo = -0.5  + 1.0 * imaps;
0920       double nmapshi = nmapslo + 1.0;
0921       
0922       int binlo = hcorr_nfake_nmaps->GetYaxis()->FindBin(nmapslo);
0923       int binhi = hcorr_nfake_nmaps->GetYaxis()->FindBin(nmapshi);
0924       
0925       TH1D *h = new TH1D("h","h",40,-1,4);
0926       hcorr_nfake_nmaps->ProjectionX("h",binlo,binhi);
0927       
0928       h->SetMarkerStyle(20);
0929       h->SetMarkerSize(1.5);
0930       h->SetMarkerColor(fake_col[imaps]);
0931       h->SetMinimum(10);      
0932       h->SetMaximum(5.0e6);   
0933       char lab[500];
0934       sprintf(lab,"%i missed maps layers",imaps);
0935       lmm->AddEntry(h,lab,"p");
0936       
0937       if(imaps == 0)
0938         {
0939           h->GetXaxis()->SetTitle("Mismatched hits");
0940           h->GetYaxis()->SetTitle("Tracks");
0941           h->DrawCopy("p");
0942           
0943         }
0944       else
0945         h->DrawCopy("same p");
0946     }
0947       
0948       lmm->Draw();
0949     }
0950     
0951   TCanvas *cvtx = new TCanvas("cvtx","cvtx",4,4,800,600);
0952   TH1D *hzevt = 0;
0953   fin->GetObject("hzevt",hzevt);
0954   if(!hzevt)
0955     {
0956       cout << "Did not get hzevt" << endl;
0957       exit(1);
0958     }
0959   hzevt->Draw();
0960 
0961   TCanvas *cvtx_res = new TCanvas("cvtx_res","cvtx_res",4,4,1200,800);
0962   cvtx_res->Divide(2,2);
0963 
0964   TH1D *hzvtx_res = 0;
0965   fin->GetObject("hzvtx_res",hzvtx_res);
0966   if(!hzvtx_res)
0967     {
0968       cout << "Did not get hzvtx_res" << endl;
0969     }
0970 
0971   if(hzvtx_res)
0972     {
0973       cvtx_res->cd(1);
0974 
0975       hzvtx_res->SetNdivisions(505);
0976       hzvtx_res->Draw();
0977       
0978       TF1 *fvtxres = new TF1("fvtxres","gaus");
0979       //hzvtx_res->Fit(fvtxres);
0980     }
0981 
0982   TH1D *hxvtx_res = 0;
0983   fin->GetObject("hxvtx_res",hxvtx_res);
0984   if(!hxvtx_res)
0985     {
0986       cout << "Did not get hxvtx_res" << endl;
0987     }
0988 
0989   if(hxvtx_res)
0990     {
0991       cvtx_res->cd(2);
0992 
0993       hxvtx_res->SetNdivisions(505);
0994       hxvtx_res->Draw();
0995       
0996       TF1 *fvtxres = new TF1("fvtxres","gaus");
0997       //hxvtx_res->Fit(fvtxres);
0998     }
0999 
1000   TH1D *hyvtx_res = 0;
1001   fin->GetObject("hyvtx_res",hyvtx_res);
1002   if(!hyvtx_res)
1003     {
1004       cout << "Did not get hyvtx_res" << endl;
1005     }
1006 
1007   if(hyvtx_res)
1008     {
1009       cvtx_res->cd(3);
1010 
1011       hyvtx_res->SetNdivisions(505);
1012       hyvtx_res->Draw();
1013       
1014       TF1 *fvtxres = new TF1("fvtxres","gaus");
1015       //hyvtx_res->Fit(fvtxres);
1016     }
1017 
1018   TH1D *hdcavtx_res = 0;
1019   fin->GetObject("hdcavtx_res",hdcavtx_res);
1020   if(!hdcavtx_res)
1021     {
1022       cout << "Did not get hdcavtx_res" << endl;
1023     }
1024 
1025   if(hdcavtx_res)
1026     {
1027       cvtx_res->cd(4);
1028 
1029       hdcavtx_res->SetNdivisions(505);
1030       hdcavtx_res->Draw();
1031       
1032       TF1 *fvtxres = new TF1("fvtxres","gaus");
1033       //hdcavtx_res->Fit(fvtxres);
1034     }
1035 
1036 
1037   TCanvas *g4ntr = new TCanvas("g4ntr","g4ntr",5,10,1200,800);
1038   g4ntr->Divide(3,1);
1039   g4ntr->cd(1);
1040   TH2D *hg4ntrack = 0;
1041   fin->GetObject("hg4ntrack",hg4ntrack);
1042   if(hg4ntrack)
1043     {
1044       hg4ntrack->GetXaxis()->SetTitle(" Tracks/event from truth");
1045       hg4ntrack->GetYaxis()->SetTitle(" Tracks/event reconstructed");
1046       hg4ntrack->Draw();
1047     }
1048   else
1049     cout << "Did not get hgn4track" << endl;
1050 
1051   g4ntr->cd(2);
1052   TH1D *hg4ntr = new TH1D("hg4ntr","hg4ntr",200, 0.0, 2000.0);
1053   hg4ntrack->ProjectionX("hg4ntr");
1054   hg4ntr->GetXaxis()->SetTitle(" Tracks/event from truth");
1055   hg4ntr->GetYaxis()->SetTitle(" Yield");
1056   hg4ntr->Draw();
1057 
1058   g4ntr->cd(3);
1059   double cent[20];
1060   double nhits[20];
1061   nhits[0] = 0;
1062   double ytot = hg4ntr->Integral(); 
1063   for(int i=1;i<20;i++)
1064     {
1065       nhits[i] = i*120.0;
1066       int binhi = hg4ntr->FindBin(nhits[i]);
1067       double y = hg4ntr->Integral(1,binhi);
1068       cent[i] = y / ytot;
1069       cout << "nhits = " << nhits[i] << " binhi = " << binhi << " y = " << y << " ytot = " << ytot << " cent = " << cent[i] << endl;
1070     }
1071   TGraph *gy = new TGraph(20,nhits,cent);
1072   gy->GetXaxis()->SetTitle(" Tracks/event from truth");
1073   gy->GetYaxis()->SetTitle(" integral fraction of events");
1074   gy->Draw();
1075 
1076   TF1 *fgy = new TF1("fgy","[0]*x+[1]*x*x+[2]*x*x*x",0,2300);
1077   gy->Fit(fgy,"R");
1078 
1079   // Output some graphs for comparison plots
1080 
1081  TFile *fout = new TFile("root_files/look_purity_out.root","recreate");
1082   gr_eff->Write();
1083   grdca2d->Write();
1084   grdcaZ->Write();
1085   grdpt->Write();
1086 
1087 }