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