ITS/UPGRADE/testDetKMC.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TCanvas.h>
   3 #include <TString.h>
   4 #include <TROOT.h>
   5 #include <TSystem.h>
   6 #include <TFile.h>
   7 #include <TObjArray.h>
   8 #include <TH1.h>
   9 #include <TH2.h>
  10 #include <TF1.h>
  11 #include "KMCDetector.h"
  12 #endif
  13
  14 TObjArray* CreateTrackConditions(const char* nm);
  15 TObjArray  summary;
  16
  17
  18 TObject* gDet=0;
  19
  20 // ptBins to generate
  21 double pts[] = {0.15,0.175,0.2,0.225,0.25,0.3,0.4,0.5,0.8,1.1,1.6,2,4,8,13,20};
  22 //double pts[] = {0.2, 1., 10.};
  23
  24 void testDetKMC(int nev=10000,  // n events to generate per bin
  25                 int setVer=0,   // optional version to pass for detector building
  26                 double mass=0.14,  // particle mass to generate
  27                 double eta=0.45,   // eta to test
  28                 Bool_t aliceNew=kTRUE,Bool_t tpc=kFALSE, Double_t eff=0.95,  // detector settings
  29                 Double_t vtxConstr=-1., // use vertex constraint in the fit
  30                 const char* sumOut="trsumDef.root" // output file name
  31                 )
  32 {
  33   if (!gROOT->GetClass("KMCDetector")) gROOT->LoadMacro("KMCDetector.cxx+");
  34   //
  35   // Detector, tracking initialization >>>
  36   //
  37   if (vtxConstr>0) KMCDetector::SetVtxConstraint(vtxConstr,vtxConstr);
  38   KMCLayer::SetDefEff(eff);
  39   //
  40   KMCDetector *itsdet = new KMCDetector((char*)"ALICE",(char*)"ITS");
  41   gDet = itsdet;
  42   KMCDetector &its = *itsdet;
  43   its.SetAvgRapidity(eta);
  44   if (aliceNew)     its.MakeAliceAllNew(tpc,1, setVer); // its sa
  45   else              its.MakeAliceCurrent(tpc,0); // its sa
  46   //its.MakeAliceCurrent(0,1); // with tpc
  47   its.SetMaxSeedToPropagate(300);
  48   its.SetUseBackground();
  49   //
  50   // min hits per track to validate it
  51   its.RequireMinITSHits( its.GetNActiveITSLayers() - 2 );//(int)TMath::Max(4.0,(3./4.)*its.GetNActiveITSLayers()));
  52   printf("Min number of hits requested: %d\n",its.GetMinITSHits());
  53   //
  54   // max cluster-track chi2
  55   its.RequireMaxChi2Cl(16.);
  56   // max chi2/NDF of the Kalman fit
  57   its.RequireMaxNormChi2NDF(5.);
  58   // penalty to chi2 from missing hit
  59   KMCProbe::SetMissingHitPenalty(6);
  60   //
  61   // Detector, tracking initialization <<<
  62
  63   its.InitMCWatchHistos();
  64   summary.Clear();
  65   //
  66   int npt = sizeof(pts)/sizeof(double);
  67   double ptMin = pts[0];
  68   double ptMax = pts[npt-1];
  69   //
  70   for (int ip=npt;ip--;) {
  71     //
  72     int nevR = nev;
  73     double pt = pts[ip];//1./(ptminI+dpt*ip);  // ptmin+ip*(ptmax-ptmin)/npt;
  74
  75     double chi2Cl = 16. - (16.-9.)*(1./pt-1./ptMax)/(1./ptMin-1./ptMax);
  76     //    if (pt<0.3) nevR = int(0.25*nevR);
  77     //    else if (pt<0.5) nevR = int(0.7*nevR);
  78     printf("Doing pt(%d)=%.3f for mass %.3f, %d ev | chi2Cl=%.2f\n",ip,pt,mass,nevR,chi2Cl);
  79     TString smn = Form("pt%d",ip);
  80     TObjArray* trSum = CreateTrackConditions(smn.Data());
  81     if (trSum) summary.AddLast(trSum);
  82     //
  83     its.SolveSingleTrack(mass, pt, eta, trSum, nevR);
  84   }
  85   //
  86   TString smout = sumOut;
  87   if (!smout.IsNull()) {
  88     gSystem->ExpandPathName(smout);
  89     TFile* fl = TFile::Open(smout.Data(),"recreate");
  90     if (!fl) {printf("Failed to open %s file\n",smout.Data()); return;}
  91     fl->WriteObject(&summary,"trSum","kSingleKey");
  92     fl->Close();
  93     delete fl;
  94   }
  95 }
  96
  97 //_____________________________________________________________
  98 TObjArray* CreateTrackConditions(const char* nm)
  99 {
 100   // create set of track criteria to extract the summary
 101   KMCTrackSummary *sm = 0;
 102   TObjArray* arr = new TObjArray();
 103   //
 104   int nlr = KMCProbe::GetNITSLayers();
 105   //-------------------- put here user code -------------------
 106   // summary for ideal correct tracks
 107   sm = new KMCTrackSummary("corrAll","corrAll",nlr);
 108   sm->SetMinMaxClITS(nlr);
 109   sm->SetMinMaxClITSFake(0,0);
 110   sm->SetNamePrefix(nm);
 111   arr->AddLast(sm);
 112   //
 113   // summary for good correct tracks
 114   sm = new KMCTrackSummary("corr4h","corr4h",nlr);
 115   sm->SetMinMaxClITS(4);
 116   sm->SetMinMaxClITSFake(0,0);
 117   sm->SetNamePrefix(nm);
 118   arr->AddLast(sm);
 119   //
 120   // summary for good correct tracks with enough points at inner single layers
 121   sm = new KMCTrackSummary("corr4hs2in","corr4hs2in",nlr);
 122   sm->SetMinMaxClITS(4);
 123   sm->SetMinMaxClITSFake(0,0);
 124   sm->SetNamePrefix(nm);
 125   sm->AddPatternITSCorr( sm->Bits(1,0));    // the innermost 2 layers must have correct cluster
 126   sm->AddPatternITSCorr( sm->Bits(0,1));
 127   arr->AddLast(sm);
 128   //
 129   // summary for good correct tracks with enough points at inner double layers
 130   sm = new KMCTrackSummary("corr4hd2in","corr5hd2in",nlr);
 131   sm->SetMinMaxClITS(4);
 132   sm->SetMinMaxClITSFake(0,0);
 133   sm->SetNamePrefix(nm);
 134   sm->AddPatternITSCorr( sm->Bits(1));    // the innermost 2 double layers must have correct cluster
 135   sm->AddPatternITSCorr( sm->Bits(0,1,1));
 136   arr->AddLast(sm);
 137   //
 138   //
 139   // summary for golden tracks with AllNew RS setup: good pattern for offset and pt measurement
 140   sm = new KMCTrackSummary("corr4hdGoodPatt","corr4hdGoodPatt",nlr);
 141   sm->SetMinMaxClITS(4);
 142   sm->SetMinMaxClITSFake(0,0);
 143   sm->SetNamePrefix(nm);
 144   sm->AddPatternITSCorr( sm->Bits(1,1));    // the innermost 2 double layers must have correct cluster
 145   sm->AddPatternITSCorr( sm->Bits(0,0,1));
 146   //
 147   sm->AddPatternITSCorr( sm->Bits(0,0,0,1,1)); // measurement in the middle (would be prefereable 2 points!)
 148   //
 149   sm->AddPatternITSCorr( sm->Bits(0,0,0, 0,0, 1,1)); // measurement at large radius
 150   //
 151   arr->AddLast(sm);
 152   //
 153   // summary for golden tracks with curr setup: good pattern for offset and pt measurement
 154   sm = new KMCTrackSummary("corr4hdGoodPatt","corr4hdGoodPatt",nlr);
 155   sm->SetMinMaxClITS(4);
 156   sm->SetMinMaxClITSFake(0,0);
 157   sm->SetNamePrefix(nm);
 158   sm->AddPatternITSCorr( sm->Bits(1,0));    // the innermost 2 double layers must have correct cluster
 159   sm->AddPatternITSCorr( sm->Bits(0,1));
 160   //
 161   sm->AddPatternITSCorr( sm->Bits(0,0,1,1)); // measurement in the middle (would be prefereable 2 points!)
 162   //
 163   sm->AddPatternITSCorr( sm->Bits(0,0, 0,0, 1,1)); // measurement at large radius
 164   //
 165   arr->AddLast(sm);
 166   //
 167   // summary for tracks with at least 1 fake cluster
 168   sm = new KMCTrackSummary("fakeAny","fakeAny",nlr);
 169   sm->SetMinMaxClITS(4);
 170   sm->SetMinMaxClITSFake(1);
 171   sm->SetNamePrefix(nm);
 172   arr->AddLast(sm);
 173   //
 174   // // summary for tracks with at least 1 fake cluster (but still acceptable?)
 175   // sm = new KMCTrackSummary("corrF","corrF",nlr);
 176   // sm->SetMinMaxClITS(4);
 177   // sm->SetMinMaxClITSCorr(3);              // at least 3 correct clusters
 178   // sm->SetMinMaxClITSFake(1);              // at least 1 fake clusters
 179   // //  /*
 180   // sm->AddPatternITSCorr( sm->Bits(1));    // the innermost 2 layers must have correct cluster
 181   // sm->AddPatternITSCorr( sm->Bits(0,1,1));
 182   // //  */
 183   // //  sm->AddPatternITSCorr( sm->Bits(1,1));    //
 184   // // sm->AddPatternITSCorr( sm->Bits(0,0,1,1,1,1));
 185
 186   // //  sm->AddPatternITSFakeExcl( sm->Bits(1,1,1)); // no fakes in innermost 3 layers
 187   // sm->SetNamePrefix(nm);
 188   // arr->AddLast(sm);
 189   //
 190   //-----------------------------------------------------------
 191   return arr->GetEntries()>0 ? arr : 0;
 192 }