ac903f1b |
1 | //____________________________________________________________________ |
2 | // |
3 | // AliITSMultReconstructor - find clusters in the pixels (theta and |
4 | // phi) and tracklets. |
5 | // |
6 | // These can be used to extract charged particles multiplcicity from the ITS. |
7 | // |
8 | // A tracklet consist of two ITS clusters, one in the first pixel |
9 | // layer and one in the second. The clusters are associates if the |
10 | // differencies in Phi (azimuth) and Zeta (longitudinal) are inside |
11 | // a fiducial volume. In case of multiple candidates it is selected the |
12 | // candidate with minimum distance in Phi. |
13 | // The parameter AssociationChoice allows to control if two clusters |
14 | // in layer 2 can be associated to the same cluster in layer 1 or not. |
15 | // |
16 | // ----------------------------------------------------------------- |
17 | // |
18 | // TODO: |
19 | // |
20 | // - Introduce a rough pt estimation from the difference in phi ? |
21 | // - Allow for a more refined selection criterium in case of multiple |
22 | // candidates (for instance by introducing weights for the difference |
23 | // in Phi and Zeta). |
24 | // |
25 | //____________________________________________________________________ |
26 | |
27 | #include "AliITSMultReconstructor.h" |
28 | |
29 | #include "TTree.h" |
30 | #include "TH1F.h" |
31 | #include "TH2F.h" |
32 | |
33 | |
b51872de |
34 | #include "AliITSRecPoint.h" |
ac903f1b |
35 | #include "AliITSgeom.h" |
36 | #include "AliLog.h" |
37 | |
38 | //____________________________________________________________________ |
0762f3a8 |
39 | ClassImp(AliITSMultReconstructor) |
ac903f1b |
40 | |
41 | //____________________________________________________________________ |
42 | AliITSMultReconstructor::AliITSMultReconstructor() { |
43 | |
44 | fGeometry =0; |
45 | |
46 | SetHistOn(); |
47 | SetPhiWindow(); |
48 | SetZetaWindow(); |
49 | SetOnlyOneTrackletPerC2(); |
50 | |
51 | fClustersLay1 = new Float_t*[300000]; |
52 | fClustersLay2 = new Float_t*[300000]; |
53 | fTracklets = new Float_t*[300000]; |
54 | fAssociationFlag = new Bool_t[300000]; |
55 | |
56 | for(Int_t i=0; i<300000; i++) { |
57 | fClustersLay1[i] = new Float_t[3]; |
58 | fClustersLay2[i] = new Float_t[3]; |
59 | fTracklets[i] = new Float_t[3]; |
60 | fAssociationFlag[i] = kFALSE; |
61 | } |
62 | |
63 | // definition of histograms |
64 | fhClustersDPhi = new TH1F("dphi", "dphi", 200,-0.1,0.1); |
65 | fhClustersDPhi->SetDirectory(0); |
66 | fhClustersDTheta = new TH1F("dtheta","dtheta",200,-0.1,0.1); |
67 | fhClustersDTheta->SetDirectory(0); |
68 | fhClustersDZeta = new TH1F("dzeta","dzeta",200,-0.2,0.2); |
69 | fhClustersDZeta->SetDirectory(0); |
70 | |
71 | fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",200,-0.1,0.1,200,-0.1,0.1); |
72 | fhDPhiVsDThetaAll->SetDirectory(0); |
73 | fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",200,-0.1,0.1,200,-0.1,0.1); |
74 | fhDPhiVsDThetaAcc->SetDirectory(0); |
75 | |
76 | } |
77 | |
78 | |
79 | //____________________________________________________________________ |
80 | void |
81 | AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) { |
82 | // |
83 | // - calls LoadClusterArray that finds the position of the clusters |
84 | // (in global coord) |
85 | // - convert the cluster coordinates to theta, phi (seen from the |
86 | // interaction vertex). The third coordinate is used for .... |
87 | // - makes an array of tracklets |
88 | // |
89 | // After this method has been called, the clusters of the two layers |
90 | // and the tracklets can be retrieved by calling the Get'er methods. |
91 | |
ac903f1b |
92 | // reset counters |
93 | fNClustersLay1 = 0; |
94 | fNClustersLay2 = 0; |
95 | fNTracklets = 0; |
96 | |
97 | // loading the clusters |
98 | LoadClusterArrays(clusterTree); |
99 | |
100 | // find the tracklets |
101 | AliDebug(1,"Looking for tracklets... "); |
102 | |
103 | //########################################################### |
104 | // Loop on layer 1 : finding theta, phi and z |
105 | for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) { |
106 | Float_t x = fClustersLay1[iC1][0] - vtx[0]; |
107 | Float_t y = fClustersLay1[iC1][1] - vtx[1]; |
108 | Float_t z = fClustersLay1[iC1][2] - vtx[2]; |
109 | |
110 | Float_t r = TMath::Sqrt(TMath::Power(x,2) + |
111 | TMath::Power(y,2) + |
112 | TMath::Power(z,2)); |
113 | |
114 | fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta |
115 | fClustersLay1[iC1][1] = TMath::ATan(y/x); // Store Phi |
116 | fClustersLay1[iC1][2] = z/r; // Store scaled z |
117 | } |
118 | |
119 | // Loop on layer 2 : finding theta, phi and r |
120 | for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) { |
121 | Float_t x = fClustersLay2[iC2][0] - vtx[0]; |
122 | Float_t y = fClustersLay2[iC2][1] - vtx[1]; |
123 | Float_t z = fClustersLay2[iC2][2] - vtx[2]; |
124 | |
125 | Float_t r = TMath::Sqrt(TMath::Power(x,2) + |
126 | TMath::Power(y,2) + |
127 | TMath::Power(z,2)); |
128 | |
129 | fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta |
130 | fClustersLay2[iC2][1] = TMath::ATan(y/x); // Store Phi |
131 | fClustersLay2[iC2][2] = z; // Store z |
132 | |
133 | // this only needs to be initialized for the fNClustersLay2 first associations |
134 | fAssociationFlag[iC2] = kFALSE; |
135 | } |
136 | |
137 | //########################################################### |
138 | // Loop on layer 1 |
139 | for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) { |
140 | |
141 | // reset of variables for multiple candidates |
142 | Int_t iC2WithBestPhi = 0; // reset |
143 | Float_t dPhimin = 100.; // just to put a huge number! |
144 | |
145 | // Loop on layer 2 |
146 | for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) { |
147 | |
148 | // The following excludes double associations |
149 | if (!fAssociationFlag[iC2]) { |
150 | |
151 | // find the difference in angles |
152 | Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0]; |
153 | Float_t dPhi = fClustersLay2[iC2][1] - fClustersLay1[iC1][1]; |
154 | |
155 | // find the difference in z (between linear projection from layer 1 |
156 | // and the actual point: Dzeta= z1/r1*r2 -z2) |
157 | Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]); |
158 | Float_t dZeta = fClustersLay2[iC1][2]*r2 - fClustersLay2[iC2][2]; |
159 | |
160 | if (fHistOn) { |
161 | fhClustersDPhi->Fill(dPhi); |
162 | fhClustersDTheta->Fill(dTheta); |
163 | fhClustersDZeta->Fill(dZeta); |
164 | fhDPhiVsDThetaAll->Fill(dTheta, dPhi); |
165 | } |
166 | // make "elliptical" cut in Phi and Zeta! |
167 | Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2)); |
168 | if (d>1) continue; |
169 | |
170 | //look for the minimum distance in Phi: the minimum is in iC2WithBestPhi |
171 | if (TMath::Abs(dPhi) < dPhimin) { |
172 | dPhimin = TMath::Abs(dPhi); |
173 | iC2WithBestPhi = iC2; |
174 | } |
175 | } |
176 | } // end of loop over clusters in layer 2 |
177 | |
178 | if (dPhimin<100) { // This means that a cluster in layer 2 was found that mathes with iC1 |
179 | |
180 | if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestPhi] = kTRUE; // flag the association |
181 | |
182 | // store the tracklet |
183 | |
184 | // use the average theta from the clusters in the two layers |
185 | fTracklets[fNTracklets][0] = 0.5*(fClustersLay1[iC1][0]+fClustersLay2[iC2WithBestPhi][0]); |
186 | // use the phi from the clusters in the first layer |
187 | fTracklets[fNTracklets][1] = fClustersLay1[iC1][1]; |
188 | // Store the difference between phi1 and phi2 |
189 | fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestPhi][1]; |
190 | fNTracklets++; |
191 | |
192 | AliDebug(1,Form(" Adding tracklet candidate %d (cluster %d of layer 1 and %d of layer 2)", fNTracklets, iC1)); |
193 | } |
194 | } // end of loop over clusters in layer 1 |
195 | |
196 | AliDebug(1,Form("%d tracklets found", fNTracklets)); |
197 | } |
198 | |
199 | //____________________________________________________________________ |
200 | void |
201 | AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) { |
202 | // This method |
203 | // - gets the clusters from the cluster tree |
204 | // - convert them into global coordinates |
205 | // - store them in the internal arrays |
206 | |
207 | AliDebug(1,"Loading clusters ..."); |
208 | |
209 | fNClustersLay1 = 0; |
210 | fNClustersLay2 = 0; |
211 | |
b51872de |
212 | TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint"); |
213 | TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints"); |
ac903f1b |
214 | itsClusterBranch->SetAddress(&itsClusters); |
215 | |
216 | Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries(); |
217 | |
218 | // loop over the its subdetectors |
219 | for (Int_t iIts=0; iIts < nItsSubs; iIts++) { |
220 | |
221 | if (!itsClusterTree->GetEvent(iIts)) |
222 | continue; |
223 | |
224 | Int_t nClusters = itsClusters->GetEntriesFast(); |
225 | |
226 | // stuff needed to get the global coordinates |
227 | Double_t rot[9]; fGeometry->GetRotMatrix(iIts,rot); |
228 | Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det); |
229 | Float_t tx,ty,tz; fGeometry->GetTrans(lay,lad,det,tx,ty,tz); |
230 | |
231 | // Below: |
232 | // "alpha" is the angle from the global X-axis to the |
233 | // local GEANT X'-axis ( rot[0]=cos(alpha) and rot[1]=sin(alpha) ) |
234 | // "phi" is the angle from the global X-axis to the |
235 | // local cluster X"-axis |
236 | |
237 | Double_t alpha = TMath::ATan2(rot[1],rot[0])+TMath::Pi(); |
238 | Double_t itsPhi = TMath::Pi()/2+alpha; |
239 | |
240 | if (lay==1) itsPhi+=TMath::Pi(); |
241 | Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi); |
242 | Double_t r=tx*cp+ty*sp; |
243 | |
244 | // loop over clusters |
245 | while(nClusters--) { |
b51872de |
246 | AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters); |
ac903f1b |
247 | |
248 | if (cluster->GetLayer()>1) |
249 | continue; |
250 | |
251 | Float_t x = r*cp - cluster->GetY()*sp; |
252 | Float_t y = r*sp + cluster->GetY()*cp; |
253 | Float_t z = cluster->GetZ(); |
254 | |
255 | if (cluster->GetLayer()==0) { |
256 | fClustersLay1[fNClustersLay1][0] = x; |
257 | fClustersLay1[fNClustersLay1][1] = y; |
258 | fClustersLay1[fNClustersLay1][2] = z; |
259 | fNClustersLay1++; |
260 | } |
261 | if (cluster->GetLayer()==1) { |
262 | fClustersLay2[fNClustersLay2][0] = x; |
263 | fClustersLay2[fNClustersLay2][1] = y; |
264 | fClustersLay2[fNClustersLay2][2] = z; |
265 | fNClustersLay2++; |
266 | } |
267 | |
268 | }// end of cluster loop |
269 | } // end of its "subdetector" loop |
270 | |
271 | AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2)); |
272 | } |
273 | //____________________________________________________________________ |
274 | void |
275 | AliITSMultReconstructor::SaveHists() { |
276 | |
277 | if (!fHistOn) |
278 | return; |
279 | |
ac903f1b |
280 | fhClustersDPhi->Write(); |
281 | fhClustersDTheta->Write(); |
282 | fhClustersDZeta->Write(); |
283 | fhDPhiVsDThetaAll->Write(); |
284 | fhDPhiVsDThetaAcc->Write(); |
285 | } |