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