StRoot: StGammaMaker/macros/RunGammaFitterDemo.C Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 void RunGammaFitterDemo(const char* filelist = "/star/u/sdhamija/StGammaMaker/StRoot/StGammaMaker/macros/photon_9_11.list")
00008 {
00009   // Load shared libraries
00010   gROOT->Macro("$STAR/StRoot/StGammaMaker/macros/loadGammaLibs.C");
00011 
00012   // Create TChain from filelist
00013   TChain* chain = new TChain("gammas");
00014   ifstream in(filelist);
00015   string line;
00016   while (getline(in, line)) chain->AddFile(line.c_str());
00017   in.close();
00018 
00019   // Create gamma event buffer
00020   StGammaEvent* event = new StGammaEvent;
00021   chain->SetBranchAddress("GammaEvent", &event);
00022 
00023   // Create gamma fitter histograms
00024   TH2* hResidualYield = new TH2F("hResidualYield", "#gamma events from pp#rightarrow#gamma+jet+X;Fit Residual Sum (R_{0,u}+R_{0,v}) [MeV];Fitted Peak Integral [MeV]", 40, 0, 200, 40, 0, 600);
00025   TH1* hVertexZ = new TH1F("hVertexZ", ";z_{vertex} [cm]", 100, -200, 200);
00026   TH2* hVertexXY = new TH2F("hVertexXY", ";x_{vertex} [cm];y_{vertex} [cm]", 40, -0.2, 0.2, 40, -0.45, -0.2);
00027   TH1* hClusterEnergy = new TH1F("hClusterEnergy", ";Cluster E [GeV]", 100, 0, 60);
00028   TH1* hClusterPt = new TH1F("hClusterPt", ";Cluster p_{T} [GeV]", 100, 0, 20);
00029   TH2* hClusterXY = new TH2F("hClusterXY", ";Cluster x [cm];Cluster y [cm]", 40, -250, 250, 40, -250, 250);
00030   TH2* hClusterEtaPhi = new TH2F("hClusterEtaPhi", ";Cluster #eta;Cluster #phi [radians]", 40, 1, 2, 40, -TMath::Pi(), TMath::Pi());
00031   TH2* hSmdPointXY = new TH2F("hSmdPointXY", ";SMD x [cm]; SMD y [cm]", 40, -250, 250, 40, -250, 250);
00032   TH1* hYieldU = new TH1F("hYieldU", ";Fitted Peak Integral [MeV]", 100, 0, 400);
00033   TH1* hYieldV = new TH1F("hYieldV", ";Fitted Peak Integral [MeV]", 100, 0, 400);
00034   TH1* hResidualU = new TH1F("hResidualU", ";Fit Residual [MeV]", 100, -40, 120);
00035   TH1* hResidualV = new TH1F("hResidualV", ";Fit Residual [MeV]", 100, -40, 120);
00036   TH1* hChi2PerNdfU = new TH1F("hChi2PerNdfU", ";#chi^{2}/ndf", 100, 0, 10);
00037   TH1* hChi2PerNdfV = new TH1F("hChi2PerNdfV", ";#chi^{2}/ndf", 100, 0, 10);
00038   TH1* hMeanU = new TH1F("hMeanU", ";Mean #mu", 100, 0, 288);
00039   TH1* hMeanV = new TH1F("hMeanV", ";Mean #mu", 100, 0, 288);
00040   TH1* hSigmaU = new TH1F("hSigmaU", ";RMS #sigma", 100, 0, 10);
00041   TH1* hSigmaV = new TH1F("hSigmaV", ";RMS #sigma", 100, 0, 10);
00042   TH1* hNhitsU = new TH1F("hNhitsU", ";Number of SMD hits", 50, 0, 50);
00043   TH1* hNhitsV = new TH1F("hNhitsV", ";Number of SMD hits", 50, 0, 50);
00044 
00045   // Function for gamma/pi0 separation
00046   TF1* fResidualCut = new TF1("fResidualCut", "pol2", 0, 200);
00047   fResidualCut->SetParameters(100, 0, 0.05);
00048 
00049   // Get number of entries in gamma tree
00050   int nevents = chain->GetEntries();
00051   cout << "nevents = " << nevents << endl;
00052 
00053   // Event loop
00054   int ncandidates = 0;
00055   for (int iEvent = 0; iEvent < nevents; ++iEvent) {
00056     chain->GetEvent(iEvent);
00057     if (iEvent % 1000 == 0) cout << "iEvent = " << iEvent << endl;
00058     // Skip events without vertex
00059     if (event->vertex() == TVector3(0,0,0)) continue;
00060     // Candidate loop
00061     for (int i = 0; i < event->numberOfCandidates(); ++i) {
00062       StGammaCandidate* candidate = event->candidate(i);
00063       assert(candidate);
00064       // EEMC candidates only
00065       if (candidate->detectorId() == StGammaCandidate::kEEmc) {
00066         // Fit SMD u & v plane
00067         StGammaFitterResult fits[2];
00068         StGammaFitterResult& u = fits[0];
00069         StGammaFitterResult& v = fits[1];
00070         StGammaFitter::instance()->fit(candidate, fits, 0);
00071         StGammaFitter::instance()->fit(candidate, fits, 1);
00072         int sector, subsector, etabin;
00073         EEmcGeomSimple::Instance().getTower(candidate->position(), sector, subsector, etabin);
00074         TVector3& position = EEmcSmdGeom::instance()->getIntersection(sector, u.mean, v.mean);
00075         if (position.z() != -999) {
00076           float residual = u.residual+v.residual;
00077           float yield = u.yield+v.yield;
00078           if (residual > 0 && yield > 0) {
00079             ++ncandidates;
00080             // Fill histograms
00081             hResidualYield->Fill(residual, yield);
00082             hVertexZ->Fill(event->vertex().z());
00083             hVertexXY->Fill(event->vertex().x(), event->vertex().y());
00084             hClusterEnergy->Fill(candidate->energy());
00085             hClusterPt->Fill(candidate->momentum().Pt());
00086             hClusterXY->Fill(candidate->position().x(), candidate->position().y());
00087             hClusterEtaPhi->Fill(candidate->position().Eta(), candidate->position().Phi());
00088             hSmdPointXY->Fill(position.x(), position.y());
00089             hYieldU->Fill(u.yield);
00090             hYieldV->Fill(v.yield);
00091             hResidualU->Fill(u.residual);
00092             hResidualV->Fill(v.residual);
00093             if (u.ndf > 0) hChi2PerNdfU->Fill(u.chiSquare / u.ndf);
00094             if (v.ndf > 0) hChi2PerNdfV->Fill(v.chiSquare / v.ndf);
00095             hMeanU->Fill(u.mean);
00096             hMeanV->Fill(v.mean);
00097             hSigmaU->Fill(u.rms);
00098             hSigmaV->Fill(v.rms);
00099             hNhitsU->Fill(candidate->numberOfSmdu());
00100             hNhitsV->Fill(candidate->numberOfSmdv());
00101           }
00102         }
00103       }
00104     } // End candidate loop
00105   } // End event loop
00106 
00107   cout << "ncandidates = " << ncandidates << endl;
00108 
00109   // Create canvas
00110   TCanvas* c1 = new TCanvas("c1", "c1", 1200, 800);
00111   c1->Divide(3, 2);
00112   gStyle->SetPalette(1);
00113   gStyle->SetOptStat(11);
00114 
00115   // Draw histograms
00116   c1->cd(1);
00117   hResidualYield->Draw("colz");
00118   fResidualCut->SetLineColor(kRed);
00119   fResidualCut->Draw("same");
00120   gPad->Print("hResidualYield.png");
00121 
00122   c1->cd(2);
00123   hVertexZ->Draw();
00124   hVertexZ->Fit("gaus");
00125   hVertexZ->GetFunction("gaus")->SetLineColor(kRed);
00126   gPad->Print("hVertexZ.png");
00127 
00128   c1->cd(3);
00129   hVertexXY->Draw("colz");
00130   gPad->Print("hVertexXY.png");
00131 
00132   c1->cd(4);
00133   hClusterEnergy->Draw();
00134   gPad->Print("hClusterEnergy.png");
00135 
00136   c1->cd(5);
00137   hClusterPt->Draw();
00138   gPad->Print("hClusterPt.png");
00139 
00140   c1->cd(6);
00141   hClusterXY->Draw("colz");
00142   gPad->Print("hClusterXY.png");
00143 
00144   c1->cd(1);
00145   hClusterEtaPhi->Draw("colz");
00146   gPad->Print("hClusterEtaPhi.png");
00147 
00148   c1->cd(2);
00149   hSmdPointXY->Draw("colz");
00150   gPad->Print("hSmdPointXY.png");
00151 
00152   c1->cd(3);
00153   hYieldU->SetLineColor(kRed);
00154   hYieldV->SetLineColor(kBlue);
00155   hYieldU->SetLineWidth(2);
00156   hYieldV->SetLineWidth(2);
00157   hYieldU->Draw();
00158   hYieldV->Draw("same");
00159   TLegend* leg1 = new TLegend(0.60, 0.60, 0.85, 0.80);
00160   leg1->AddEntry(hYieldU, "SMD-u", "L");
00161   leg1->AddEntry(hYieldV, "SMD-v", "L");
00162   leg1->Draw();
00163   gPad->Print("hYield.png");
00164 
00165   c1->cd(4);
00166   hResidualU->SetLineColor(kRed);
00167   hResidualV->SetLineColor(kBlue);
00168   hResidualU->SetLineWidth(2);
00169   hResidualV->SetLineWidth(2);
00170   hResidualU->Draw();
00171   hResidualV->Draw("same");
00172   TLegend* leg2 = new TLegend(0.60, 0.60, 0.85, 0.80);
00173   leg2->AddEntry(hResidualU, "SMD-u", "L");
00174   leg2->AddEntry(hResidualV, "SMD-v", "L");
00175   leg2->Draw();
00176   gPad->Print("hResidual.png");
00177 
00178   c1->cd(5);
00179   hChi2PerNdfU->SetLineColor(kRed);
00180   hChi2PerNdfV->SetLineColor(kBlue);
00181   hChi2PerNdfU->SetLineWidth(2);
00182   hChi2PerNdfV->SetLineWidth(2);
00183   hChi2PerNdfU->Draw();
00184   hChi2PerNdfV->Draw("same");
00185   TLegend* leg3 = new TLegend(0.60, 0.60, 0.85, 0.80);
00186   leg3->AddEntry(hChi2PerNdfU, "SMD-u", "L");
00187   leg3->AddEntry(hChi2PerNdfV, "SMD-v", "L");
00188   leg3->Draw();
00189   gPad->Print("hChi2PerNdf.png");
00190 
00191   c1->cd(6);
00192   hMeanU->SetLineColor(kRed);
00193   hMeanV->SetLineColor(kBlue);
00194   hMeanU->SetLineWidth(2);
00195   hMeanV->SetLineWidth(2);
00196   hMeanU->Draw();
00197   hMeanV->Draw("same");
00198   TLegend* leg4 = new TLegend(0.60, 0.60, 0.85, 0.80);
00199   leg4->AddEntry(hMeanU, "SMD-u", "L");
00200   leg4->AddEntry(hMeanV, "SMD-v", "L");
00201   leg4->Draw();
00202   gPad->Print("hMean.png");
00203 
00204   c1->cd(1);
00205   hSigmaU->SetLineColor(kRed);
00206   hSigmaV->SetLineColor(kBlue);
00207   hSigmaU->SetLineWidth(2);
00208   hSigmaV->SetLineWidth(2);
00209   hSigmaU->Draw();
00210   hSigmaV->Draw("same");
00211   TLegend* leg5 = new TLegend(0.60, 0.60, 0.85, 0.80);
00212   leg5->AddEntry(hSigmaU, "SMD-u", "L");
00213   leg5->AddEntry(hSigmaV, "SMD-v", "L");
00214   leg5->Draw();
00215   gPad->Print("hSigma.png");
00216 
00217   c1->cd(2);
00218   hNhitsU->SetLineColor(kRed);
00219   hNhitsV->SetLineColor(kBlue);
00220   hNhitsU->SetLineWidth(2);
00221   hNhitsV->SetLineWidth(2);
00222   hNhitsU->Draw();
00223   hNhitsV->Draw("same");
00224   TLegend* leg6 = new TLegend(0.60, 0.60, 0.85, 0.80);
00225   leg6->AddEntry(hNhitsU, "SMD-u", "L");
00226   leg6->AddEntry(hNhitsV, "SMD-v", "L");
00227   leg6->Draw();
00228   gPad->Print("hNhits.png");
00229 }