ITS/UPGRADE/PrepSummaryKMC.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TROOT.h>
   3 #include <TFile.h>
   4 #include <TDirectory.h>
   5 #include <TF1.h>
   6 #include <TH1.h>
   7 #include <TObjArray.h>
   8 #include <TGraphErrors.h>
   9 #include "KMCDetector.h"
  10 #endif
  11
  12 // Process the output file of testDetKMC.C macro to prepare summary
  13 // for the "summary class" icl (definition of acceptable track).
  14 // See ProcessSummary comment for details.
  15 // The names of the graphs are prefixed by `pref`.
  16 // If outF is provided, the graphs will be stored there
  17
  18
  19 TObjArray* PrepSummaryKMC(const char* sumf, int cls, const char* pref=0, const char* outF=0);
  20 TObjArray* ProcessSummary(TObjArray* sums, int icl, const char* pref=0);
  21
  22 enum {kSigAD,kSigAZ,kSigAPt,kSigD,kSigZ,kSigPt,kEff,kUpd};
  23 TF1* gs = 0;
  24
  25 TObjArray* PrepSummaryKMC(const char* sumf, int cls, const char* pref, const char* outF)
  26 {
  27   if (!gROOT->GetClass("KMCDetector")) gROOT->LoadMacro("KMCDetector.cxx+");
  28   TFile* fl = TFile::Open(sumf);
  29   if (!fl) {printf("No file %s\n",sumf); return 0;}
  30   TObjArray* arrs = (TObjArray*)gDirectory->Get("trSum");
  31   if (!arrs) {printf("No summary in file %s\n",sumf); return 0;}
  32   //
  33   TObjArray* sums =  ProcessSummary(arrs,cls,pref);
  34   if (!sums) return 0;
  35   //
  36   if (outF) {
  37     TFile* flOut = TFile::Open(outF,"update");
  38     if (!flOut) {printf("Failed to open output file %s\n",outF);}
  39     else {
  40       flOut->WriteObject(sums,sums->GetName(),"kSingleKey");
  41       flOut->Close();
  42       delete flOut;
  43       printf("Stored array %s in %s\n",sums->GetName(), outF);
  44     }
  45   }
  46   //
  47   return sums;
  48 }
  49
  50
  51 TObjArray* ProcessSummary(TObjArray* arrs, int icl, const char* pref)
  52 {
  53   // Process TObjArray (e.g. for set of pt bins) of TObjArray of KMCTrackSummary objects:
  54   // pick the KMCTrackSummary for "summary class" icl (definition of acceptable track) and create
  55   // graphs vs bin.
  56   // These graphs are returned in new TObjArray
  57   //
  58   TString prefs = pref;
  59   if (!gs) gs = new TF1("gs","gaus",-1,1);
  60   //
  61   int nb = arrs->GetEntriesFast();
  62   TObjArray* sums = (TObjArray*) arrs->At(0);
  63   int nclass = sums->GetEntriesFast();
  64   if (icl>=nclass) {printf("summary set has %d classes only, %d requested\n",nclass,icl);return 0;}
  65   //
  66   KMCTrackSummary* sm = (KMCTrackSummary*)sums->At(icl);
  67   //
  68   TH1* h;
  69   //
  70   h = sm->GetHMCSigDCARPhi();  // MC resolution in transverse DCA
  71   TGraphErrors * grSigD = 0;
  72   if (h) {
  73     grSigD = new TGraphErrors(nb);
  74     grSigD->SetName(Form("%s%s",prefs.Data(),h->GetName()));
  75     grSigD->SetTitle(Form("%s%s",prefs.Data(),h->GetTitle()));
  76   }
  77   //
  78   TGraphErrors * grSigZ = 0;
  79   h = sm->GetHMCSigDCAZ();    // MC resolution in Z DCA
  80   if (h) {
  81     grSigZ = new TGraphErrors(nb);
  82     grSigZ->SetName(Form("%s%s",prefs.Data(),h->GetName()));
  83     grSigZ->SetTitle(Form("%s%s",prefs.Data(),h->GetTitle()));
  84   }
  85   //
  86   TGraphErrors * grSigAD = 0; // anaitical estimate for resolution in transverse DCA
  87   {
  88     grSigAD = new TGraphErrors(nb);
  89     grSigAD->SetName(Form("%s%s",prefs.Data(),"sigmaDan"));
  90     grSigAD->SetTitle(Form("%s%s",prefs.Data(),"#sigmaD an"));
  91   }
  92   //
  93   TGraphErrors * grSigAZ = 0; // anaitical estimate for resolution in Z DCA
  94   {
  95     grSigAZ = new TGraphErrors(nb);
  96     grSigAZ->SetName(Form("%s%s",prefs.Data(),"sigmaZan"));
  97     grSigAZ->SetTitle(Form("%s%s",prefs.Data(),"#sigmaZ an"));
  98   }
  99   //
 100   //
 101   TGraphErrors * grSigPt = 0; // MC res. in pt
 102   {
 103     grSigPt = new TGraphErrors(nb);
 104     grSigPt->SetName(Form("%s%s",prefs.Data(),"sigmaPt"));
 105     grSigPt->SetTitle(Form("%s%s",prefs.Data(),"#sigmaPt"));
 106   }
 107   //
 108   TGraphErrors * grSigAPt = 0; // analitycal res. in pt
 109   {
 110     grSigAPt = new TGraphErrors(nb);
 111     grSigAPt->SetName(Form("%s%s",prefs.Data(),"sigmaPtan"));
 112     grSigAPt->SetTitle(Form("%s%s",prefs.Data(),"#sigmaPt an"));
 113   }
 114
 115   //
 116   TGraphErrors * grEff = 0; // MC efficiency
 117   {
 118     grEff = new TGraphErrors(nb);
 119     grEff->SetName(Form("%s_rate",prefs.Data()));
 120     grEff->SetTitle(Form("%s Rate",prefs.Data()));
 121   }
 122   //
 123   TGraphErrors * grUpd = 0; // number of Kalman track updates
 124   {
 125     grUpd = new TGraphErrors(nb);
 126     grUpd->SetName(Form("%s_updCalls",prefs.Data()));
 127     grUpd->SetTitle(Form("%s Updates",prefs.Data()));
 128   }
 129   //
 130   for (int ib=0;ib<nb;ib++) {
 131     sums = (TObjArray*) arrs->At(ib);
 132     sm = (KMCTrackSummary*)sums->At(icl);
 133     KMCProbe& prbRef = sm->GetRefProbe();
 134     KMCProbe& prbAn  = sm->GetAnProbe();
 135
 136     double pt = prbRef.Pt();
 137     //
 138     if (grSigAD) {
 139       grSigAD->SetPoint(ib, pt,prbAn.GetSigmaY2()>0 ? TMath::Sqrt(prbAn.GetSigmaY2()) : 0.);
 140     }
 141     //
 142     if (grSigAZ) {
 143       grSigAZ->SetPoint(ib, pt,prbAn.GetSigmaZ2()>0 ? TMath::Sqrt(prbAn.GetSigmaZ2()) : 0.);
 144     }
 145     //
 146     if (grSigAPt) {
 147       double pts = TMath::Sqrt(prbAn.GetSigma1Pt2());
 148       grSigAPt->SetPoint(ib, pt,pts>0 ? pts*pt : 0.);
 149     }
 150     //
 151     if (grSigPt) {
 152       h = sm->GetHMCSigPt();
 153       h->Fit(gs,"0q");
 154       grSigPt->SetPoint(ib, pt, gs->GetParameter(2));
 155       grSigPt->SetPointError(ib, 0, gs->GetParError(2));
 156     }
 157     //
 158      if (grSigD) {
 159       h = sm->GetHMCSigDCARPhi();
 160       h->Fit(gs,"0q");
 161       grSigD->SetPoint(ib, pt,gs->GetParameter(2));
 162       grSigD->SetPointError(ib, 0,gs->GetParError(2));
 163     }
 164     //
 165     if (grSigZ) {
 166       h = sm->GetHMCSigDCAZ();
 167       h->Fit(gs,"0q");
 168       grSigZ->SetPoint(ib, pt,gs->GetParameter(2));
 169       grSigZ->SetPointError(ib, 0,gs->GetParError(2));
 170     }
 171     //
 172     if (grEff) {
 173       grEff->SetPoint(ib, pt,sm->GetEff());
 174       grEff->SetPointError(ib, 0,sm->GetEffErr());
 175     }
 176     //
 177     if (grUpd) {
 178       grUpd->SetPoint(ib, pt,sm->GetUpdCalls());
 179       grUpd->SetPointError(ib, 0, 0);
 180     }
 181   }
 182   //
 183   TObjArray* dest = new TObjArray();
 184   dest->AddAtAndExpand(grSigAD,kSigAD);
 185   dest->AddAtAndExpand(grSigAZ,kSigAZ);
 186   dest->AddAtAndExpand(grSigAPt,kSigAPt);
 187   dest->AddAtAndExpand(grSigD,kSigD);
 188   dest->AddAtAndExpand(grSigZ,kSigZ);
 189   dest->AddAtAndExpand(grSigPt,kSigPt);
 190   dest->AddAtAndExpand(grEff,kEff);
 191   dest->AddAtAndExpand(grUpd,kUpd);
 192   //
 193   if (!prefs.IsNull()) dest->SetName(pref);
 194   return dest;
 195 }