]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ITS/AliITSTrackleterSPDEff.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Copy and =operator implemented
[u/mrichter/AliRoot.git] / ITS / AliITSTrackleterSPDEff.cxx
... / ...
CommitLineData
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//
17// AliITSTrackleterSPDEff - find SPD chips efficiencies by using tracklets.
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.
22// The trackleter prediction is built using the vertex and 1 cluster.
23//
24//
25// Author : Giuseppe Eugenio Bruno, based on the skeleton of Reconstruct method provided by Tiziano Virgili
26// email: giuseppe.bruno@ba.infn.it
27//
28//____________________________________________________________________
29
30/* $Id$ */
31
32#include <TFile.h>
33#include <TTree.h>
34#include <TParticle.h>
35#include <TSystem.h>
36#include <Riostream.h>
37#include <TClonesArray.h>
38
39#include "AliITSMultReconstructor.h"
40#include "AliITSTrackleterSPDEff.h"
41#include "AliITSgeomTGeo.h"
42#include "AliLog.h"
43#include "AliITSPlaneEffSPD.h"
44#include "AliStack.h"
45#include "AliTrackReference.h"
46
47//____________________________________________________________________
48ClassImp(AliITSTrackleterSPDEff)
49
50
51//____________________________________________________________________
52AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
53AliITSMultReconstructor(),
54fAssociationFlag1(0),
55fChipPredOnLay2(0),
56fChipPredOnLay1(0),
57fNTracklets1(0),
58fPhiWindowL1(0),
59fZetaWindowL1(0),
60fOnlyOneTrackletPerC1(0),
61fUpdateOncePerEventPlaneEff(0),
62fChipUpdatedInEvent(0),
63fPlaneEffSPD(0),
64fMC(0),
65fUseOnlyPrimaryForPred(0),
66fUseOnlySecondaryForPred(0),
67fUseOnlySameParticle(0),
68fUseOnlyDifferentParticle(0),
69fUseOnlyStableParticle(0),
70fPredictionPrimary(0),
71fPredictionSecondary(0),
72fClusterPrimary(0),
73fClusterSecondary(0),
74fSuccessPP(0),
75fSuccessTT(0),
76fSuccessS(0),
77fSuccessP(0),
78fFailureS(0),
79fFailureP(0),
80fRecons(0),
81fNonRecons(0),
82fhClustersDPhiInterpAcc(0),
83fhClustersDThetaInterpAcc(0),
84fhClustersDZetaInterpAcc(0),
85fhClustersDPhiInterpAll(0),
86fhClustersDThetaInterpAll(0),
87fhClustersDZetaInterpAll(0),
88fhDPhiVsDThetaInterpAll(0),
89fhDPhiVsDThetaInterpAcc(0),
90fhDPhiVsDZetaInterpAll(0),
91fhDPhiVsDZetaInterpAcc(0),
92fhetaClustersLay2(0),
93fhphiClustersLay2(0)
94{
95 // default constructor
96
97 SetPhiWindowL1();
98 SetZetaWindowL1();
99 SetOnlyOneTrackletPerC1();
100
101 fAssociationFlag1 = new Bool_t[300000];
102 fChipPredOnLay2 = new UInt_t[300000];
103 fChipPredOnLay1 = new UInt_t[300000];
104 fChipUpdatedInEvent = new Bool_t[1200];
105
106 for(Int_t i=0; i<300000; i++) {
107 fAssociationFlag1[i] = kFALSE;
108 }
109 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
110
111 if (GetHistOn()) BookHistos();
112
113 fPlaneEffSPD = new AliITSPlaneEffSPD();
114}
115//______________________________________________________________________
116AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) : AliITSMultReconstructor(mr),
117fAssociationFlag1(mr.fAssociationFlag1),
118fChipPredOnLay2(mr.fChipPredOnLay2),
119fChipPredOnLay1(mr.fChipPredOnLay1),
120fNTracklets1(mr.fNTracklets1),
121fPhiWindowL1(mr.fPhiWindowL1),
122fZetaWindowL1(mr.fZetaWindowL1),
123fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
124fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
125fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
126fPlaneEffSPD(mr.fPlaneEffSPD),
127fMC(mr.fMC),
128fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
129fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
130fUseOnlySameParticle(mr.fUseOnlySameParticle),
131fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
132fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
133fPredictionPrimary(mr.fPredictionPrimary),
134fPredictionSecondary(mr.fPredictionSecondary),
135fClusterPrimary(mr.fClusterPrimary),
136fClusterSecondary(mr.fClusterSecondary),
137fSuccessPP(mr.fSuccessPP),
138fSuccessTT(mr.fSuccessTT),
139fSuccessS(mr.fSuccessS),
140fSuccessP(mr.fSuccessP),
141fFailureS(mr.fFailureS),
142fFailureP(mr.fFailureP),
143fRecons(mr.fRecons),
144fNonRecons(mr.fNonRecons),
145fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
146fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
147fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
148fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
149fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
150fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
151fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
152fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
153fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
154fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
155fhetaClustersLay2(mr.fhetaClustersLay2),
156fhphiClustersLay2(mr.fhphiClustersLay2)
157{
158 // Copy constructor
159}
160
161//______________________________________________________________________
162AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
163 // Assignment operator
164 this->~AliITSTrackleterSPDEff();
165 new(this) AliITSTrackleterSPDEff(mr);
166 return *this;
167}
168//______________________________________________________________________
169AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
170 // Destructor
171
172 // delete histograms
173 DeleteHistos();
174
175 delete [] fAssociationFlag1;
176
177 delete [] fChipPredOnLay2;
178 delete [] fChipPredOnLay1;
179
180 delete [] fChipUpdatedInEvent;
181
182 delete [] fPredictionPrimary;
183 delete [] fPredictionSecondary;
184 delete [] fClusterPrimary;
185 delete [] fClusterSecondary;
186 delete [] fSuccessPP;
187 delete [] fSuccessTT;
188 delete [] fSuccessS;
189 delete [] fSuccessP;
190 delete [] fFailureS;
191 delete [] fFailureP;
192 delete [] fRecons;
193 delete [] fNonRecons;
194
195 // delete PlaneEff
196 delete fPlaneEffSPD;
197}
198//____________________________________________________________________
199void
200AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t*, AliStack *pStack, TTree *tRef) {
201 //
202 // - calls LoadClusterArray that finds the position of the clusters
203 // (in global coord)
204 // - convert the cluster coordinates to theta, phi (seen from the
205 // interaction vertex). Find the extrapolation/interpolation point.
206 // - Find the chip corresponding to that
207 // - Check if there is a cluster near that point
208 //
209
210 // reset counters
211 fNClustersLay1 = 0;
212 fNClustersLay2 = 0;
213 fNTracklets = 0;
214 fNSingleCluster = 0;
215 // loading the clusters
216 LoadClusterArrays(clusterTree);
217 if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
218 if(fMC && !tRef) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
219 Bool_t found;
220 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
221 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
222 Int_t nfClu1=0; Int_t ntClu1=0;
223 Int_t nfClu2=0; Int_t ntClu2=0;
224
225 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
226 // for this new event)
227 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
228
229 // find the tracklets
230 AliDebug(1,"Looking for tracklets... ");
231 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
232
233 //###########################################################
234 // Loop on layer 1 : finding theta, phi and z
235 UInt_t key;
236 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
237 Float_t x = fClustersLay1[iC1][0] - vtx[0];
238 Float_t y = fClustersLay1[iC1][1] - vtx[1];
239 Float_t z = fClustersLay1[iC1][2] - vtx[2];
240
241 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
242
243 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
244 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
245 fClustersLay1[iC1][2] = z; // Store z
246
247 // find the Radius and the chip corresponding to the extrapolation point
248
249 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
250 if (!found) {
251 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
252 key=999999; // also some other actions should be taken if not Found
253 }
254 nfTraPred2+=(Int_t)found; // this for debugging purpose
255 ntTraPred2++; // to check efficiency of the method FindChip
256 fChipPredOnLay2[iC1] = key;
257 fAssociationFlag1[iC1] = kFALSE;
258
259 if (fHistOn) {
260 Float_t eta=fClustersLay1[iC1][0];
261 eta= TMath::Tan(eta/2.);
262 eta=-TMath::Log(eta);
263 fhetaClustersLay1->Fill(eta);
264 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
265 }
266 }
267 // Loop on layer 2 : finding theta, phi and r
268 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
269 Float_t x = fClustersLay2[iC2][0] - vtx[0];
270 Float_t y = fClustersLay2[iC2][1] - vtx[1];
271 Float_t z = fClustersLay2[iC2][2] - vtx[2];
272
273 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
274
275 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
276 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
277 fClustersLay2[iC2][2] = z; // Store z
278
279 // find the Radius and the chip corresponding to the extrapolation point
280
281 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
282 if (!found) {
283 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
284 key=999999;
285 }
286 nfTraPred1+=(Int_t)found; // this for debugging purpose
287 ntTraPred1++; // to check efficiency of the method FindChip
288 fChipPredOnLay1[iC2] = key;
289 fAssociationFlag[iC2] = kFALSE;
290
291 if (fHistOn) {
292 Float_t eta=fClustersLay2[iC2][0];
293 eta= TMath::Tan(eta/2.);
294 eta=-TMath::Log(eta);
295 fhetaClustersLay2->Fill(eta);
296 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
297 }
298 }
299
300 //###########################################################
301
302 // First part : Extrapolation to Layer 2
303
304 // Loop on layer 1
305 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
306
307 // here the control to check whether the efficiency of the chip traversed by this tracklet
308 // prediction has already been updated in this event using another tracklet prediction
309 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
310
311 // reset of variables for multiple candidates
312 Int_t iC2WithBestDist = 0; // reset
313 Float_t distmin = 100.; // just to put a huge number!
314 Float_t dPhimin = 0.; // Used for histograms only!
315 Float_t dThetamin = 0.; // Used for histograms only!
316 Float_t dZetamin = 0.; // Used for histograms only!
317
318 // in any case, if MC has been required, store statistics of primaries and secondaries
319 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
320 if (fMC) {
321 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
322 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
323 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
324 // do it always as a function of the chip number used to built the prediction
325 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
326 if (!found) {AliWarning(
327 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
328 else {
329 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
330 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
331 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
332 { // this cluster is from a primary particle
333 fClusterPrimary[key]++;
334 primary=kTRUE;
335 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
336 } else { // this cluster is from a secondary particle
337 fClusterSecondary[key]++;
338 secondary=kTRUE;
339 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
340 }
341 }
342 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
343 // (in case the prediction is reliable)
344 if( fChipPredOnLay2[iC1]<1200) {
345 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
346 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
347 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
348 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
349 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
350 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
351 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
352 else fNonRecons[fChipPredOnLay2[iC1]]++;
353 }
354 }
355
356 // Loop on layer 2
357 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
358
359 // The following excludes double associations
360 if (!fAssociationFlag[iC2]) {
361
362 // find the difference in angles
363 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
364 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
365 // take into account boundary condition
366 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
367
368 // find the difference in z (between linear projection from layer 1
369 // and the actual point: Dzeta= z1/r1*r2 -z2)
370 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
371 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
372
373 if (fHistOn) {
374 fhClustersDPhiAll->Fill(dPhi);
375 fhClustersDThetaAll->Fill(dTheta);
376 fhClustersDZetaAll->Fill(dZeta);
377 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
378 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
379 }
380
381 // make "elliptical" cut in Phi and Zeta!
382 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
383 dZeta*dZeta/fZetaWindow/fZetaWindow);
384
385 if (d>1) continue;
386
387 //look for the minimum distance: the minimum is in iC2WithBestDist
388 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
389 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
390 dPhimin = dPhi;
391 dThetamin = dTheta;
392 dZetamin = dZeta;
393 iC2WithBestDist = iC2;
394 }
395 }
396 } // end of loop over clusters in layer 2
397
398 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
399
400 if (fHistOn) {
401 fhClustersDPhiAcc->Fill(dPhimin);
402 fhClustersDThetaAcc->Fill(dThetamin);
403 fhClustersDZetaAcc->Fill(dZetamin);
404 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
405 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
406 }
407
408 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
409 // flag the association
410
411 // store the tracklet
412
413 // use the theta from the clusters in the first layer
414 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
415 // use the phi from the clusters in the first layer
416 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
417 // Store the difference between phi1 and phi2
418 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
419
420 // find labels
421 Int_t label1 = 0;
422 Int_t label2 = 0;
423 while (label2 < 3)
424 {
425 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
426 break;
427 label1++;
428 if (label1 == 3)
429 {
430 label1 = 0;
431 label2++;
432 }
433 }
434
435 if (label2 < 3)
436 {
437 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
438 }
439 else
440 {
441 fTracklets[fNTracklets][3] = -2;
442 }
443
444 if (fHistOn) {
445 Float_t eta=fTracklets[fNTracklets][0];
446 eta= TMath::Tan(eta/2.);
447 eta=-TMath::Log(eta);
448 fhetaTracklets->Fill(eta);
449 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
450 }
451
452// Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
453 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
454 if(!found){
455 AliWarning(
456 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
457 key=999999;
458 }
459 nfClu2+=(Int_t)found; // this for debugging purpose
460 ntClu2++; // to check efficiency of the method FindChip
461 if(key<1200) { // the Chip has been found
462 if(fMC) { // this part only for MC
463 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
464 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
465 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
466 if (label2 < 3) {
467 fSuccessTT[key]++;
468 if(primary) fSuccessPP[key]++;
469 }
470 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
471 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
472 }
473
474 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
475 // OK, success
476 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
477 fChipUpdatedInEvent[key]=kTRUE;
478 if(fMC) {
479 if(primary) fSuccessP[key]++;
480 if(secondary) fSuccessS[key]++;
481 }
482 }
483 else {
484 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
485 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
486 if(fMC) {
487 if(primary) fSuccessP[key]++;
488 if(secondary) fSuccessS[key]++;
489 }
490 }
491 }
492
493 fNTracklets++;
494
495 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
496 else if (fChipPredOnLay2[iC1]<1200) {
497 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
498 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
499 if(fMC) {
500 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
501 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
502 }
503 }
504 } // end of loop over clusters in layer 1
505
506 fNTracklets1=fNTracklets;
507
508//###################################################################
509
510 // Second part : Interpolation to Layer 1
511
512 // Loop on layer 2
513 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
514
515 // here the control to check whether the efficiency of the chip traversed by this tracklet
516 // prediction has already been updated in this event using another tracklet prediction
517 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
518
519 // reset of variables for multiple candidates
520 Int_t iC1WithBestDist = 0; // reset
521 Float_t distmin = 100.; // just to put a huge number!
522 Float_t dPhimin = 0.; // Used for histograms only!
523 Float_t dThetamin = 0.; // Used for histograms only!
524 Float_t dZetamin = 0.; // Used for histograms only!
525
526 // in any case, if MC has been required, store statistics of primaries and secondaries
527 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
528 if (fMC) {
529 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
530 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
531 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
532 // do it always as a function of the chip number used to built the prediction
533 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
534 if (!found) {AliWarning(
535 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
536 else {
537 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
538 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
539 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
540 { // this cluster is from a primary particle
541 fClusterPrimary[key]++;
542 primary=kTRUE;
543 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
544 } else { // this cluster is from a secondary particle
545 fClusterSecondary[key]++;
546 secondary=kTRUE;
547 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
548 }
549 }
550 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
551 // (in case the prediction is reliable)
552 if( fChipPredOnLay1[iC2]<1200) {
553 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
554 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
555 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
556 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
557 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
558 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
559 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
560 else fNonRecons[fChipPredOnLay1[iC2]]++;
561 }
562 }
563
564 // Loop on layer 1
565 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
566
567 // The following excludes double associations
568 if (!fAssociationFlag1[iC1]) {
569
570 // find the difference in angles
571 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
572 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
573 // take into account boundary condition
574 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
575
576
577 // find the difference in z (between linear projection from layer 2
578 // and the actual point: Dzeta= z2/r2*r1 -z1)
579 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
580 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
581
582
583 if (fHistOn) {
584 fhClustersDPhiInterpAll->Fill(dPhi);
585 fhClustersDThetaInterpAll->Fill(dTheta);
586 fhClustersDZetaInterpAll->Fill(dZeta);
587 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
588 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
589 }
590 // make "elliptical" cut in Phi and Zeta!
591 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
592 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
593
594 if (d>1) continue;
595
596 //look for the minimum distance: the minimum is in iC1WithBestDist
597 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
598 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
599 dPhimin = dPhi;
600 dThetamin = dTheta;
601 dZetamin = dZeta;
602 iC1WithBestDist = iC1;
603 }
604 }
605 } // end of loop over clusters in layer 1
606
607 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
608
609 if (fHistOn) {
610 fhClustersDPhiInterpAcc->Fill(dPhimin);
611 fhClustersDThetaInterpAcc->Fill(dThetamin);
612 fhClustersDZetaInterpAcc->Fill(dZetamin);
613 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
614 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
615 }
616
617 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
618 // flag the association
619
620 // store the tracklet
621
622 // use the theta from the clusters in the first layer
623 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
624 // use the phi from the clusters in the first layer
625 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
626 // Store the difference between phi1 and phi2
627 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
628
629 // find labels
630 Int_t label1 = 0;
631 Int_t label2 = 0;
632 while (label2 < 3)
633 {
634 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
635 break;
636 label1++;
637 if (label1 == 3)
638 {
639 label1 = 0;
640 label2++;
641 }
642 }
643
644 if (label2 < 3)
645 {
646 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
647 }
648 else
649 {
650 fTracklets[fNTracklets][3] = -2;
651 }
652
653// Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
654 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
655 if(!found){
656 AliWarning(
657 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
658 key=999999;
659 }
660 nfClu1+=(Int_t)found; // this for debugging purpose
661 ntClu1++; // to check efficiency of the method FindChip
662 if(key<1200) {
663 if(fMC) { // this part only for MC
664 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
665 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
666 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
667 if (label2 < 3) { // same label
668 fSuccessTT[key]++;
669 if(primary) fSuccessPP[key]++;
670 }
671 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
672 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
673 }
674
675 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
676 // OK, success
677 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
678 fChipUpdatedInEvent[key]=kTRUE;
679 if(fMC) {
680 if(primary) fSuccessP[key]++;
681 if(secondary) fSuccessS[key]++;
682 }
683 } else {
684 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
685 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
686 if(fMC) {
687 if(primary) fSuccessP[key]++;
688 if(secondary) fSuccessS[key]++;
689 }
690 }
691 }
692
693 fNTracklets++;
694
695 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
696 else if (fChipPredOnLay1[iC2]<1200) {
697 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
698 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
699 if(fMC) {
700 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
701 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
702 }
703 }
704 } // end of loop over clusters in layer 2
705
706 AliDebug(1,Form("%d tracklets found", fNTracklets));
707 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
708 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
709 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
710 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
711}
712//____________________________________________________________________
713Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
714 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
715//
716// Input: a) layer number in the range [0,1]
717// b) vtx[3]: actual vertex
718// c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
719// d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
720// e) phiVtx /
721// Output: Unique key to locate a chip
722// return: kTRUE if succesfull
723
724 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
725 Double_t r=GetRLayer(layer);
726 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
727
728 // set phiVtx in the range [0,2pi]
729 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
730
731 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
732 // of the intersection of the tracklet with the pixel layer.
733 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
734 else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
735 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
736 Double_t vtxy[2]={vtx[0],vtx[1]};
737 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
738 // this method gives you two interceptions
739 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
740 // set phiAbs in the range [0,2pi]
741 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
742 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
743 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
744 // taking the closest one to phiVtx
745 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
746 } else phiAbs=phiVtx;
747 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
748
749 // now you need to locate the chip within the idet detector,
750 // starting from the local coordinates in such a detector
751
752 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
753 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
754 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
755
756 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
757 return kTRUE;
758}
759//______________________________________________________________________________
760Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
761//
762// Return the average radius of a layer from Geometry
763//
764 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
765 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
766
767 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
768 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
769 Double_t r=TMath::Sqrt(x*x + y*y);
770
771 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
772 r += TMath::Sqrt(x*x + y*y);
773 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
774 r += TMath::Sqrt(x*x + y*y);
775 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
776 r += TMath::Sqrt(x*x + y*y);
777 r*=0.25;
778 return r;
779}
780//______________________________________________________________________________
781Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
782 Float_t &xloc, Float_t &zloc) {
783 // this method transform Global Cilindrical coordinates into local (i.e. module)
784 // cartesian coordinates
785 //
786 //Compute Cartesian Global Coordinate
787 Double_t xyzGlob[3],xyzLoc[3];
788 xyzGlob[2]=z;
789 xyzGlob[0]=r*TMath::Cos(phi);
790 xyzGlob[1]=r*TMath::Sin(phi);
791
792 xloc=0.;
793 zloc=0.;
794
795 if(idet<0) return kFALSE;
796
797 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
798 Int_t lad = Int_t(idet/ndet) + 1;
799 Int_t det = idet - (lad-1)*ndet + 1;
800
801 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
802
803 xloc = (Float_t)xyzLoc[0];
804 zloc = (Float_t)xyzLoc[2];
805
806return kTRUE;
807}
808//______________________________________________________________________________
809Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
810 //--------------------------------------------------------------------
811 // This function finds the detector crossed by the track
812 // Input: layer in range [0,1]
813 // phi in ALICE absolute reference system
814 // z " " " " "
815 //--------------------------------------------------------------------
816 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
817 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
818 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
819 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
820
821 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
822 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
823 Double_t phiOffset=TMath::ATan2(y,x);
824 Double_t zOffset=z2;
825
826 Double_t dphi;
827 if (zOffset<0) // old geometry
828 dphi = -(phi-phiOffset);
829 else // new geometry
830 dphi = phi-phiOffset;
831
832 if (dphi < 0) dphi += 2*TMath::Pi();
833 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
834 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
835 if (np>=nladders) np-=nladders;
836 if (np<0) np+=nladders;
837
838 Double_t dz=zOffset-z;
839 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
840 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
841 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
842 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
843
844 return np*ndetectors + nz;
845}
846//____________________________________________________________
847Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
848// this method find the intersection in xy between a tracklet (straight line) and
849// a circonference (r=R), using polar coordinates.
850/*
851Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
852 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
853 - R: radius of the circle
854Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
855
856Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
857r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
858In the same system, the equation of a semi-line is: phi=phiVtx;
859Hence you get one interception only: P=(r,phiVtx)
860Finally you want P in the ABSOLUTE ALICE reference system.
861*/
862Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
863Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
864Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
865Double_t cC=rO*rO-R*R;
866Double_t dDelta=bB*bB-4*cC;
867if(dDelta<0) return kFALSE;
868Double_t r1,r2;
869r1=(-bB-TMath::Sqrt(dDelta))/2;
870r2=(-bB+TMath::Sqrt(dDelta))/2;
871if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
872Double_t r=TMath::Max(r1,r2); // take the positive
873Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
874Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
875pvtx[0]=r*TMath::Cos(phiVtx);
876pvtx[1]=r*TMath::Sin(phiVtx);
877pP[0]=vtx[0]+pvtx[0];
878pP[1]=vtx[1]+pvtx[1];
879phi=TMath::ATan2(pP[1],pP[0]);
880return kTRUE;
881}
882//___________________________________________________________
883Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
884//
885// simple method to reduce all angles (in rad)
886// in range [0,2pi[
887//
888//
889while(angle >=2*TMath::Pi() || angle<0) {
890 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
891 if(angle < 0) angle+=2*TMath::Pi();
892}
893return kTRUE;
894}
895//___________________________________________________________
896Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
897//
898// This method check if a particle is primary; i.e.
899// it comes from the main vertex and it is a "stable" particle, according to
900// AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
901// a stable particle: it has no effect on this analysis).
902// This method can be called only for MC events, where Kinematics is available.
903// if fUseOnlyStableParticle is kTRUE (via SetseOnlyStableParticle) then it
904// returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
905// The latter (see below) try to verify if a primary particle is also "detectable".
906//
907if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
908if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
909if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
910// return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
911 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
912 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
913 TParticle* part = stack->Particle(ipart);
914 TParticle* part0 = stack->Particle(0); // first primary
915 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
916 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
917 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
918
919 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
920 TParticlePDG* pdgPart = part->GetPDG();
921 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
922
923 Double_t distx = part->Vx() - part0->Vx();
924 Double_t disty = part->Vy() - part0->Vy();
925 Double_t distz = part->Vz() - part0->Vz();
926 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
927
928 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
929 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
930 return kFALSE; }// primary if within 500 microns from true Vertex
931
932 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
933 return kTRUE;
934}
935//_____________________________________________________________________________________________
936Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
937//
938// This private method can be applied on MC particles (if stack is available),
939// provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
940//
941// It define "detectable" a primary particle according to the following criteria:
942//
943// - if no decay products can be found in the stack (note that this does not
944// means it is stable, since a particle is stored in stack if it has at least 1 hit in a
945// sensitive detector)
946// - if it has at least one decay daughter produced outside or just on the outer pixel layer
947// - if the last decay particle is an electron (or a muon) which is not produced in-between
948// the two pixel layers (this is likely to be a kink).
949if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
950if(!stack) {AliError("null pointer to MC stack"); return 0;}
951if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
952
953TParticle* part = stack->Particle(ipart);
954//TParticle* part0 = stack->Particle(0); // first primary
955
956 Int_t nret=0;
957 TParticle* dau = 0;
958 Int_t nDau = 0;
959 Int_t pdgDau;
960 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
961 // its real fate ! But you have to take it !
962 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
963 Int_t lastDau = part->GetLastDaughter();
964 nDau = lastDau - firstDau + 1;
965 Double_t distMax=0.;
966 Int_t jmax=0;
967 for(Int_t j=firstDau; j<=lastDau; j++) {
968 dau = stack->Particle(j);
969 Double_t distx = dau->Vx();
970 Double_t disty = dau->Vy();
971 //Double_t distz = dau->Vz();
972 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
973 if(distR<distMax) continue; // considere only the daughter produced at largest radius
974 distMax=distR;
975 jmax=j;
976 }
977 dau = stack->Particle(jmax);
978 pdgDau=dau->GetPdgCode();
979 if (pdgDau == 11 || pdgDau == 13 ) {
980 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
981 else nret =0; // delta-ray emission in material (keep it)
982 }
983 else {// not ele or muon
984 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
985 }
986return nret;
987}
988//_________________________________________________________________
989void AliITSTrackleterSPDEff::InitPredictionMC() {
990//
991// this method allocate memory for the MC related informations
992// all the counters are set to 0
993//
994//
995if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
996fPredictionPrimary = new Int_t[1200];
997fPredictionSecondary = new Int_t[1200];
998fClusterPrimary = new Int_t[1200];
999fClusterSecondary = new Int_t[1200];
1000fSuccessPP = new Int_t[1200];
1001fSuccessTT = new Int_t[1200];
1002fSuccessS = new Int_t[1200];
1003fSuccessP = new Int_t[1200];
1004fFailureS = new Int_t[1200];
1005fFailureP = new Int_t[1200];
1006fRecons = new Int_t[1200];
1007fNonRecons = new Int_t[1200];
1008for(Int_t i=0; i<1200; i++) {
1009 fPredictionPrimary[i]=0;
1010 fPredictionSecondary[i]=0;
1011 fPredictionSecondary[i]=0;
1012 fClusterSecondary[i]=0;
1013 fSuccessPP[i]=0;
1014 fSuccessTT[i]=0;
1015 fSuccessS[i]=0;
1016 fSuccessP[i]=0;
1017 fFailureS[i]=0;
1018 fFailureP[i]=0;
1019 fRecons[i]=0;
1020 fNonRecons[i]=0;
1021}
1022return;
1023}
1024//______________________________________________________________________
1025Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
1026//
1027// This method return the Data menmber fPredictionPrimary [1200].
1028// You can call it only for MC events.
1029// fPredictionPrimary[key] contains the number of tracklet predictions on the
1030// given chip key built using a cluster on the other layer produced (at least)
1031// from a primary particle.
1032// Key refers to the chip crossed by the prediction
1033//
1034//
1035if (!fMC) {CallWarningMC(); return 0;}
1036if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1037return fPredictionPrimary[(Int_t)key];
1038}
1039//______________________________________________________________________
1040Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
1041//
1042// This method return the Data menmber fPredictionSecondary [1200].
1043// You can call it only for MC events.
1044// fPredictionSecondary[key] contains the number of tracklet predictions on the
1045// given chip key built using a cluster on the other layer produced (only)
1046// from a secondary particle
1047// Key refers to the chip crossed by the prediction
1048//
1049//
1050if (!fMC) {CallWarningMC(); return 0;}
1051if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1052return fPredictionSecondary[(Int_t)key];
1053}
1054//______________________________________________________________________
1055Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
1056//
1057// This method return the Data menmber fClusterPrimary [1200].
1058// You can call it only for MC events.
1059// fClusterPrimary[key] contains the number of tracklet predictions
1060// built using a cluster on that layer produced (only)
1061// from a primary particle
1062// Key refers to the chip used to build the prediction
1063//
1064//
1065if (!fMC) {CallWarningMC(); return 0;}
1066if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1067return fClusterPrimary[(Int_t)key];
1068}
1069//______________________________________________________________________
1070Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
1071//
1072// This method return the Data menmber fClusterSecondary [1200].
1073// You can call it only for MC events.
1074// fClusterSecondary[key] contains the number of tracklet predictions
1075// built using a cluster on that layer produced (only)
1076// from a secondary particle
1077// Key refers to the chip used to build the prediction
1078//
1079if (!fMC) {CallWarningMC(); return 0;}
1080if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1081return fClusterSecondary[(Int_t)key];
1082}
1083//______________________________________________________________________
1084Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1085//
1086// This method return the Data menmber fSuccessPP [1200].
1087// You can call it only for MC events.
1088// fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1089// with a cluster on the other layer) built by using the same primary particle
1090// the unique chip key refers to the chip which get updated its efficiency
1091//
1092if (!fMC) {CallWarningMC(); return 0;}
1093if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1094return fSuccessPP[(Int_t)key];
1095}
1096//______________________________________________________________________
1097Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1098//
1099// This method return the Data menmber fSuccessTT [1200].
1100// You can call it only for MC events.
1101// fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1102// with a cluster on the other layer) built by using the same particle (whatever)
1103// the unique chip key refers to the chip which get updated its efficiency
1104//
1105if (!fMC) {CallWarningMC(); return 0;}
1106if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1107return fSuccessTT[(Int_t)key];
1108}
1109//______________________________________________________________________
1110Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1111//
1112// This method return the Data menmber fSuccessS [1200].
1113// You can call it only for MC events.
1114// fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1115// with a cluster on the other layer) built by using a secondary particle
1116// the unique chip key refers to the chip which get updated its efficiency
1117//
1118if (!fMC) {CallWarningMC(); return 0;}
1119if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1120return fSuccessS[(Int_t)key];
1121}
1122//______________________________________________________________________
1123Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1124//
1125// This method return the Data menmber fSuccessP [1200].
1126// You can call it only for MC events.
1127// fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1128// with a cluster on the other layer) built by using a primary particle
1129// the unique chip key refers to the chip which get updated its efficiency
1130//
1131if (!fMC) {CallWarningMC(); return 0;}
1132if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1133return fSuccessP[(Int_t)key];
1134}
1135//______________________________________________________________________
1136Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1137//
1138// This method return the Data menmber fFailureS [1200].
1139// You can call it only for MC events.
1140// fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1141// with a cluster on the other layer) built by using a secondary particle
1142// the unique chip key refers to the chip which get updated its efficiency
1143//
1144if (!fMC) {CallWarningMC(); return 0;}
1145if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1146return fFailureS[(Int_t)key];
1147}
1148//______________________________________________________________________
1149Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1150//
1151// This method return the Data menmber fFailureP [1200].
1152// You can call it only for MC events.
1153// fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1154// with a cluster on the other layer) built by using a primary particle
1155// the unique chip key refers to the chip which get updated its efficiency
1156//
1157if (!fMC) {CallWarningMC(); return 0;}
1158if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1159return fFailureP[(Int_t)key];
1160}
1161//_____________________________________________________________________
1162Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1163//
1164// This method return the Data menmber fRecons [1200].
1165// You can call it only for MC events.
1166// fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1167// has an hit in the detector)
1168// the unique chip key refers to the chip where fall the prediction
1169//
1170if (!fMC) {CallWarningMC(); return 0;}
1171if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1172return fRecons[(Int_t)key];
1173}
1174//_____________________________________________________________________
1175Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1176//
1177// This method return the Data menmber fNonRecons [1200].
1178// You can call it only for MC events.
1179// fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1180// has not any hit in the detector)
1181// the unique chip key refers to the chip where fall the prediction
1182//
1183if (!fMC) {CallWarningMC(); return 0;}
1184if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1185return fNonRecons[(Int_t)key];
1186}
1187//______________________________________________________________________
1188void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1189 // Print out some class data values in Ascii Form to output stream
1190 // Inputs:
1191 // ostream *os Output stream where Ascii data is to be writen
1192 // Outputs:
1193 // none.
1194 // Return:
1195 // none.
1196 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow
1197 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1198 << " " << fUpdateOncePerEventPlaneEff ;
1199 *os << " " << fMC;
1200 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1201 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1202 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1203 << " " << fUseOnlyStableParticle ;
1204 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1205 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1206 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1207 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
1208 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1209 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1210 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1211 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1212 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1213 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1214 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1215 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
1216 return;
1217}
1218//______________________________________________________________________
1219void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1220 // Read in some class data values in Ascii Form to output stream
1221 // Inputs:
1222 // istream *is Input stream where Ascii data is to be read in from
1223 // Outputs:
1224 // none.
1225 // Return:
1226 // none.
1227
1228 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow
1229 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1230 >> fUpdateOncePerEventPlaneEff ;
1231 *is >> fMC;
1232 if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
1233 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1234 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1235 >> fUseOnlyStableParticle;
1236 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1237 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1238 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1239 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1240 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1241 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1242 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1243 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1244 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1245 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1246 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1247 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1248 return;
1249}
1250//______________________________________________________________________
1251ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1252 // Standard output streaming function
1253 // Inputs:
1254 // ostream &os output steam
1255 // AliITSTrackleterSPDEff &s class to be streamed.
1256 // Output:
1257 // none.
1258 // Return:
1259 // ostream &os The stream pointer
1260
1261 s.PrintAscii(&os);
1262 return os;
1263}
1264//______________________________________________________________________
1265istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1266 // Standard inputput streaming function
1267 // Inputs:
1268 // istream &is input steam
1269 // AliITSTrackleterSPDEff &s class to be streamed.
1270 // Output:
1271 // none.
1272 // Return:
1273 // ostream &os The stream pointer
1274
1275 //printf("prova %d \n", (Int_t)s.GetMC());
1276 s.ReadAscii(&is);
1277 return is;
1278}
1279//______________________________________________________________________
1280void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
1281//
1282// This Method write into an asci file (do not know why binary does not work)
1283// the used cuts and the statistics of the MC related quantities
1284// The method SetMC() has to be called before
1285// Input TString filename: name of file for output (it deletes already existing
1286// file)
1287// Output: none
1288//
1289//
1290 if(!fMC) {CallWarningMC(); return;}
1291 ofstream out(filename.Data(),ios::out | ios::binary);
1292 out << *this;
1293 out.close();
1294return;
1295}
1296//____________________________________________________________________
1297void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
1298//
1299// This Method read from an asci file (do not know why binary does not work)
1300// the cuts to be used and the statistics of the MC related quantities
1301// Input TString filename: name of input file for output
1302// The method SetMC() has to be called before
1303// Output: none
1304//
1305//
1306 if(!fMC) {CallWarningMC(); return;}
1307 if( gSystem->AccessPathName( filename.Data() ) ) {
1308 AliError( Form( "file (%s) not found", filename.Data() ) );
1309 return;
1310 }
1311
1312 ifstream in(filename.Data(),ios::in | ios::binary);
1313 in >> *this;
1314 in.close();
1315 return;
1316}
1317//____________________________________________________________________
1318Bool_t AliITSTrackleterSPDEff::SaveHists() {
1319 // This (private) method save the histograms on the output file
1320 // (only if fHistOn is TRUE).
1321 // Also the histograms from the base class are saved through the
1322 // AliITSMultReconstructor::SaveHists() call
1323
1324 if (!GetHistOn()) return kFALSE;
1325
1326 AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1327
1328 fhClustersDPhiInterpAll->Write();
1329 fhClustersDThetaInterpAll->Write();
1330 fhClustersDZetaInterpAll->Write();
1331 fhDPhiVsDThetaInterpAll->Write();
1332 fhDPhiVsDZetaInterpAll->Write();
1333
1334 fhClustersDPhiInterpAcc->Write();
1335 fhClustersDThetaInterpAcc->Write();
1336 fhClustersDZetaInterpAcc->Write();
1337 fhDPhiVsDThetaInterpAcc->Write();
1338 fhDPhiVsDZetaInterpAcc->Write();
1339
1340 fhetaClustersLay2->Write();
1341 fhphiClustersLay2->Write();
1342 return kTRUE;
1343}
1344//__________________________________________________________
1345Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1346 //
1347 // Saves the histograms into a tree and saves the trees into a file
1348 // Also the histograms from the base class are saved
1349 //
1350 if (!GetHistOn()) return kFALSE;
1351 if (filename.Data()=="") {
1352 AliWarning("WriteHistosToFile: null output filename!");
1353 return kFALSE;
1354 }
1355 TFile *hFile=new TFile(filename.Data(),option,
1356 "The File containing the histos for SPD efficiency studies with tracklets");
1357 if(!SaveHists()) return kFALSE;
1358 hFile->Write();
1359 hFile->Close();
1360 return kTRUE;
1361}
1362//____________________________________________________________
1363void AliITSTrackleterSPDEff::BookHistos() {
1364//
1365// This method books addtitional histograms
1366// w.r.t. those of the base class.
1367// In particular, the differences of cluster coordinate between the two SPD
1368// layers are computed in the interpolation phase
1369//
1370 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1371 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1372 fhClustersDPhiInterpAcc->SetDirectory(0);
1373 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1374 fhClustersDThetaInterpAcc->SetDirectory(0);
1375 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1376 fhClustersDZetaInterpAcc->SetDirectory(0);
1377
1378 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1379 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1380 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1381 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1382
1383 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1384 fhClustersDPhiInterpAll->SetDirectory(0);
1385 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1386 fhClustersDThetaInterpAll->SetDirectory(0);
1387 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1388 fhClustersDZetaInterpAll->SetDirectory(0);
1389
1390 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1391 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1392 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1393 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1394
1395 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1396 fhetaClustersLay2->SetDirectory(0);
1397 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1398 fhphiClustersLay2->SetDirectory(0);
1399 return;
1400}
1401//____________________________________________________________
1402void AliITSTrackleterSPDEff::DeleteHistos() {
1403//
1404// Private method to delete Histograms from memory
1405// it is called. e.g., by the destructor.
1406//
1407 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1408 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1409 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1410 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1411 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1412 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1413 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1414 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1415 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1416 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1417 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1418 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1419}
1420//_______________________________________________________________
1421Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1422 Float_t* vtx, AliStack *stack, TTree *ref) {
1423// This (private) method can be used only for MC events, where both AliStack and the TrackReference
1424// are available.
1425// It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1426// Input:
1427// - Int_t layer (either 0 or 1): layer which you want to chech if the tracklete can be
1428// reconstructed at
1429// - Int_t iC : cluster index used to build the tracklet prediction
1430// if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1431// - Float_t* vtx: actual event vertex
1432// - stack: pointer to Stack
1433// - ref: pointer to TTRee of TrackReference
1434Bool_t ret=kFALSE; // returned value
1435Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1436if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1437if(!stack) {AliError("null pointer to MC stack"); return ret;}
1438if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1439if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1440if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1441
1442AliTrackReference *tref=0x0;
1443Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1444Int_t nentries = (Int_t)ref->GetEntries();
1445TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1446TBranch *br = ref->GetBranch("TrackReferences");
1447br->SetAddress(&tcaRef);
1448for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1449 br->GetEntry(itrack);
1450 Int_t nref=tcaRef->GetEntriesFast();
1451 if(nref>0) { //it is enough to look at the first one
1452 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1453 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1454 }
1455}
1456if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1457br->GetEntry(imatch); // redundant, nevertheless ...
1458Int_t nref=tcaRef->GetEntriesFast();
1459for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1460 tref=(AliTrackReference*)tcaRef->At(iref);
1461 if(tref->R()>10) continue; // not SPD ref
1462 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1463 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1464
1465// compute the proper quantities for this tref, as was done for fClustersLay1/2
1466 Float_t x = tref->X() - vtx[0];
1467 Float_t y = tref->Y() - vtx[1];
1468 Float_t z = tref->Z() - vtx[2];
1469
1470 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1471
1472 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1473 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1474 trefLayExtr[2] = z; // Store z
1475
1476 if(layer==1) { // try to see if it is reconstructable at the outer layer
1477// find the difference in angles
1478 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1479 // take into account boundary condition
1480 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1481
1482 // find the difference in z (between linear projection from layer 1
1483 // and the actual point: Dzeta= z1/r1*r2 -z2)
1484 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1485 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1486
1487 // make "elliptical" cut in Phi and Zeta!
1488 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
1489 dZeta*dZeta/fZetaWindow/fZetaWindow);
1490 if (d<1) {ret=kTRUE; break;}
1491 }
1492 if(layer==0) { // try to see if it is reconstructable at the inner layer
1493
1494 // find the difference in angles
1495 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1496 // take into account boundary condition
1497 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1498
1499 // find the difference in z (between linear projection from layer 2
1500 // and the actual point: Dzeta= z2/r2*r1 -z1)
1501 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1502 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1503
1504 // make "elliptical" cut in Phi and Zeta!
1505 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1506 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1507 if (d<1) {ret=kTRUE; break;};
1508 }
1509}
1510delete tcaRef;
1511return ret;
1512}