don't sort clusters after local reco, do this in AliITSUTrackerGlo
[u/mrichter/AliRoot.git] / 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 }