]>
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 | ||
34 | #include "AliITSclusterV2.h" | |
35 | #include "AliITSgeom.h" | |
36 | #include "AliLog.h" | |
37 | ||
38 | //____________________________________________________________________ | |
39 | ClassImp(AliITSMultReconstructor); | |
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 | ||
92 | cout << " HEEEEEEEEEEEEEEEEEEEEE " << flush << endl; | |
93 | ||
94 | // reset counters | |
95 | fNClustersLay1 = 0; | |
96 | fNClustersLay2 = 0; | |
97 | fNTracklets = 0; | |
98 | ||
99 | // loading the clusters | |
100 | LoadClusterArrays(clusterTree); | |
101 | ||
102 | // find the tracklets | |
103 | AliDebug(1,"Looking for tracklets... "); | |
104 | ||
105 | //########################################################### | |
106 | // Loop on layer 1 : finding theta, phi and z | |
107 | for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) { | |
108 | Float_t x = fClustersLay1[iC1][0] - vtx[0]; | |
109 | Float_t y = fClustersLay1[iC1][1] - vtx[1]; | |
110 | Float_t z = fClustersLay1[iC1][2] - vtx[2]; | |
111 | ||
112 | Float_t r = TMath::Sqrt(TMath::Power(x,2) + | |
113 | TMath::Power(y,2) + | |
114 | TMath::Power(z,2)); | |
115 | ||
116 | fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta | |
117 | fClustersLay1[iC1][1] = TMath::ATan(y/x); // Store Phi | |
118 | fClustersLay1[iC1][2] = z/r; // Store scaled z | |
119 | } | |
120 | ||
121 | // Loop on layer 2 : finding theta, phi and r | |
122 | for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) { | |
123 | Float_t x = fClustersLay2[iC2][0] - vtx[0]; | |
124 | Float_t y = fClustersLay2[iC2][1] - vtx[1]; | |
125 | Float_t z = fClustersLay2[iC2][2] - vtx[2]; | |
126 | ||
127 | Float_t r = TMath::Sqrt(TMath::Power(x,2) + | |
128 | TMath::Power(y,2) + | |
129 | TMath::Power(z,2)); | |
130 | ||
131 | fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta | |
132 | fClustersLay2[iC2][1] = TMath::ATan(y/x); // Store Phi | |
133 | fClustersLay2[iC2][2] = z; // Store z | |
134 | ||
135 | // this only needs to be initialized for the fNClustersLay2 first associations | |
136 | fAssociationFlag[iC2] = kFALSE; | |
137 | } | |
138 | ||
139 | //########################################################### | |
140 | // Loop on layer 1 | |
141 | for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) { | |
142 | ||
143 | // reset of variables for multiple candidates | |
144 | Int_t iC2WithBestPhi = 0; // reset | |
145 | Float_t dPhimin = 100.; // just to put a huge number! | |
146 | ||
147 | // Loop on layer 2 | |
148 | for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) { | |
149 | ||
150 | // The following excludes double associations | |
151 | if (!fAssociationFlag[iC2]) { | |
152 | ||
153 | // find the difference in angles | |
154 | Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0]; | |
155 | Float_t dPhi = fClustersLay2[iC2][1] - fClustersLay1[iC1][1]; | |
156 | ||
157 | // find the difference in z (between linear projection from layer 1 | |
158 | // and the actual point: Dzeta= z1/r1*r2 -z2) | |
159 | Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]); | |
160 | Float_t dZeta = fClustersLay2[iC1][2]*r2 - fClustersLay2[iC2][2]; | |
161 | ||
162 | if (fHistOn) { | |
163 | fhClustersDPhi->Fill(dPhi); | |
164 | fhClustersDTheta->Fill(dTheta); | |
165 | fhClustersDZeta->Fill(dZeta); | |
166 | fhDPhiVsDThetaAll->Fill(dTheta, dPhi); | |
167 | } | |
168 | // make "elliptical" cut in Phi and Zeta! | |
169 | Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2)); | |
170 | if (d>1) continue; | |
171 | ||
172 | //look for the minimum distance in Phi: the minimum is in iC2WithBestPhi | |
173 | if (TMath::Abs(dPhi) < dPhimin) { | |
174 | dPhimin = TMath::Abs(dPhi); | |
175 | iC2WithBestPhi = iC2; | |
176 | } | |
177 | } | |
178 | } // end of loop over clusters in layer 2 | |
179 | ||
180 | if (dPhimin<100) { // This means that a cluster in layer 2 was found that mathes with iC1 | |
181 | ||
182 | if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestPhi] = kTRUE; // flag the association | |
183 | ||
184 | // store the tracklet | |
185 | ||
186 | // use the average theta from the clusters in the two layers | |
187 | fTracklets[fNTracklets][0] = 0.5*(fClustersLay1[iC1][0]+fClustersLay2[iC2WithBestPhi][0]); | |
188 | // use the phi from the clusters in the first layer | |
189 | fTracklets[fNTracklets][1] = fClustersLay1[iC1][1]; | |
190 | // Store the difference between phi1 and phi2 | |
191 | fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestPhi][1]; | |
192 | fNTracklets++; | |
193 | ||
194 | AliDebug(1,Form(" Adding tracklet candidate %d (cluster %d of layer 1 and %d of layer 2)", fNTracklets, iC1)); | |
195 | } | |
196 | } // end of loop over clusters in layer 1 | |
197 | ||
198 | AliDebug(1,Form("%d tracklets found", fNTracklets)); | |
199 | } | |
200 | ||
201 | //____________________________________________________________________ | |
202 | void | |
203 | AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) { | |
204 | // This method | |
205 | // - gets the clusters from the cluster tree | |
206 | // - convert them into global coordinates | |
207 | // - store them in the internal arrays | |
208 | ||
209 | AliDebug(1,"Loading clusters ..."); | |
210 | ||
211 | fNClustersLay1 = 0; | |
212 | fNClustersLay2 = 0; | |
213 | ||
214 | TClonesArray* itsClusters = new TClonesArray("AliITSclusterV2"); | |
215 | TBranch* itsClusterBranch=itsClusterTree->GetBranch("Clusters"); | |
216 | itsClusterBranch->SetAddress(&itsClusters); | |
217 | ||
218 | Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries(); | |
219 | ||
220 | // loop over the its subdetectors | |
221 | for (Int_t iIts=0; iIts < nItsSubs; iIts++) { | |
222 | ||
223 | if (!itsClusterTree->GetEvent(iIts)) | |
224 | continue; | |
225 | ||
226 | Int_t nClusters = itsClusters->GetEntriesFast(); | |
227 | ||
228 | // stuff needed to get the global coordinates | |
229 | Double_t rot[9]; fGeometry->GetRotMatrix(iIts,rot); | |
230 | Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det); | |
231 | Float_t tx,ty,tz; fGeometry->GetTrans(lay,lad,det,tx,ty,tz); | |
232 | ||
233 | // Below: | |
234 | // "alpha" is the angle from the global X-axis to the | |
235 | // local GEANT X'-axis ( rot[0]=cos(alpha) and rot[1]=sin(alpha) ) | |
236 | // "phi" is the angle from the global X-axis to the | |
237 | // local cluster X"-axis | |
238 | ||
239 | Double_t alpha = TMath::ATan2(rot[1],rot[0])+TMath::Pi(); | |
240 | Double_t itsPhi = TMath::Pi()/2+alpha; | |
241 | ||
242 | if (lay==1) itsPhi+=TMath::Pi(); | |
243 | Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi); | |
244 | Double_t r=tx*cp+ty*sp; | |
245 | ||
246 | // loop over clusters | |
247 | while(nClusters--) { | |
248 | AliITSclusterV2* cluster = (AliITSclusterV2*)itsClusters->UncheckedAt(nClusters); | |
249 | ||
250 | if (cluster->GetLayer()>1) | |
251 | continue; | |
252 | ||
253 | Float_t x = r*cp - cluster->GetY()*sp; | |
254 | Float_t y = r*sp + cluster->GetY()*cp; | |
255 | Float_t z = cluster->GetZ(); | |
256 | ||
257 | if (cluster->GetLayer()==0) { | |
258 | fClustersLay1[fNClustersLay1][0] = x; | |
259 | fClustersLay1[fNClustersLay1][1] = y; | |
260 | fClustersLay1[fNClustersLay1][2] = z; | |
261 | fNClustersLay1++; | |
262 | } | |
263 | if (cluster->GetLayer()==1) { | |
264 | fClustersLay2[fNClustersLay2][0] = x; | |
265 | fClustersLay2[fNClustersLay2][1] = y; | |
266 | fClustersLay2[fNClustersLay2][2] = z; | |
267 | fNClustersLay2++; | |
268 | } | |
269 | ||
270 | }// end of cluster loop | |
271 | } // end of its "subdetector" loop | |
272 | ||
273 | AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2)); | |
274 | } | |
275 | //____________________________________________________________________ | |
276 | void | |
277 | AliITSMultReconstructor::SaveHists() { | |
278 | ||
279 | if (!fHistOn) | |
280 | return; | |
281 | ||
282 | cout << "Saving histograms" << endl; | |
283 | ||
284 | fhClustersDPhi->Write(); | |
285 | fhClustersDTheta->Write(); | |
286 | fhClustersDZeta->Write(); | |
287 | fhDPhiVsDThetaAll->Write(); | |
288 | fhDPhiVsDThetaAcc->Write(); | |
289 | } |