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