changed sixe of fOutputContainer
[u/mrichter/AliRoot.git] / ITS / AliITSMultReconstructor.cxx
CommitLineData
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//
3ef75756 18// NOTE: The cuts on phi and zeta depends on the interacting system (p-p
19// or Pb-Pb). Please, check the file AliITSMultReconstructor.h and be
20// sure that SetPhiWindow and SetZetaWindow are defined accordingly.
ac903f1b 21//
3ef75756 22//
23//
ac903f1b 24//
25//____________________________________________________________________
26
27#include "AliITSMultReconstructor.h"
28
29#include "TTree.h"
30#include "TH1F.h"
31#include "TH2F.h"
32
b51872de 33#include "AliITSRecPoint.h"
ac903f1b 34#include "AliITSgeom.h"
35#include "AliLog.h"
36
37//____________________________________________________________________
0762f3a8 38ClassImp(AliITSMultReconstructor)
ac903f1b 39
3ef75756 40
ac903f1b 41//____________________________________________________________________
7537d03c 42AliITSMultReconstructor::AliITSMultReconstructor():
43fGeometry(0),
44fClustersLay1(0),
45fClustersLay2(0),
46fTracklets(0),
47fAssociationFlag(0),
48fNClustersLay1(0),
49fNClustersLay2(0),
50fNTracklets(0),
51fPhiWindow(0),
52fZetaWindow(0),
53fOnlyOneTrackletPerC2(0),
54fHistOn(0),
55fhClustersDPhiAcc(0),
56fhClustersDThetaAcc(0),
57fhClustersDZetaAcc(0),
58fhClustersDPhiAll(0),
59fhClustersDThetaAll(0),
60fhClustersDZetaAll(0),
61fhDPhiVsDThetaAll(0),
62fhDPhiVsDThetaAcc(0),
63fhDPhiVsDZetaAll(0),
64fhDPhiVsDZetaAcc(0),
65fhetaTracklets(0),
66fhphiTracklets(0),
67fhetaClustersLay1(0),
68fhphiClustersLay1(0){
3ef75756 69 // Method to reconstruct the charged particles multiplicity with the
70 // SPD (tracklets).
ac903f1b 71
72 fGeometry =0;
73
74 SetHistOn();
75 SetPhiWindow();
76 SetZetaWindow();
77 SetOnlyOneTrackletPerC2();
78
79 fClustersLay1 = new Float_t*[300000];
80 fClustersLay2 = new Float_t*[300000];
81 fTracklets = new Float_t*[300000];
82 fAssociationFlag = new Bool_t[300000];
83
84 for(Int_t i=0; i<300000; i++) {
85 fClustersLay1[i] = new Float_t[3];
86 fClustersLay2[i] = new Float_t[3];
87 fTracklets[i] = new Float_t[3];
88 fAssociationFlag[i] = kFALSE;
89 }
90
91 // definition of histograms
ddced3c8 92 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,-0.1,0.1);
93 fhClustersDPhiAcc->SetDirectory(0);
94 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
95 fhClustersDThetaAcc->SetDirectory(0);
96 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
97 fhClustersDZetaAcc->SetDirectory(0);
98
99 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,-0.1,0.1);
100 fhDPhiVsDZetaAcc->SetDirectory(0);
101 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
ac903f1b 102 fhDPhiVsDThetaAcc->SetDirectory(0);
103
ddced3c8 104 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,-0.5,0.5);
105 fhClustersDPhiAll->SetDirectory(0);
106 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
107 fhClustersDThetaAll->SetDirectory(0);
108 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
109 fhClustersDZetaAll->SetDirectory(0);
110
111 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,-0.5,0.5);
112 fhDPhiVsDZetaAll->SetDirectory(0);
113 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,-0.5,0.5);
114 fhDPhiVsDThetaAll->SetDirectory(0);
115
116 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
ddced3c8 117 fhphiTracklets = new TH1F("phiTracklets", "phi", 100,-3.14159,3.14159);
ddced3c8 118 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
ddced3c8 119 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100,-3.141,3.141);
3ef75756 120
ac903f1b 121}
ddced3c8 122
3ef75756 123//______________________________________________________________________
7537d03c 124AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
125fGeometry(mr.fGeometry),
126fClustersLay1(mr.fClustersLay1),
127fClustersLay2(mr.fClustersLay2),
128fTracklets(mr.fTracklets),
129fAssociationFlag(mr.fAssociationFlag),
130fNClustersLay1(mr.fNClustersLay1),
131fNClustersLay2(mr.fNClustersLay2),
132fNTracklets(mr.fNTracklets),
133fPhiWindow(mr.fPhiWindow),
134fZetaWindow(mr.fZetaWindow),
135fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
136fHistOn(mr.fHistOn),
137fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
138fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
139fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
140fhClustersDPhiAll(mr.fhClustersDPhiAll),
141fhClustersDThetaAll(mr.fhClustersDThetaAll),
142fhClustersDZetaAll(mr.fhClustersDZetaAll),
143fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
144fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
145fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
146fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
147fhetaTracklets(mr.fhetaTracklets),
148fhphiTracklets(mr.fhphiTracklets),
149fhetaClustersLay1(mr.fhetaClustersLay1),
150fhphiClustersLay1(mr.fhphiClustersLay1) {
3ef75756 151 // Copy constructor
7537d03c 152
3ef75756 153}
154
155//______________________________________________________________________
7537d03c 156AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
3ef75756 157 // Assignment operator
7537d03c 158 this->~AliITSMultReconstructor();
159 new(this) AliITSMultReconstructor(mr);
3ef75756 160 return *this;
161}
162
163//______________________________________________________________________
164AliITSMultReconstructor::~AliITSMultReconstructor(){
165 // Destructor
1ba5b31c 166
167 // delete histograms
168 delete fhClustersDPhiAcc;
169 delete fhClustersDThetaAcc;
170 delete fhClustersDZetaAcc;
171 delete fhClustersDPhiAll;
172 delete fhClustersDThetaAll;
173 delete fhClustersDZetaAll;
174 delete fhDPhiVsDThetaAll;
175 delete fhDPhiVsDThetaAcc;
176 delete fhDPhiVsDZetaAll;
177 delete fhDPhiVsDZetaAcc;
178 delete fhetaTracklets;
179 delete fhphiTracklets;
180 delete fhetaClustersLay1;
181 delete fhphiClustersLay1;
182
183 // delete arrays
184 for(Int_t i=0; i<300000; i++) {
185 delete [] fClustersLay1[i];
186 delete [] fClustersLay2[i];
187 delete [] fTracklets[i];
ddced3c8 188 }
1ba5b31c 189 delete [] fClustersLay1;
190 delete [] fClustersLay2;
191 delete [] fTracklets;
192
193 delete [] fAssociationFlag;
ddced3c8 194}
ac903f1b 195
196//____________________________________________________________________
197void
198AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
199 //
200 // - calls LoadClusterArray that finds the position of the clusters
201 // (in global coord)
202 // - convert the cluster coordinates to theta, phi (seen from the
203 // interaction vertex). The third coordinate is used for ....
204 // - makes an array of tracklets
205 //
206 // After this method has been called, the clusters of the two layers
207 // and the tracklets can be retrieved by calling the Get'er methods.
208
ac903f1b 209 // reset counters
210 fNClustersLay1 = 0;
211 fNClustersLay2 = 0;
212 fNTracklets = 0;
213
214 // loading the clusters
215 LoadClusterArrays(clusterTree);
3ef75756 216
ac903f1b 217 // find the tracklets
218 AliDebug(1,"Looking for tracklets... ");
219
220 //###########################################################
221 // Loop on layer 1 : finding theta, phi and z
222 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
223 Float_t x = fClustersLay1[iC1][0] - vtx[0];
224 Float_t y = fClustersLay1[iC1][1] - vtx[1];
225 Float_t z = fClustersLay1[iC1][2] - vtx[2];
ddced3c8 226
ac903f1b 227 Float_t r = TMath::Sqrt(TMath::Power(x,2) +
228 TMath::Power(y,2) +
229 TMath::Power(z,2));
230
231 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
ddced3c8 232 fClustersLay1[iC1][1] = TMath::ATan2(x,y); // Store Phi
ac903f1b 233 fClustersLay1[iC1][2] = z/r; // Store scaled z
ddced3c8 234 if (fHistOn) {
235 Float_t eta=fClustersLay1[iC1][0];
236 eta= TMath::Tan(eta/2.);
237 eta=-TMath::Log(eta);
238 fhetaClustersLay1->Fill(eta);
239 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
240 }
241}
ac903f1b 242
243 // Loop on layer 2 : finding theta, phi and r
244 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
245 Float_t x = fClustersLay2[iC2][0] - vtx[0];
246 Float_t y = fClustersLay2[iC2][1] - vtx[1];
247 Float_t z = fClustersLay2[iC2][2] - vtx[2];
ddced3c8 248
ac903f1b 249 Float_t r = TMath::Sqrt(TMath::Power(x,2) +
250 TMath::Power(y,2) +
251 TMath::Power(z,2));
252
253 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
ddced3c8 254 fClustersLay2[iC2][1] = TMath::ATan2(x,y); // Store Phi
ac903f1b 255 fClustersLay2[iC2][2] = z; // Store z
256
ddced3c8 257 // this only needs to be initialized for the fNClustersLay2 first associations
ac903f1b 258 fAssociationFlag[iC2] = kFALSE;
259 }
260
261 //###########################################################
262 // Loop on layer 1
263 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
264
265 // reset of variables for multiple candidates
ddced3c8 266 Int_t iC2WithBestDist = 0; // reset
3ef75756 267 Float_t distmin = 100.; // just to put a huge number!
ddced3c8 268 Float_t dPhimin = 0.; // Used for histograms only!
269 Float_t dThetamin = 0.; // Used for histograms only!
270 Float_t dZetamin = 0.; // Used for histograms only!
ac903f1b 271
272 // Loop on layer 2
273 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
274
275 // The following excludes double associations
276 if (!fAssociationFlag[iC2]) {
277
278 // find the difference in angles
279 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
280 Float_t dPhi = fClustersLay2[iC2][1] - fClustersLay1[iC1][1];
281
282 // find the difference in z (between linear projection from layer 1
283 // and the actual point: Dzeta= z1/r1*r2 -z2)
ddced3c8 284 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
285 Float_t dZeta = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
286
287 if (fHistOn) {
288 fhClustersDPhiAll->Fill(dPhi);
289 fhClustersDThetaAll->Fill(dTheta);
290 fhClustersDZetaAll->Fill(dZeta);
ac903f1b 291 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
ddced3c8 292 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
ac903f1b 293 }
294 // make "elliptical" cut in Phi and Zeta!
295 Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
3ef75756 296
ac903f1b 297 if (d>1) continue;
298
ddced3c8 299 //look for the minimum distance: the minimum is in iC2WithBestDist
3ef75756 300 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
301 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
ddced3c8 302 dPhimin = dPhi;
303 dThetamin = dTheta;
304 dZetamin = dZeta;
305 iC2WithBestDist = iC2;
ac903f1b 306 }
307 }
308 } // end of loop over clusters in layer 2
309
3ef75756 310 if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
311
312 if (fHistOn) {
313 fhClustersDPhiAcc->Fill(dPhimin);
314 fhClustersDThetaAcc->Fill(dThetamin);
315 fhClustersDZetaAcc->Fill(dZetamin);
316 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
317 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
318 }
ac903f1b 319
ddced3c8 320 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
ac903f1b 321
322 // store the tracklet
323
ddced3c8 324 // use the theta from the clusters in the first layer
325 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
ac903f1b 326 // use the phi from the clusters in the first layer
327 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
328 // Store the difference between phi1 and phi2
ddced3c8 329 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
330
3ef75756 331 if (fHistOn) {
332 Float_t eta=fTracklets[fNTracklets][0];
333 eta= TMath::Tan(eta/2.);
334 eta=-TMath::Log(eta);
335 fhetaTracklets->Fill(eta);
336 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
337 }
ac903f1b 338
3ef75756 339 AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
340 AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1,
341 iC2WithBestDist));
342 fNTracklets++;
ac903f1b 343 }
3ef75756 344
345 // Delete the following else if you do not want to save Clusters!
346
347 else { // This means that the cluster has not been associated
348
349 // store the cluster
350
351 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
352 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
353 // Store a flag. This will indicate that the "tracklet"
354 // was indeed a single cluster!
355 fTracklets[fNTracklets][2] = -999999.;
356 AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
357 fNTracklets, iC1));
358 fNTracklets++;
359 }
360
ac903f1b 361 } // end of loop over clusters in layer 1
362
363 AliDebug(1,Form("%d tracklets found", fNTracklets));
364}
365
366//____________________________________________________________________
367void
368AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
369 // This method
370 // - gets the clusters from the cluster tree
371 // - convert them into global coordinates
372 // - store them in the internal arrays
373
374 AliDebug(1,"Loading clusters ...");
375
376 fNClustersLay1 = 0;
377 fNClustersLay2 = 0;
378
b51872de 379 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
380 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
ddced3c8 381
ac903f1b 382 itsClusterBranch->SetAddress(&itsClusters);
ddced3c8 383
ac903f1b 384 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
ddced3c8 385
ac903f1b 386 // loop over the its subdetectors
387 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
388
389 if (!itsClusterTree->GetEvent(iIts))
390 continue;
391
392 Int_t nClusters = itsClusters->GetEntriesFast();
393
394 // stuff needed to get the global coordinates
395 Double_t rot[9]; fGeometry->GetRotMatrix(iIts,rot);
396 Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det);
397 Float_t tx,ty,tz; fGeometry->GetTrans(lay,lad,det,tx,ty,tz);
398
399 // Below:
400 // "alpha" is the angle from the global X-axis to the
401 // local GEANT X'-axis ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
402 // "phi" is the angle from the global X-axis to the
403 // local cluster X"-axis
404
405 Double_t alpha = TMath::ATan2(rot[1],rot[0])+TMath::Pi();
406 Double_t itsPhi = TMath::Pi()/2+alpha;
407
408 if (lay==1) itsPhi+=TMath::Pi();
409 Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi);
410 Double_t r=tx*cp+ty*sp;
411
412 // loop over clusters
413 while(nClusters--) {
b51872de 414 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
ac903f1b 415
416 if (cluster->GetLayer()>1)
417 continue;
418
419 Float_t x = r*cp - cluster->GetY()*sp;
420 Float_t y = r*sp + cluster->GetY()*cp;
421 Float_t z = cluster->GetZ();
422
423 if (cluster->GetLayer()==0) {
424 fClustersLay1[fNClustersLay1][0] = x;
425 fClustersLay1[fNClustersLay1][1] = y;
426 fClustersLay1[fNClustersLay1][2] = z;
427 fNClustersLay1++;
428 }
429 if (cluster->GetLayer()==1) {
430 fClustersLay2[fNClustersLay2][0] = x;
431 fClustersLay2[fNClustersLay2][1] = y;
432 fClustersLay2[fNClustersLay2][2] = z;
433 fNClustersLay2++;
434 }
435
436 }// end of cluster loop
437 } // end of its "subdetector" loop
438
439 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
440}
441//____________________________________________________________________
442void
443AliITSMultReconstructor::SaveHists() {
3ef75756 444 // This method save the histograms on the output file
445 // (only if fHistOn is TRUE).
ac903f1b 446
447 if (!fHistOn)
448 return;
449
ddced3c8 450 fhClustersDPhiAll->Write();
451 fhClustersDThetaAll->Write();
452 fhClustersDZetaAll->Write();
ac903f1b 453 fhDPhiVsDThetaAll->Write();
ddced3c8 454 fhDPhiVsDZetaAll->Write();
455
456 fhClustersDPhiAcc->Write();
457 fhClustersDThetaAcc->Write();
458 fhClustersDZetaAcc->Write();
ac903f1b 459 fhDPhiVsDThetaAcc->Write();
ddced3c8 460 fhDPhiVsDZetaAcc->Write();
461
462 fhetaTracklets->Write();
463 fhphiTracklets->Write();
464 fhetaClustersLay1->Write();
465 fhphiClustersLay1->Write();
ac903f1b 466}