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