]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSTrackleterSPDEff.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Overlap fixed (M. Sitta)
[u/mrichter/AliRoot.git] / ITS / AliITSTrackleterSPDEff.cxx
CommitLineData
275a301c 1/**************************************************************************
2 * Copyright(c) 2007-2009, 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 **************************************************************************/
275a301c 15//____________________________________________________________________
84161aec 16//
275a301c 17// AliITSTrackleterSPDEff - find SPD chips efficiencies by using tracklets.
84161aec 18//
19// This class has been developed from AliITSMultReconstructor (see
20// it for more details). It is the class for the Trackleter used to estimate
21// SPD plane efficiency.
275a301c 22// The trackleter prediction is built using the vertex and 1 cluster.
23//
84161aec 24//
275a301c 25// Author : Giuseppe Eugenio Bruno, based on the skeleton of Reconstruct method provided by Tiziano Virgili
26// email: giuseppe.bruno@ba.infn.it
84161aec 27//
275a301c 28//____________________________________________________________________
29
84161aec 30/* $Id$ */
31
275a301c 32#include <TFile.h>
33#include <TParticle.h>
34#include <TSystem.h>
35#include <Riostream.h>
36
37#include "AliITSMultReconstructor.h"
38#include "AliITSTrackleterSPDEff.h"
39#include "AliITSgeomTGeo.h"
40#include "AliLog.h"
41#include "AliITSPlaneEffSPD.h"
42#include "AliStack.h"
43
44//____________________________________________________________________
45ClassImp(AliITSTrackleterSPDEff)
46
47
48//____________________________________________________________________
49AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
50AliITSMultReconstructor(),
51fAssociationFlag1(0),
52fChipPredOnLay2(0),
53fChipPredOnLay1(0),
54fNTracklets1(0),
55fPhiWindowL1(0),
56fZetaWindowL1(0),
57fOnlyOneTrackletPerC1(0),
58fPlaneEffSPD(0),
59fMC(0),
60fUseOnlyPrimaryForPred(0),
61fUseOnlySecondaryForPred(0),
62fUseOnlySameParticle(0),
63fUseOnlyDifferentParticle(0),
64fUseOnlyStableParticle(0),
65fPredictionPrimary(0),
66fPredictionSecondary(0),
67fClusterPrimary(0),
68fClusterSecondary(0),
69fhClustersDPhiInterpAcc(0),
70fhClustersDThetaInterpAcc(0),
71fhClustersDZetaInterpAcc(0),
72fhClustersDPhiInterpAll(0),
73fhClustersDThetaInterpAll(0),
74fhClustersDZetaInterpAll(0),
75fhDPhiVsDThetaInterpAll(0),
76fhDPhiVsDThetaInterpAcc(0),
77fhDPhiVsDZetaInterpAll(0),
78fhDPhiVsDZetaInterpAcc(0),
79fhetaClustersLay2(0),
80fhphiClustersLay2(0)
81{
84161aec 82 // default constructor
275a301c 83
84 SetPhiWindowL1();
85 SetZetaWindowL1();
86 SetOnlyOneTrackletPerC1();
87
88 fAssociationFlag1 = new Bool_t[300000];
89 fChipPredOnLay2 = new UInt_t[300000];
90 fChipPredOnLay1 = new UInt_t[300000];
91
92 for(Int_t i=0; i<300000; i++) {
93 fAssociationFlag1[i] = kFALSE;
94 }
95
96 if (GetHistOn()) BookHistos();
97
98 fPlaneEffSPD = new AliITSPlaneEffSPD();
99}
100//______________________________________________________________________
101AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) : AliITSMultReconstructor(mr),
102fAssociationFlag1(mr.fAssociationFlag1),
103fChipPredOnLay2(mr.fChipPredOnLay2),
104fChipPredOnLay1(mr.fChipPredOnLay1),
105fNTracklets1(mr.fNTracklets1),
106fPhiWindowL1(mr.fPhiWindowL1),
107fZetaWindowL1(mr.fZetaWindowL1),
108fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
109fPlaneEffSPD(mr.fPlaneEffSPD),
110fMC(mr.fMC),
111fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
112fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
113fUseOnlySameParticle(mr.fUseOnlySameParticle),
114fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
115fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
116fPredictionPrimary(mr.fPredictionPrimary),
117fPredictionSecondary(mr.fPredictionSecondary),
118fClusterPrimary(mr.fClusterPrimary),
119fClusterSecondary(mr.fClusterSecondary),
120fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
121fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
122fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
123fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
124fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
125fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
126fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
127fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
128fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
129fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
130fhetaClustersLay2(mr.fhetaClustersLay2),
131fhphiClustersLay2(mr.fhphiClustersLay2)
132{
133 // Copy constructor
134}
135
136//______________________________________________________________________
137AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
138 // Assignment operator
139 this->~AliITSTrackleterSPDEff();
140 new(this) AliITSTrackleterSPDEff(mr);
141 return *this;
142}
143//______________________________________________________________________
144AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
145 // Destructor
146
147 // delete histograms
148 DeleteHistos();
149
150 delete [] fAssociationFlag1;
151
152 delete [] fChipPredOnLay2;
153 delete [] fChipPredOnLay1;
154
155 delete [] fPredictionPrimary;
156 delete [] fPredictionSecondary;
157 delete [] fClusterPrimary;
158 delete [] fClusterSecondary;
159
160 // delete PlaneEff
161 delete fPlaneEffSPD;
162}
163//____________________________________________________________________
164void
165AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/,AliStack *pStack) {
166 //
167 // - calls LoadClusterArray that finds the position of the clusters
168 // (in global coord)
169 // - convert the cluster coordinates to theta, phi (seen from the
170 // interaction vertex). Find the extrapolation/interpolation point.
171 // - Find the chip corresponding to that
172 // - Check if there is a cluster near that point
173 //
174
175 // reset counters
176 fNClustersLay1 = 0;
177 fNClustersLay2 = 0;
178 fNTracklets = 0;
179 fNSingleCluster = 0;
180 // loading the clusters
181 LoadClusterArrays(clusterTree);
182 if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
183 Bool_t found;
184 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
185 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
186 Int_t nfClu1=0; Int_t ntClu1=0;
187 Int_t nfClu2=0; Int_t ntClu2=0;
188
189
190 // find the tracklets
191 AliDebug(1,"Looking for tracklets... ");
192 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
193
194 //###########################################################
195 // Loop on layer 1 : finding theta, phi and z
196 UInt_t key;
197 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
198 Float_t x = fClustersLay1[iC1][0] - vtx[0];
199 Float_t y = fClustersLay1[iC1][1] - vtx[1];
200 Float_t z = fClustersLay1[iC1][2] - vtx[2];
201
202 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
203
204 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
205 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
206 fClustersLay1[iC1][2] = z; // Store z
207
208 // find the Radius and the chip corresponding to the extrapolation point
209
210 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
211 if (!found) {
212 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
213 key=999999; // also some other actions should be taken if not Found
214 }
215 nfTraPred2+=(Int_t)found; // this for debugging purpose
216 ntTraPred2++; // to check efficiency of the method FindChip
217 fChipPredOnLay2[iC1] = key;
218 fAssociationFlag1[iC1] = kFALSE;
219
220 if (fHistOn) {
221 Float_t eta=fClustersLay1[iC1][0];
222 eta= TMath::Tan(eta/2.);
223 eta=-TMath::Log(eta);
224 fhetaClustersLay1->Fill(eta);
225 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
226 }
227 }
228
229 // Loop on layer 2 : finding theta, phi and r
230 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
231 Float_t x = fClustersLay2[iC2][0] - vtx[0];
232 Float_t y = fClustersLay2[iC2][1] - vtx[1];
233 Float_t z = fClustersLay2[iC2][2] - vtx[2];
234
235 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
236
237 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
238 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
239 fClustersLay2[iC2][2] = z; // Store z
240
241 // find the Radius and the chip corresponding to the extrapolation point
242
243 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
244 if (!found) {
245 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
246 key=999999;
247 }
248 nfTraPred1+=(Int_t)found; // this for debugging purpose
249 ntTraPred1++; // to check efficiency of the method FindChip
250 fChipPredOnLay1[iC2] = key;
251 fAssociationFlag[iC2] = kFALSE;
252
253 if (fHistOn) {
254 Float_t eta=fClustersLay2[iC2][0];
255 eta= TMath::Tan(eta/2.);
256 eta=-TMath::Log(eta);
257 fhetaClustersLay2->Fill(eta);
258 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
259 }
260 }
261
262 //###########################################################
263
264 // First part : Extrapolation to Layer 2
265
266 // Loop on layer 1
267 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
268
269 // reset of variables for multiple candidates
270 Int_t iC2WithBestDist = 0; // reset
271 Float_t distmin = 100.; // just to put a huge number!
272 Float_t dPhimin = 0.; // Used for histograms only!
273 Float_t dThetamin = 0.; // Used for histograms only!
274 Float_t dZetamin = 0.; // Used for histograms only!
275
276 // in any case, if MC has been required, store statistics of primaries and secondaries
277 if (fMC) {
278 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
279 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
280 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
281 // do it always as a function of the chip number used to built the prediction
282 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
283 if (!found) {AliWarning(
284 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
285 else {
286 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
287 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
288 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
289 { // this cluster is from a primary particle
290 fClusterPrimary[key]++;
291 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
292 } else { // this cluster is from a secondary particle
293 fClusterSecondary[key]++;
294 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
295 }
296 }
297 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
298 // (in case the prediction is reliable)
299 if( fChipPredOnLay2[iC1]<1200) {
300 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
301 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
302 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
303 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
304 }
305 }
306
307 // Loop on layer 2
308 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
309
310 // The following excludes double associations
311 if (!fAssociationFlag[iC2]) {
312
313 // find the difference in angles
314 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
315 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
316 // take into account boundary condition
317 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
318
319 // find the difference in z (between linear projection from layer 1
320 // and the actual point: Dzeta= z1/r1*r2 -z2)
321 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
322 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
323
324 if (fHistOn) {
325 fhClustersDPhiAll->Fill(dPhi);
326 fhClustersDThetaAll->Fill(dTheta);
327 fhClustersDZetaAll->Fill(dZeta);
328 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
329 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
330 }
331
332 // make "elliptical" cut in Phi and Zeta!
333 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
334 dZeta*dZeta/fZetaWindow/fZetaWindow);
335
336 if (d>1) continue;
337
338 //look for the minimum distance: the minimum is in iC2WithBestDist
339 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
340 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
341 dPhimin = dPhi;
342 dThetamin = dTheta;
343 dZetamin = dZeta;
344 iC2WithBestDist = iC2;
345 }
346 }
347 } // end of loop over clusters in layer 2
348
349 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
350
351 if (fHistOn) {
352 fhClustersDPhiAcc->Fill(dPhimin);
353 fhClustersDThetaAcc->Fill(dThetamin);
354 fhClustersDZetaAcc->Fill(dZetamin);
355 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
356 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
357 }
358
359 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
360 // flag the association
361
362 // store the tracklet
363
364 // use the theta from the clusters in the first layer
365 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
366 // use the phi from the clusters in the first layer
367 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
368 // Store the difference between phi1 and phi2
369 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
370
371 // find labels
372 Int_t label1 = 0;
373 Int_t label2 = 0;
374 while (label2 < 3)
375 {
376 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
377 break;
378 label1++;
379 if (label1 == 3)
380 {
381 label1 = 0;
382 label2++;
383 }
384 }
385
386 if (label2 < 3)
387 {
388 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
389 }
390 else
391 {
392 fTracklets[fNTracklets][3] = -2;
393 }
394
395 if (fHistOn) {
396 Float_t eta=fTracklets[fNTracklets][0];
397 eta= TMath::Tan(eta/2.);
398 eta=-TMath::Log(eta);
399 fhetaTracklets->Fill(eta);
400 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
401 }
402
403// Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
404 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
405 if(!found){
406 AliWarning(
407 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
408 key=999999;
409 }
410 nfClu2+=(Int_t)found; // this for debugging purpose
411 ntClu2++; // to check efficiency of the method FindChip
412 if(key<1200) { // the Chip has been found
413 if(fMC) { // this part only for MC
414 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
415 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
416 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
417 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
418 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
419 }
420
421 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
422 // OK, success
423 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
424 }
425 else {
426 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
427 // (might be in the tracking tollerance)
428 }
429 }
430
431 fNTracklets++;
432
433 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
434 else if (fChipPredOnLay2[iC1]<1200) fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
435
436 } // end of loop over clusters in layer 1
437
438 fNTracklets1=fNTracklets;
439
440//###################################################################
441
442 // Second part : Interpolation to Layer 1
443
444 // Loop on layer 2
445 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
446
447 // reset of variables for multiple candidates
448 Int_t iC1WithBestDist = 0; // reset
449 Float_t distmin = 100.; // just to put a huge number!
450 Float_t dPhimin = 0.; // Used for histograms only!
451 Float_t dThetamin = 0.; // Used for histograms only!
452 Float_t dZetamin = 0.; // Used for histograms only!
453
454 // in any case, if MC has been required, store statistics of primaries and secondaries
455 if (fMC) {
456 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
457 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
458 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
459 // do it always as a function of the chip number used to built the prediction
460 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
461 if (!found) {AliWarning(
462 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
463 else {
464 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
465 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
466 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
467 { // this cluster is from a primary particle
468 fClusterPrimary[key]++;
469 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
470 } else { // this cluster is from a secondary particle
471 fClusterSecondary[key]++;
472 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
473 }
474 }
475 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
476 // (in case the prediction is reliable)
477 if( fChipPredOnLay1[iC2]<1200) {
478 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
479 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
480 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
481 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
482 }
483 }
484
485 // Loop on layer 1
486 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
487
488 // The following excludes double associations
489 if (!fAssociationFlag1[iC1]) {
490
491 // find the difference in angles
492 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
493 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
494 // take into account boundary condition
495 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
496
497
498 // find the difference in z (between linear projection from layer 2
499 // and the actual point: Dzeta= z2/r2*r1 -z1)
500 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
501 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
502
503
504 if (fHistOn) {
505 fhClustersDPhiInterpAll->Fill(dPhi);
506 fhClustersDThetaInterpAll->Fill(dTheta);
507 fhClustersDZetaInterpAll->Fill(dZeta);
508 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
509 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
510 }
511 // make "elliptical" cut in Phi and Zeta!
512 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
513 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
514
515 if (d>1) continue;
516
517 //look for the minimum distance: the minimum is in iC1WithBestDist
518 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
519 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
520 dPhimin = dPhi;
521 dThetamin = dTheta;
522 dZetamin = dZeta;
523 iC1WithBestDist = iC1;
524 }
525 }
526 } // end of loop over clusters in layer 1
527
528 if (distmin<100) { // This means that a cluster in layer 1 was found that mathes with iC2
529
530 if (fHistOn) {
531 fhClustersDPhiInterpAcc->Fill(dPhimin);
532 fhClustersDThetaInterpAcc->Fill(dThetamin);
533 fhClustersDZetaInterpAcc->Fill(dZetamin);
534 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
535 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
536 }
537
538 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
539 // flag the association
540
541 // store the tracklet
542
543 // use the theta from the clusters in the first layer
544 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
545 // use the phi from the clusters in the first layer
546 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
547 // Store the difference between phi1 and phi2
548 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
549
550 // find labels
551 Int_t label1 = 0;
552 Int_t label2 = 0;
553 while (label2 < 3)
554 {
555 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
556 break;
557 label1++;
558 if (label1 == 3)
559 {
560 label1 = 0;
561 label2++;
562 }
563 }
564
565 if (label2 < 3)
566 {
567 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
568 }
569 else
570 {
571 fTracklets[fNTracklets][3] = -2;
572 }
573
574// Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
575 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
576 if(!found){
577 AliWarning(
578 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
579 key=999999;
580 }
581 nfClu1+=(Int_t)found; // this for debugging purpose
582 ntClu1++; // to check efficiency of the method FindChip
583 if(key<1200) {
584 if(fMC) { // this part only for MC
585 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
586 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
587 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
588 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
589 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
590 }
591
592 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
593 // OK, success
594 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
595 } else {
596 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
597 // (might be in the tracking tollerance)
598 }
599 }
600
601 fNTracklets++;
602
603 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
604 else if (fChipPredOnLay1[iC2]<1200) fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
605
606 } // end of loop over clusters in layer 2
607
608 AliDebug(1,Form("%d tracklets found", fNTracklets));
609 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
610 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
611 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
612 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
613}
614//____________________________________________________________________
615Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
616 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
617//
618// Input: a) layer number in the range [0,1]
619// b) vtx[3]: actual vertex
620// c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
621// d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
622// e) phiVtx /
623// Output: Unique key to locate a chip
624// return: kTRUE if succesfull
625
626 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
627 Double_t r=GetRLayer(layer);
628 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
629
630 // set phiVtx in the range [0,2pi]
631 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
632
633 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
634 // of the intersection of the tracklet with the pixel layer.
635 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
636 else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
637 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
638 Double_t vtxy[2]={vtx[0],vtx[1]};
639 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
640 // this method gives you two interceptions
641 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
642 // set phiAbs in the range [0,2pi]
643 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
644 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
645 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
646 // taking the closest one to phiVtx
647 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
648 } else phiAbs=phiVtx;
649 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
650
651 // now you need to locate the chip within the idet detector,
652 // starting from the local coordinates in such a detector
653
654 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
655 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
656 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
657
658 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
659 return kTRUE;
660}
661//______________________________________________________________________________
662Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
663 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
664 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
665
666 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
667 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
668 Double_t r=TMath::Sqrt(x*x + y*y);
669
670 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
671 r += TMath::Sqrt(x*x + y*y);
672 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
673 r += TMath::Sqrt(x*x + y*y);
674 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
675 r += TMath::Sqrt(x*x + y*y);
676 r*=0.25;
677 return r;
678}
679//______________________________________________________________________________
680Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
681 Float_t &xloc, Float_t &zloc) {
682 // this method transform Global Cilindrical coordinates into local (i.e. module)
683 // cartesian coordinates
684 //
685 //Compute Cartesian Global Coordinate
686 Double_t xyzGlob[3],xyzLoc[3];
687 xyzGlob[2]=z;
688 xyzGlob[0]=r*TMath::Cos(phi);
689 xyzGlob[1]=r*TMath::Sin(phi);
690
691 xloc=0.;
692 zloc=0.;
693
694 if(idet<0) return kFALSE;
695
696 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
697 Int_t lad = Int_t(idet/ndet) + 1;
698 Int_t det = idet - (lad-1)*ndet + 1;
699
700 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
701
702 xloc = (Float_t)xyzLoc[0];
703 zloc = (Float_t)xyzLoc[2];
704
705return kTRUE;
706}
707//______________________________________________________________________________
708Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
709 //--------------------------------------------------------------------
710 //This function finds the detector crossed by the track
711 //--------------------------------------------------------------------
712 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
713 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
714 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
715 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
716
717 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
718 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
719 Double_t phiOffset=TMath::ATan2(y,x);
720 Double_t zOffset=z2;
721
722 Double_t dphi;
723 if (zOffset<0) // old geometry
724 dphi = -(phi-phiOffset);
725 else // new geometry
726 dphi = phi-phiOffset;
727
728 if (dphi < 0) dphi += 2*TMath::Pi();
729 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
730 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
731 if (np>=nladders) np-=nladders;
732 if (np<0) np+=nladders;
733
734 Double_t dz=zOffset-z;
735 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
736 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
737 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
738 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
739
740 return np*ndetectors + nz;
741}
742//____________________________________________________________
743Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
744// this method find the intersection in xy between a tracklet (straight line) and
745// a circonference (r=R), using polar coordinates.
746/*
747Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
748 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
749 - R: radius of the circle
750Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
751
752Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
753r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
754In the same system, the equation of a semi-line is: phi=phiVtx;
755Hence you get one interception only: P=(r,phiVtx)
756Finally you want P in the ABSOLUTE ALICE system.
757*/
758Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
759Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
760Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
761Double_t cC=rO*rO-R*R;
762Double_t dDelta=bB*bB-4*cC;
763if(dDelta<0) return kFALSE;
764Double_t r1,r2;
765r1=(-bB-TMath::Sqrt(dDelta))/2;
766r2=(-bB+TMath::Sqrt(dDelta))/2;
767if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
768Double_t r=TMath::Max(r1,r2); // take the positive
769Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
770Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
771pvtx[0]=r*TMath::Cos(phiVtx);
772pvtx[1]=r*TMath::Sin(phiVtx);
773pP[0]=vtx[0]+pvtx[0];
774pP[1]=vtx[1]+pvtx[1];
775phi=TMath::ATan2(pP[1],pP[0]);
776return kTRUE;
777}
778//___________________________________________________________
779Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
780while(angle >=2*TMath::Pi() || angle<0) {
781 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
782 if(angle < 0) angle+=2*TMath::Pi();
783}
784return kTRUE;
785}
786//___________________________________________________________
787Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
788if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
789if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
790if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
791// return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
792 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
793 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
794 TParticle* part = stack->Particle(ipart);
795 TParticle* part0 = stack->Particle(0); // first primary
796 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
797 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
798 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
799
800 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
801 TParticlePDG* pdgPart = part->GetPDG();
802 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
803
804 Double_t distx = part->Vx() - part0->Vx();
805 Double_t disty = part->Vy() - part0->Vy();
806 Double_t distz = part->Vz() - part0->Vz();
807 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
808
809 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
810 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
811 return kFALSE; }// primary if within 500 microns from true Vertex
812
813 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)<2) return kFALSE;
814 return kTRUE;
815}
816//_____________________________________________________________________________________________
817Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
818if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
819if(!stack) {AliError("null pointer to MC stack"); return 0;}
820if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
821
822TParticle* part = stack->Particle(ipart);
823//TParticle* part0 = stack->Particle(0); // first primary
824
825 Int_t nret=0;
826 TParticle* dau = 0;
827 Int_t nDau = 0;
828 Int_t firstDau = part->GetFirstDaughter();
829 if (firstDau > 0) {
830 Int_t lastDau = part->GetLastDaughter();
831 nDau = lastDau - firstDau + 1;
832 //printf("number of daugthers %d \n",nDau);
833 if (nDau > 0) {
834 //for(Int_t j=firstDau; j<=lastDau; j++)
835 for(Int_t j=firstDau; j<=firstDau; j++)
836 { // only first one
837 dau = stack->Particle(j);
838 Double_t distx = dau->Vx()-part->Vx();
839 Double_t disty = dau->Vy()-part->Vy();
840 Double_t distz = dau->Vz()-part->Vz();
841 Double_t distR = TMath::Sqrt(distx*distx+disty*disty+distz*distz);
842 if (distR > GetRLayer(0)+0.5) nret=1; // decay after first pixel layer
843 if (distR > GetRLayer(1)+0.5) nret=2; // decay after second pixel layer
844 }
845 }
846 } else nret = 3; // stable particle
847return nret;
848}
849//_________________________________________________________________
850void AliITSTrackleterSPDEff::InitPredictionMC() {
851if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
852fPredictionPrimary = new Int_t[1200];
853fPredictionSecondary = new Int_t[1200];
854fClusterPrimary = new Int_t[1200];
855fClusterSecondary = new Int_t[1200];
856for(Int_t i=0; i<1200; i++) {
857 fPredictionPrimary[i]=0;
858 fPredictionSecondary[i]=0;
859 fPredictionSecondary[i]=0;
860 fClusterSecondary[i]=0;
861}
862return;
863}
864//______________________________________________________________________
865Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
84161aec 866//
867// This method return the Data menmber fPredictionPrimary [1200].
868// You can call it only for MC events.
869// fPredictionPrimary[key] contains the number of tracklet predictions on the
870// given chip key built using a cluster on the other layer produced (at least)
871// from a primary particle.
872// Key refers to the chip crossed by the prediction
873//
874//
275a301c 875if (!fMC) {CallWarningMC(); return 0;}
876if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
877return fPredictionPrimary[(Int_t)key];
878}
879//______________________________________________________________________
880Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
84161aec 881//
882// This method return the Data menmber fPredictionSecondary [1200].
883// You can call it only for MC events.
884// fPredictionSecondary[key] contains the number of tracklet predictions on the
885// given chip key built using a cluster on the other layer produced (only)
886// from a secondary particle
887// Key refers to the chip crossed by the prediction
888//
889//
275a301c 890if (!fMC) {CallWarningMC(); return 0;}
891if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
892return fPredictionSecondary[(Int_t)key];
893}
894//______________________________________________________________________
895Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
84161aec 896//
897// This method return the Data menmber fClusterPrimary [1200].
898// You can call it only for MC events.
899// fClusterPrimary[key] contains the number of tracklet predictions
900// built using a cluster on that layer produced (only)
901// from a primary particle
902// Key refers to the chip used to build the prediction
903//
904//
275a301c 905if (!fMC) {CallWarningMC(); return 0;}
906if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
907return fClusterPrimary[(Int_t)key];
908}
909//______________________________________________________________________
910Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
84161aec 911//
912// This method return the Data menmber fClusterSecondary [1200].
913// You can call it only for MC events.
914// fClusterSecondary[key] contains the number of tracklet predictions
915// built using a cluster on that layer produced (only)
916// from a secondary particle
917// Key refers to the chip used to build the prediction
918//
275a301c 919if (!fMC) {CallWarningMC(); return 0;}
920if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
921return fClusterSecondary[(Int_t)key];
922}
923//______________________________________________________________________
924void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
925 // Print out some class data values in Ascii Form to output stream
926 // Inputs:
927 // ostream *os Output stream where Ascii data is to be writen
928 // Outputs:
929 // none.
930 // Return:
931 // none.
932 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow ;
933 *os << " " << fMC;
934 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
935 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
936 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
937 << " " << fUseOnlyStableParticle ;
938 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
939 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
940 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
941 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
942 return;
943}
944//______________________________________________________________________
945void AliITSTrackleterSPDEff::ReadAscii(istream *is){
946 // Read in some class data values in Ascii Form to output stream
947 // Inputs:
948 // istream *is Input stream where Ascii data is to be read in from
949 // Outputs:
950 // none.
951 // Return:
952 // none.
953
954 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow;
955 *is >> fMC;
956 if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
957 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
958 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
959 >> fUseOnlyStableParticle;
960 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
961 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
962 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
963 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
964 return;
965}
966//______________________________________________________________________
967ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
968 // Standard output streaming function
969 // Inputs:
970 // ostream &os output steam
971 // AliITSTrackleterSPDEff &s class to be streamed.
972 // Output:
973 // none.
974 // Return:
975 // ostream &os The stream pointer
976
977 s.PrintAscii(&os);
978 return os;
979}
980//______________________________________________________________________
981istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
982 // Standard inputput streaming function
983 // Inputs:
984 // istream &is input steam
985 // AliITSTrackleterSPDEff &s class to be streamed.
986 // Output:
987 // none.
988 // Return:
989 // ostream &os The stream pointer
990
991 //printf("prova %d \n", (Int_t)s.GetMC());
992 s.ReadAscii(&is);
993 return is;
994}
995//______________________________________________________________________
996void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
84161aec 997//
998// This Method write into an asci file (do not know why binary does not work)
999// the used cuts and the statistics of the MC related quantities
1000// The method SetMC() has to be called before
1001// Input TString filename: name of file for output (it deletes already existing
1002// file)
1003// Output: none
1004//
1005//
275a301c 1006 if(!fMC) {CallWarningMC(); return;}
1007 ofstream out(filename.Data(),ios::out | ios::binary);
1008 out << *this;
1009 out.close();
1010return;
1011}
1012//____________________________________________________________________
1013void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
84161aec 1014//
1015// This Method read from an asci file (do not know why binary does not work)
1016// the cuts to be used and the statistics of the MC related quantities
1017// Input TString filename: name of input file for output
1018// The method SetMC() has to be called before
1019// Output: none
1020//
1021//
275a301c 1022 if(!fMC) {CallWarningMC(); return;}
1023 if( gSystem->AccessPathName( filename.Data() ) ) {
1024 AliError( Form( "file (%s) not found", filename.Data() ) );
1025 return;
1026 }
1027
1028 ifstream in(filename.Data(),ios::in | ios::binary);
1029 in >> *this;
1030 in.close();
1031 return;
1032}
1033//____________________________________________________________________
1034Bool_t AliITSTrackleterSPDEff::SaveHists() {
84161aec 1035 // This (private) method save the histograms on the output file
275a301c 1036 // (only if fHistOn is TRUE).
84161aec 1037 // Also the histograms from the base class are saved through the
1038 // AliITSMultReconstructor::SaveHists() call
275a301c 1039
1040 if (!GetHistOn()) return kFALSE;
1041
1042 AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1043
1044 fhClustersDPhiInterpAll->Write();
1045 fhClustersDThetaInterpAll->Write();
1046 fhClustersDZetaInterpAll->Write();
1047 fhDPhiVsDThetaInterpAll->Write();
1048 fhDPhiVsDZetaInterpAll->Write();
1049
1050 fhClustersDPhiInterpAcc->Write();
1051 fhClustersDThetaInterpAcc->Write();
1052 fhClustersDZetaInterpAcc->Write();
1053 fhDPhiVsDThetaInterpAcc->Write();
1054 fhDPhiVsDZetaInterpAcc->Write();
1055
1056 fhetaClustersLay2->Write();
1057 fhphiClustersLay2->Write();
1058 return kTRUE;
1059}
1060//__________________________________________________________
1061Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1062 //
1063 // Saves the histograms into a tree and saves the trees into a file
84161aec 1064 // Also the histograms from the base class are saved
275a301c 1065 //
1066 if (!GetHistOn()) return kFALSE;
1067 if (filename.Data()=="") {
1068 AliWarning("WriteHistosToFile: null output filename!");
1069 return kFALSE;
1070 }
1071 TFile *hFile=new TFile(filename.Data(),option,
1072 "The File containing the histos for SPD efficiency studies with tracklets");
1073 if(!SaveHists()) return kFALSE;
1074 hFile->Write();
1075 hFile->Close();
1076 return kTRUE;
1077}
1078//____________________________________________________________
1079void AliITSTrackleterSPDEff::BookHistos() {
84161aec 1080//
1081// This method books addtitional histograms
1082// w.r.t. those of the base class.
1083// In particular, the differences of cluster coordinate between the two SPD
1084// layers are computed in the interpolation phase
1085//
275a301c 1086 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1087 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1088 fhClustersDPhiInterpAcc->SetDirectory(0);
1089 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1090 fhClustersDThetaInterpAcc->SetDirectory(0);
1091 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1092 fhClustersDZetaInterpAcc->SetDirectory(0);
1093
1094 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1095 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1096 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1097 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1098
1099 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1100 fhClustersDPhiInterpAll->SetDirectory(0);
1101 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1102 fhClustersDThetaInterpAll->SetDirectory(0);
1103 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1104 fhClustersDZetaInterpAll->SetDirectory(0);
1105
1106 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1107 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1108 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1109 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1110
1111 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1112 fhetaClustersLay2->SetDirectory(0);
1113 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1114 fhphiClustersLay2->SetDirectory(0);
1115 return;
1116}
1117//____________________________________________________________
1118void AliITSTrackleterSPDEff::DeleteHistos() {
84161aec 1119//
1120// Private method to delete Histograms from memory
1121// it is called. e.g., by the destructor.
1122//
275a301c 1123 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1124 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1125 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1126 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1127 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1128 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1129 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1130 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1131 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1132 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1133 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1134 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1135}
1136//_______________________________________________________________