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