]>
Commit | Line | Data |
---|---|---|
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 | } |