Copy and =operator implemented
[u/mrichter/AliRoot.git] / ITS / AliITSTrackleterSPDEff.cxx
CommitLineData
275a301c 1/**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
275a301c 15//____________________________________________________________________
84161aec 16//
275a301c 17// AliITSTrackleterSPDEff - find SPD chips efficiencies by using tracklets.
84161aec 18//
19// This class has been developed from AliITSMultReconstructor (see
20// it for more details). It is the class for the Trackleter used to estimate
21// SPD plane efficiency.
275a301c 22// The trackleter prediction is built using the vertex and 1 cluster.
23//
84161aec 24//
275a301c 25// Author : Giuseppe Eugenio Bruno, based on the skeleton of Reconstruct method provided by Tiziano Virgili
26// email: giuseppe.bruno@ba.infn.it
84161aec 27//
275a301c 28//____________________________________________________________________
29
84161aec 30/* $Id$ */
31
275a301c 32#include <TFile.h>
a3b31967 33#include <TTree.h>
275a301c 34#include <TParticle.h>
35#include <TSystem.h>
36#include <Riostream.h>
a3b31967 37#include <TClonesArray.h>
275a301c 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"
a3b31967 45#include "AliTrackReference.h"
275a301c 46
47//____________________________________________________________________
48ClassImp(AliITSTrackleterSPDEff)
49
50
51//____________________________________________________________________
52AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
53AliITSMultReconstructor(),
54fAssociationFlag1(0),
55fChipPredOnLay2(0),
56fChipPredOnLay1(0),
57fNTracklets1(0),
58fPhiWindowL1(0),
59fZetaWindowL1(0),
60fOnlyOneTrackletPerC1(0),
a3b31967 61fUpdateOncePerEventPlaneEff(0),
62fChipUpdatedInEvent(0),
275a301c 63fPlaneEffSPD(0),
64fMC(0),
65fUseOnlyPrimaryForPred(0),
66fUseOnlySecondaryForPred(0),
67fUseOnlySameParticle(0),
68fUseOnlyDifferentParticle(0),
69fUseOnlyStableParticle(0),
70fPredictionPrimary(0),
71fPredictionSecondary(0),
72fClusterPrimary(0),
73fClusterSecondary(0),
a3b31967 74fSuccessPP(0),
75fSuccessTT(0),
76fSuccessS(0),
77fSuccessP(0),
78fFailureS(0),
79fFailureP(0),
80fRecons(0),
81fNonRecons(0),
275a301c 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{
84161aec 95 // default constructor
275a301c 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];
a3b31967 104 fChipUpdatedInEvent = new Bool_t[1200];
275a301c 105
106 for(Int_t i=0; i<300000; i++) {
107 fAssociationFlag1[i] = kFALSE;
108 }
a3b31967 109 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
275a301c 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),
a3b31967 124fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
125fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
275a301c 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),
a3b31967 137fSuccessPP(mr.fSuccessPP),
138fSuccessTT(mr.fSuccessTT),
139fSuccessS(mr.fSuccessS),
140fSuccessP(mr.fSuccessP),
141fFailureS(mr.fFailureS),
142fFailureP(mr.fFailureP),
143fRecons(mr.fRecons),
144fNonRecons(mr.fNonRecons),
275a301c 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
a3b31967 180 delete [] fChipUpdatedInEvent;
181
275a301c 182 delete [] fPredictionPrimary;
183 delete [] fPredictionSecondary;
184 delete [] fClusterPrimary;
185 delete [] fClusterSecondary;
a3b31967 186 delete [] fSuccessPP;
187 delete [] fSuccessTT;
188 delete [] fSuccessS;
189 delete [] fSuccessP;
190 delete [] fFailureS;
191 delete [] fFailureP;
192 delete [] fRecons;
193 delete [] fNonRecons;
275a301c 194
195 // delete PlaneEff
196 delete fPlaneEffSPD;
197}
198//____________________________________________________________________
199void
a3b31967 200AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t*, AliStack *pStack, TTree *tRef) {
275a301c 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;}
a3b31967 218 if(fMC && !tRef) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
275a301c 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
a3b31967 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;
275a301c 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 }
275a301c 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
a3b31967 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
275a301c 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
a3b31967 319 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
275a301c 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]++;
a3b31967 334 primary=kTRUE;
275a301c 335 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
336 } else { // this cluster is from a secondary particle
337 fClusterSecondary[key]++;
a3b31967 338 secondary=kTRUE;
275a301c 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]]++;
a3b31967 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]]++;
275a301c 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];
a3b31967 466 if (label2 < 3) {
467 fSuccessTT[key]++;
468 if(primary) fSuccessPP[key]++;
469 }
275a301c 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 !
a3b31967 475 // OK, success
275a301c 476 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
a3b31967 477 fChipUpdatedInEvent[key]=kTRUE;
478 if(fMC) {
479 if(primary) fSuccessP[key]++;
480 if(secondary) fSuccessS[key]++;
481 }
275a301c 482 }
483 else {
484 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
a3b31967 485 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
486 if(fMC) {
487 if(primary) fSuccessP[key]++;
488 if(secondary) fSuccessS[key]++;
489 }
275a301c 490 }
491 }
492
493 fNTracklets++;
494
495 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
a3b31967 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 }
275a301c 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
a3b31967 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
275a301c 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
a3b31967 527 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
275a301c 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]++;
a3b31967 542 primary=kTRUE;
275a301c 543 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
544 } else { // this cluster is from a secondary particle
545 fClusterSecondary[key]++;
a3b31967 546 secondary=kTRUE;
275a301c 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]]++;
a3b31967 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]]++;
275a301c 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
a3b31967 607 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
275a301c 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];
a3b31967 667 if (label2 < 3) { // same label
668 fSuccessTT[key]++;
669 if(primary) fSuccessPP[key]++;
670 }
275a301c 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 !
a3b31967 676 // OK, success
275a301c 677 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
a3b31967 678 fChipUpdatedInEvent[key]=kTRUE;
679 if(fMC) {
680 if(primary) fSuccessP[key]++;
681 if(secondary) fSuccessS[key]++;
682 }
275a301c 683 } else {
684 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
a3b31967 685 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
686 if(fMC) {
687 if(primary) fSuccessP[key]++;
688 if(secondary) fSuccessS[key]++;
689 }
275a301c 690 }
691 }
692
693 fNTracklets++;
694
695 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
a3b31967 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 }
275a301c 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) {
a3b31967 761//
762// Return the average radius of a layer from Geometry
763//
275a301c 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 //--------------------------------------------------------------------
a3b31967 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 " " " " "
275a301c 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)
a3b31967 860Finally you want P in the ABSOLUTE ALICE reference system.
275a301c 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) {
a3b31967 884//
885// simple method to reduce all angles (in rad)
886// in range [0,2pi[
887//
888//
275a301c 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) {
a3b31967 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//
275a301c 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
a3b31967 932 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
275a301c 933 return kTRUE;
934}
935//_____________________________________________________________________________________________
936Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
a3b31967 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).
275a301c 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;
a3b31967 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
275a301c 963 Int_t lastDau = part->GetLastDaughter();
964 nDau = lastDau - firstDau + 1;
a3b31967 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)
275a301c 985 }
a3b31967 986return nret;
275a301c 987}
988//_________________________________________________________________
989void AliITSTrackleterSPDEff::InitPredictionMC() {
a3b31967 990//
991// this method allocate memory for the MC related informations
992// all the counters are set to 0
993//
994//
275a301c 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];
a3b31967 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];
275a301c 1008for(Int_t i=0; i<1200; i++) {
1009 fPredictionPrimary[i]=0;
1010 fPredictionSecondary[i]=0;
1011 fPredictionSecondary[i]=0;
1012 fClusterSecondary[i]=0;
a3b31967 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;
275a301c 1021}
1022return;
1023}
1024//______________________________________________________________________
1025Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
84161aec 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//
275a301c 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 {
84161aec 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//
275a301c 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 {
84161aec 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//
275a301c 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 {
84161aec 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//
275a301c 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//______________________________________________________________________
a3b31967 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//______________________________________________________________________
275a301c 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.
a3b31967 1196 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow
1197 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1198 << " " << fUpdateOncePerEventPlaneEff ;
275a301c 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) ;
a3b31967 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) ;
275a301c 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
a3b31967 1228 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow
1229 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1230 >> fUpdateOncePerEventPlaneEff ;
275a301c 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] ;
a3b31967 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] ;
275a301c 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 {
84161aec 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//
275a301c 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) {
84161aec 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//
275a301c 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() {
84161aec 1319 // This (private) method save the histograms on the output file
275a301c 1320 // (only if fHistOn is TRUE).
84161aec 1321 // Also the histograms from the base class are saved through the
1322 // AliITSMultReconstructor::SaveHists() call
275a301c 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
84161aec 1348 // Also the histograms from the base class are saved
275a301c 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() {
84161aec 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//
275a301c 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() {
84161aec 1403//
1404// Private method to delete Histograms from memory
1405// it is called. e.g., by the destructor.
1406//
275a301c 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//_______________________________________________________________
a3b31967 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}