-
Notifications
You must be signed in to change notification settings - Fork 1
/
testDetectorR5.C
129 lines (104 loc) · 4.87 KB
/
testDetectorR5.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "R5Detector.h"
#include "AliESDUtils.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include <TBits.h>
#include <iostream>
#endif
// You have to load the class before ... ;-)
// .L DetectorK.cxx++
//void standardPlots() {
R5Detector* its = 0;
R5Detector* CreateDetector();
AliESDEvent* esdEv = 0;
void testDetectorR5(int nev = 10, int ntracks=10) {
its = CreateDetector();
esdEv = new AliESDEvent();
esdEv->CreateStdContent();
// create dummy vertex for beam diamond
const double kBeamSig = 50e-4;
const double beamPos[3] = {0.,0.,0.};
const double beamSig[3] = {kBeamSig, kBeamSig, 6. };
AliESDVertex diamond((double*)beamPos,(double*)beamSig,"diamond");
for (int iev=0;iev<nev;iev++) {
esdEv->Reset();
esdEv->SetMagneticField(its->GetBField()*10.);
esdEv->SetDiamond(&diamond);
// generate random vertex
double vx = gRandom->Gaus(diamond.GetX(), diamond.GetXRes());
double vy = gRandom->Gaus(diamond.GetY(), diamond.GetYRes());
double vz = gRandom->Gaus(diamond.GetZ(), diamond.GetZRes());
for (int i=0;i<ntracks;i++) {
//
double eta = 2.*gRandom->Rndm()-1.; // random eta between -1:1
double pT = 2.*gRandom->Rndm()+0.2; // random pT between 0.2 and 2.2
double phi = gRandom->Rndm()*TMath::Pi()*2;
int charge = gRandom->Rndm()>0.5 ? 1:-1;
Bool_t res = its->ProcessTrack(pT, eta, 0.139, charge, phi, vx, vy, vz);
printf("inp %e %e %d %e | NHits assigned: %d\n",pT,eta,charge,phi, its->GetNHitsAssigned());
if (!res) continue; // track is not reconstructed
// transfer the track to ESDevent
AliESDtrack* esdTr = (AliESDtrack*)its->GetProbeTrackInwardAsESDTrack();
// if needed, add fake TPC flags, though this may create problems in AOD filtering
esdTr->SetStatus(AliESDtrack::kTPCin|AliESDtrack::kTPCout|AliESDtrack::kTPCrefit);
// add extra info if needed, for instance, MC label
esdTr->SetLabel(i);
its->AddESDTrackToEvent(esdEv, esdTr); // Call this for proper transfer of info to ESD event
}
// fit vertexTracks
AliESDUtils::RefitESDVertexTracks(esdEv);
int ntr = esdEv->GetNumberOfTracks();
printf("Event %d: ntr = %d\n",iev,ntr);
printf("Generated vertex: %+e %+e %+e\n",vx,vy,vz);
esdEv->GetPrimaryVertexTracks()->Print();
for (int itr=0;itr<ntr;itr++) {
AliESDtrack* esdtr = esdEv->GetTrack(itr);
printf("Tr#%2d Pt: %5.2f Eta: %+4.1f Phi: %+5.2f",itr,esdtr->Pt(), esdtr->Eta(), esdtr->Phi());
// since esd track has can account for 7 ITS layers pattern at mose, I stored
// the hits pattern info in the TPC data...
const TBits &hits = esdtr->GetTPCClusterMap();
//TBits &fakes = esdtr->GetTPCSharedMap(); // here we store fakes, but at the moment they are not set
printf(" Hits: ");
for (int ilr=0;ilr<its->GetNActiveLayers();ilr++) {
printf("%c", hits.TestBitNumber(ilr) ? '+':'-');
}
printf("\n");
}
}
}
R5Detector* CreateDetector()
{
AliESDtrack::OnlineMode(kTRUE); // to avoid friend track creation
R5Detector* det = new R5Detector("ALICE","ITS");
det->SetPropagateToOrigin(kTRUE); // if we want all tracks to be propagated to DCA to 0/0/0.
det->SetBField(1.);
// new ideal Pixel properties?
Double_t x0IB = 0.001;
Double_t x0OB = 0.005;
Double_t xRho = 0.;
Double_t resRPhiIB = 0.0001;
Double_t resZIB = 0.0001;
Double_t resRPhiOB = 0.0005;
Double_t resZOB = 0.0005;
Double_t eff = 0.98;
//
// select Z span in such a way to have +-1 unit eta coverage for vertex at 2sigmaZ (~12cm) from nominal IP
// i.e. Zmax >= 12 + R/tan( 2*atan(exp(-1.)) ) = 12 + R/0.851
det->AddLayer((char*)"vertex", 0.0, 0.1, 0, 0); // dummy vertex for matrix calculation
det->AddLayer((char*)"bpipe", 1.6, 200., 0.0022);
det->AddLayer((char*)"ddd1", 1.8, 21.0, x0IB, xRho, resRPhiIB, resZIB,eff);
det->AddLayer((char*)"ddd2", 2.8, 21.0, x0IB, xRho, resRPhiIB, resZIB,eff);
det->AddLayer((char*)"ddd3", 3.8, 21.0, x0IB, xRho, resRPhiIB, resZIB,eff);
det->AddLayer((char*)"ddd3a", 8.0, 21.0, x0IB, xRho, resRPhiOB, resZOB,eff);
det->AddLayer((char*)"ddd4", 20.0, 42.0, x0OB, xRho, resRPhiOB, resZOB,eff);
det->AddLayer((char*)"ddd5", 25.0, 42.0, x0OB, xRho, resRPhiOB, resZOB,eff);
//det->AddLayer((char*)"ddd6", 35.0, 80.0, x0OB, xRho, resRPhiOB, resZOB,eff);
det->AddLayer((char*)"ddd7", 40.0, 80.0, x0OB, xRho, resRPhiOB, resZOB,eff);
det->AddLayer((char*)"ddd8", 55.0, 80.0, x0OB, xRho, resRPhiOB, resZOB,eff);
// det->AddLayer((char*)"dddZ", 90., 130., x0OB, xRho, resRPhiOB, resZOB,eff);
det->AddLayer((char*)"dddY", 80.0, 130.0, x0OB, xRho, resRPhiOB, resZOB,eff);
det->AddLayer((char*)"dddX", 100.0, 130.0, x0OB, xRho, resRPhiOB, resZOB,eff);
return det;
}