]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSTrackleterSPDEff.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix warning
[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
58e8dc31 39#include "AliTracker.h"
275a301c 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"
58e8dc31 46#include "AliRunLoader.h"
47#include "AliITSReconstructor.h"
48#include "AliITSRecPoint.h"
275a301c 49//____________________________________________________________________
50ClassImp(AliITSTrackleterSPDEff)
51
52
53//____________________________________________________________________
54AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
58e8dc31 55AliTracker(),
56//
57fClustersLay1(0),
58fClustersLay2(0),
59fTracklets(0),
60fAssociationFlag(0),
61fNClustersLay1(0),
62fNClustersLay2(0),
63fNTracklets(0),
64fOnlyOneTrackletPerC2(0),
65fPhiWindow(0),
66fZetaWindow(0),
67fPhiOverlapCut(0),
68fZetaOverlapCut(0),
69fHistOn(0),
70fhClustersDPhiAcc(0),
71fhClustersDThetaAcc(0),
72fhClustersDZetaAcc(0),
73fhClustersDPhiAll(0),
74fhClustersDThetaAll(0),
75fhClustersDZetaAll(0),
76fhDPhiVsDThetaAll(0),
77fhDPhiVsDThetaAcc(0),
78fhDPhiVsDZetaAll(0),
79fhDPhiVsDZetaAcc(0),
80fhetaTracklets(0),
81fhphiTracklets(0),
82fhetaClustersLay1(0),
83fhphiClustersLay1(0),
84//
275a301c 85fAssociationFlag1(0),
86fChipPredOnLay2(0),
87fChipPredOnLay1(0),
88fNTracklets1(0),
89fPhiWindowL1(0),
90fZetaWindowL1(0),
91fOnlyOneTrackletPerC1(0),
a3b31967 92fUpdateOncePerEventPlaneEff(0),
93fChipUpdatedInEvent(0),
275a301c 94fPlaneEffSPD(0),
0fce916f 95fReflectClusterAroundZAxisForLayer0(kFALSE),
96fReflectClusterAroundZAxisForLayer1(kFALSE),
275a301c 97fMC(0),
58e8dc31 98fUseOnlyPrimaryForPred(1),
275a301c 99fUseOnlySecondaryForPred(0),
100fUseOnlySameParticle(0),
101fUseOnlyDifferentParticle(0),
102fUseOnlyStableParticle(0),
103fPredictionPrimary(0),
104fPredictionSecondary(0),
105fClusterPrimary(0),
106fClusterSecondary(0),
a3b31967 107fSuccessPP(0),
108fSuccessTT(0),
109fSuccessS(0),
110fSuccessP(0),
111fFailureS(0),
112fFailureP(0),
113fRecons(0),
114fNonRecons(0),
275a301c 115fhClustersDPhiInterpAcc(0),
116fhClustersDThetaInterpAcc(0),
117fhClustersDZetaInterpAcc(0),
118fhClustersDPhiInterpAll(0),
119fhClustersDThetaInterpAll(0),
120fhClustersDZetaInterpAll(0),
121fhDPhiVsDThetaInterpAll(0),
122fhDPhiVsDThetaInterpAcc(0),
123fhDPhiVsDZetaInterpAll(0),
124fhDPhiVsDZetaInterpAcc(0),
125fhetaClustersLay2(0),
126fhphiClustersLay2(0)
127{
84161aec 128 // default constructor
58e8dc31 129// from AliITSMultReconstructor
130 SetPhiWindow();
131 SetZetaWindow();
132 SetOnlyOneTrackletPerC2();
133 fClustersLay1 = new Float_t*[300000];
134 fClustersLay2 = new Float_t*[300000];
135 fTracklets = new Float_t*[300000];
136 fAssociationFlag = new Bool_t[300000];
137//
275a301c 138 SetPhiWindowL1();
139 SetZetaWindowL1();
140 SetOnlyOneTrackletPerC1();
141
142 fAssociationFlag1 = new Bool_t[300000];
143 fChipPredOnLay2 = new UInt_t[300000];
144 fChipPredOnLay1 = new UInt_t[300000];
a3b31967 145 fChipUpdatedInEvent = new Bool_t[1200];
275a301c 146
147 for(Int_t i=0; i<300000; i++) {
58e8dc31 148 // from AliITSMultReconstructor
149 fClustersLay1[i] = new Float_t[6];
150 fClustersLay2[i] = new Float_t[6];
151 fTracklets[i] = new Float_t[5];
152 fAssociationFlag[i] = kFALSE;
153 //
275a301c 154 fAssociationFlag1[i] = kFALSE;
155 }
a3b31967 156 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
275a301c 157
158 if (GetHistOn()) BookHistos();
159
160 fPlaneEffSPD = new AliITSPlaneEffSPD();
161}
162//______________________________________________________________________
58e8dc31 163AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
164AliTracker(mr),
165// from AliITSMultReconstructor
166fClustersLay1(mr.fClustersLay1),
167fClustersLay2(mr.fClustersLay2),
168fTracklets(mr.fTracklets),
169fAssociationFlag(mr.fAssociationFlag),
170fNClustersLay1(mr.fNClustersLay1),
171fNClustersLay2(mr.fNClustersLay2),
172fNTracklets(mr.fNTracklets),
173fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
174fPhiWindow(mr.fPhiWindow),
175fZetaWindow(mr.fZetaWindow),
176fPhiOverlapCut(mr.fPhiOverlapCut),
177fZetaOverlapCut(mr.fZetaOverlapCut),
178fHistOn(mr.fHistOn),
179fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
180fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
181fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
182fhClustersDPhiAll(mr.fhClustersDPhiAll),
183fhClustersDThetaAll(mr.fhClustersDThetaAll),
184fhClustersDZetaAll(mr.fhClustersDZetaAll),
185fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
186fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
187fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
188fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
189fhetaTracklets(mr.fhetaTracklets),
190fhphiTracklets(mr.fhphiTracklets),
191fhetaClustersLay1(mr.fhetaClustersLay1),
192fhphiClustersLay1(mr.fhphiClustersLay1),
193//
275a301c 194fAssociationFlag1(mr.fAssociationFlag1),
195fChipPredOnLay2(mr.fChipPredOnLay2),
196fChipPredOnLay1(mr.fChipPredOnLay1),
197fNTracklets1(mr.fNTracklets1),
198fPhiWindowL1(mr.fPhiWindowL1),
199fZetaWindowL1(mr.fZetaWindowL1),
200fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
a3b31967 201fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
202fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
275a301c 203fPlaneEffSPD(mr.fPlaneEffSPD),
0fce916f 204fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
205fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
275a301c 206fMC(mr.fMC),
207fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
208fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
209fUseOnlySameParticle(mr.fUseOnlySameParticle),
210fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
211fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
212fPredictionPrimary(mr.fPredictionPrimary),
213fPredictionSecondary(mr.fPredictionSecondary),
214fClusterPrimary(mr.fClusterPrimary),
215fClusterSecondary(mr.fClusterSecondary),
a3b31967 216fSuccessPP(mr.fSuccessPP),
217fSuccessTT(mr.fSuccessTT),
218fSuccessS(mr.fSuccessS),
219fSuccessP(mr.fSuccessP),
220fFailureS(mr.fFailureS),
221fFailureP(mr.fFailureP),
222fRecons(mr.fRecons),
223fNonRecons(mr.fNonRecons),
275a301c 224fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
225fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
226fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
227fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
228fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
229fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
230fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
231fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
232fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
233fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
234fhetaClustersLay2(mr.fhetaClustersLay2),
235fhphiClustersLay2(mr.fhphiClustersLay2)
236{
237 // Copy constructor
238}
239
240//______________________________________________________________________
241AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
242 // Assignment operator
243 this->~AliITSTrackleterSPDEff();
244 new(this) AliITSTrackleterSPDEff(mr);
245 return *this;
246}
247//______________________________________________________________________
248AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
249 // Destructor
58e8dc31 250// from AliITSMultReconstructor
251 // delete arrays
252 for(Int_t i=0; i<300000; i++) {
253 delete [] fClustersLay1[i];
254 delete [] fClustersLay2[i];
255 delete [] fTracklets[i];
256 }
257 delete [] fClustersLay1;
258 delete [] fClustersLay2;
259 delete [] fTracklets;
260 delete [] fAssociationFlag;
261//
275a301c 262 // delete histograms
263 DeleteHistos();
264
265 delete [] fAssociationFlag1;
266
267 delete [] fChipPredOnLay2;
268 delete [] fChipPredOnLay1;
269
a3b31967 270 delete [] fChipUpdatedInEvent;
271
275a301c 272 delete [] fPredictionPrimary;
273 delete [] fPredictionSecondary;
274 delete [] fClusterPrimary;
275 delete [] fClusterSecondary;
a3b31967 276 delete [] fSuccessPP;
277 delete [] fSuccessTT;
278 delete [] fSuccessS;
279 delete [] fSuccessP;
280 delete [] fFailureS;
281 delete [] fFailureP;
282 delete [] fRecons;
283 delete [] fNonRecons;
275a301c 284
285 // delete PlaneEff
286 delete fPlaneEffSPD;
287}
288//____________________________________________________________________
289void
58e8dc31 290//AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t*, AliStack *pStack, TTree *tRef) {
291AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef) {
275a301c 292 //
58e8dc31 293 // - you have to take care of the following, before of using Reconstruct
294 // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
295 // and convert the cluster coordinates to theta, phi (seen from the
296 // interaction vertex).
297 // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
298 // - Find the extrapolation/interpolation point.
275a301c 299 // - Find the chip corresponding to that
300 // - Check if there is a cluster near that point
301 //
275a301c 302 // reset counters
275a301c 303 fNTracklets = 0;
58e8dc31 304 // retrieve the vertex position
305 Float_t vtx[3];
306 vtx[0]=(Float_t)GetX();
307 vtx[1]=(Float_t)GetY();
308 vtx[2]=(Float_t)GetZ();
0fce916f 309 // to study residual background (i.e. contribution from TT' to measured efficiency)
310 if(fReflectClusterAroundZAxisForLayer0) ReflectClusterAroundZAxisForLayer(0);
311 if(fReflectClusterAroundZAxisForLayer1) ReflectClusterAroundZAxisForLayer(1);
312 //
275a301c 313 if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
a3b31967 314 if(fMC && !tRef) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
275a301c 315 Bool_t found;
316 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
317 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
318 Int_t nfClu1=0; Int_t ntClu1=0;
319 Int_t nfClu2=0; Int_t ntClu2=0;
320
a3b31967 321 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
322 // for this new event)
323 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
275a301c 324
325 // find the tracklets
326 AliDebug(1,"Looking for tracklets... ");
327 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
328
329 //###########################################################
330 // Loop on layer 1 : finding theta, phi and z
331 UInt_t key;
332 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
333 Float_t x = fClustersLay1[iC1][0] - vtx[0];
334 Float_t y = fClustersLay1[iC1][1] - vtx[1];
335 Float_t z = fClustersLay1[iC1][2] - vtx[2];
336
337 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
338
339 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
340 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
341 fClustersLay1[iC1][2] = z; // Store z
342
343 // find the Radius and the chip corresponding to the extrapolation point
344
345 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
346 if (!found) {
347 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
348 key=999999; // also some other actions should be taken if not Found
349 }
350 nfTraPred2+=(Int_t)found; // this for debugging purpose
351 ntTraPred2++; // to check efficiency of the method FindChip
352 fChipPredOnLay2[iC1] = key;
353 fAssociationFlag1[iC1] = kFALSE;
354
355 if (fHistOn) {
356 Float_t eta=fClustersLay1[iC1][0];
357 eta= TMath::Tan(eta/2.);
358 eta=-TMath::Log(eta);
359 fhetaClustersLay1->Fill(eta);
360 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
361 }
362 }
275a301c 363 // Loop on layer 2 : finding theta, phi and r
364 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
365 Float_t x = fClustersLay2[iC2][0] - vtx[0];
366 Float_t y = fClustersLay2[iC2][1] - vtx[1];
367 Float_t z = fClustersLay2[iC2][2] - vtx[2];
368
369 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
370
371 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
372 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
373 fClustersLay2[iC2][2] = z; // Store z
374
375 // find the Radius and the chip corresponding to the extrapolation point
376
377 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
378 if (!found) {
379 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
380 key=999999;
381 }
382 nfTraPred1+=(Int_t)found; // this for debugging purpose
383 ntTraPred1++; // to check efficiency of the method FindChip
384 fChipPredOnLay1[iC2] = key;
385 fAssociationFlag[iC2] = kFALSE;
386
387 if (fHistOn) {
388 Float_t eta=fClustersLay2[iC2][0];
389 eta= TMath::Tan(eta/2.);
390 eta=-TMath::Log(eta);
391 fhetaClustersLay2->Fill(eta);
392 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
393 }
394 }
395
396 //###########################################################
397
398 // First part : Extrapolation to Layer 2
399
400 // Loop on layer 1
401 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
402
a3b31967 403 // here the control to check whether the efficiency of the chip traversed by this tracklet
404 // prediction has already been updated in this event using another tracklet prediction
405 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
406
275a301c 407 // reset of variables for multiple candidates
408 Int_t iC2WithBestDist = 0; // reset
409 Float_t distmin = 100.; // just to put a huge number!
410 Float_t dPhimin = 0.; // Used for histograms only!
411 Float_t dThetamin = 0.; // Used for histograms only!
412 Float_t dZetamin = 0.; // Used for histograms only!
413
414 // in any case, if MC has been required, store statistics of primaries and secondaries
a3b31967 415 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
275a301c 416 if (fMC) {
417 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
418 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
419 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
420 // do it always as a function of the chip number used to built the prediction
421 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
422 if (!found) {AliWarning(
423 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
424 else {
425 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
426 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
427 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
428 { // this cluster is from a primary particle
429 fClusterPrimary[key]++;
a3b31967 430 primary=kTRUE;
275a301c 431 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
432 } else { // this cluster is from a secondary particle
433 fClusterSecondary[key]++;
a3b31967 434 secondary=kTRUE;
275a301c 435 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
436 }
437 }
438 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
439 // (in case the prediction is reliable)
440 if( fChipPredOnLay2[iC1]<1200) {
441 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
442 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
443 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
444 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
a3b31967 445 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
446 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
447 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
448 else fNonRecons[fChipPredOnLay2[iC1]]++;
275a301c 449 }
450 }
451
452 // Loop on layer 2
453 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
454
455 // The following excludes double associations
456 if (!fAssociationFlag[iC2]) {
457
458 // find the difference in angles
459 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
460 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
461 // take into account boundary condition
462 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
463
464 // find the difference in z (between linear projection from layer 1
465 // and the actual point: Dzeta= z1/r1*r2 -z2)
466 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
467 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
468
469 if (fHistOn) {
470 fhClustersDPhiAll->Fill(dPhi);
471 fhClustersDThetaAll->Fill(dTheta);
472 fhClustersDZetaAll->Fill(dZeta);
473 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
474 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
475 }
476
477 // make "elliptical" cut in Phi and Zeta!
478 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
479 dZeta*dZeta/fZetaWindow/fZetaWindow);
480
481 if (d>1) continue;
482
483 //look for the minimum distance: the minimum is in iC2WithBestDist
484 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
485 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
486 dPhimin = dPhi;
487 dThetamin = dTheta;
488 dZetamin = dZeta;
489 iC2WithBestDist = iC2;
490 }
491 }
492 } // end of loop over clusters in layer 2
493
494 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
495
496 if (fHistOn) {
497 fhClustersDPhiAcc->Fill(dPhimin);
498 fhClustersDThetaAcc->Fill(dThetamin);
499 fhClustersDZetaAcc->Fill(dZetamin);
500 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
501 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
502 }
503
504 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
505 // flag the association
506
507 // store the tracklet
508
509 // use the theta from the clusters in the first layer
510 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
511 // use the phi from the clusters in the first layer
512 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
513 // Store the difference between phi1 and phi2
514 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
515
516 // find labels
517 Int_t label1 = 0;
518 Int_t label2 = 0;
519 while (label2 < 3)
520 {
521 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
522 break;
523 label1++;
524 if (label1 == 3)
525 {
526 label1 = 0;
527 label2++;
528 }
529 }
530
531 if (label2 < 3)
532 {
533 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
534 }
535 else
536 {
537 fTracklets[fNTracklets][3] = -2;
538 }
539
540 if (fHistOn) {
541 Float_t eta=fTracklets[fNTracklets][0];
542 eta= TMath::Tan(eta/2.);
543 eta=-TMath::Log(eta);
544 fhetaTracklets->Fill(eta);
545 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
546 }
547
548// Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
549 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
550 if(!found){
551 AliWarning(
552 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
553 key=999999;
554 }
555 nfClu2+=(Int_t)found; // this for debugging purpose
556 ntClu2++; // to check efficiency of the method FindChip
557 if(key<1200) { // the Chip has been found
558 if(fMC) { // this part only for MC
559 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
560 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
561 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
a3b31967 562 if (label2 < 3) {
563 fSuccessTT[key]++;
564 if(primary) fSuccessPP[key]++;
565 }
275a301c 566 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
567 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
568 }
569
570 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
a3b31967 571 // OK, success
275a301c 572 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
a3b31967 573 fChipUpdatedInEvent[key]=kTRUE;
574 if(fMC) {
575 if(primary) fSuccessP[key]++;
576 if(secondary) fSuccessS[key]++;
577 }
275a301c 578 }
579 else {
580 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
a3b31967 581 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
582 if(fMC) {
583 if(primary) fSuccessP[key]++;
584 if(secondary) fSuccessS[key]++;
585 }
275a301c 586 }
587 }
588
589 fNTracklets++;
590
591 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
a3b31967 592 else if (fChipPredOnLay2[iC1]<1200) {
593 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
594 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
595 if(fMC) {
596 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
597 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
598 }
599 }
275a301c 600 } // end of loop over clusters in layer 1
601
602 fNTracklets1=fNTracklets;
603
604//###################################################################
605
606 // Second part : Interpolation to Layer 1
607
608 // Loop on layer 2
609 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
610
a3b31967 611 // here the control to check whether the efficiency of the chip traversed by this tracklet
612 // prediction has already been updated in this event using another tracklet prediction
613 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
614
275a301c 615 // reset of variables for multiple candidates
616 Int_t iC1WithBestDist = 0; // reset
617 Float_t distmin = 100.; // just to put a huge number!
618 Float_t dPhimin = 0.; // Used for histograms only!
619 Float_t dThetamin = 0.; // Used for histograms only!
620 Float_t dZetamin = 0.; // Used for histograms only!
621
622 // in any case, if MC has been required, store statistics of primaries and secondaries
a3b31967 623 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
275a301c 624 if (fMC) {
625 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
626 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
627 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
628 // do it always as a function of the chip number used to built the prediction
629 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
630 if (!found) {AliWarning(
631 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
632 else {
633 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
634 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
635 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
636 { // this cluster is from a primary particle
637 fClusterPrimary[key]++;
a3b31967 638 primary=kTRUE;
275a301c 639 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
640 } else { // this cluster is from a secondary particle
641 fClusterSecondary[key]++;
a3b31967 642 secondary=kTRUE;
275a301c 643 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
644 }
645 }
646 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
647 // (in case the prediction is reliable)
648 if( fChipPredOnLay1[iC2]<1200) {
649 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
650 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
651 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
652 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
a3b31967 653 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
654 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
655 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
656 else fNonRecons[fChipPredOnLay1[iC2]]++;
275a301c 657 }
658 }
659
660 // Loop on layer 1
661 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
662
663 // The following excludes double associations
664 if (!fAssociationFlag1[iC1]) {
665
666 // find the difference in angles
667 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
668 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
669 // take into account boundary condition
670 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
671
672
673 // find the difference in z (between linear projection from layer 2
674 // and the actual point: Dzeta= z2/r2*r1 -z1)
675 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
676 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
677
678
679 if (fHistOn) {
680 fhClustersDPhiInterpAll->Fill(dPhi);
681 fhClustersDThetaInterpAll->Fill(dTheta);
682 fhClustersDZetaInterpAll->Fill(dZeta);
683 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
684 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
685 }
686 // make "elliptical" cut in Phi and Zeta!
687 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
688 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
689
690 if (d>1) continue;
691
692 //look for the minimum distance: the minimum is in iC1WithBestDist
693 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
694 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
695 dPhimin = dPhi;
696 dThetamin = dTheta;
697 dZetamin = dZeta;
698 iC1WithBestDist = iC1;
699 }
700 }
701 } // end of loop over clusters in layer 1
702
a3b31967 703 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
275a301c 704
705 if (fHistOn) {
706 fhClustersDPhiInterpAcc->Fill(dPhimin);
707 fhClustersDThetaInterpAcc->Fill(dThetamin);
708 fhClustersDZetaInterpAcc->Fill(dZetamin);
709 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
710 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
711 }
712
713 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
714 // flag the association
715
716 // store the tracklet
717
718 // use the theta from the clusters in the first layer
719 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
720 // use the phi from the clusters in the first layer
721 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
722 // Store the difference between phi1 and phi2
723 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
724
725 // find labels
726 Int_t label1 = 0;
727 Int_t label2 = 0;
728 while (label2 < 3)
729 {
730 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
731 break;
732 label1++;
733 if (label1 == 3)
734 {
735 label1 = 0;
736 label2++;
737 }
738 }
739
740 if (label2 < 3)
741 {
742 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
743 }
744 else
745 {
746 fTracklets[fNTracklets][3] = -2;
747 }
748
749// Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
750 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
751 if(!found){
752 AliWarning(
753 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
754 key=999999;
755 }
756 nfClu1+=(Int_t)found; // this for debugging purpose
757 ntClu1++; // to check efficiency of the method FindChip
758 if(key<1200) {
759 if(fMC) { // this part only for MC
760 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
761 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
762 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
a3b31967 763 if (label2 < 3) { // same label
764 fSuccessTT[key]++;
765 if(primary) fSuccessPP[key]++;
766 }
275a301c 767 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
768 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
769 }
770
771 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
a3b31967 772 // OK, success
275a301c 773 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
a3b31967 774 fChipUpdatedInEvent[key]=kTRUE;
775 if(fMC) {
776 if(primary) fSuccessP[key]++;
777 if(secondary) fSuccessS[key]++;
778 }
275a301c 779 } else {
780 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
a3b31967 781 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
782 if(fMC) {
783 if(primary) fSuccessP[key]++;
784 if(secondary) fSuccessS[key]++;
785 }
275a301c 786 }
787 }
788
789 fNTracklets++;
790
791 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
a3b31967 792 else if (fChipPredOnLay1[iC2]<1200) {
793 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
794 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
795 if(fMC) {
796 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
797 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
798 }
799 }
275a301c 800 } // end of loop over clusters in layer 2
801
802 AliDebug(1,Form("%d tracklets found", fNTracklets));
803 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
804 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
805 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
806 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
807}
808//____________________________________________________________________
809Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
810 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
811//
812// Input: a) layer number in the range [0,1]
813// b) vtx[3]: actual vertex
814// c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
815// d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
816// e) phiVtx /
817// Output: Unique key to locate a chip
818// return: kTRUE if succesfull
819
820 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
821 Double_t r=GetRLayer(layer);
822 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
823
824 // set phiVtx in the range [0,2pi]
825 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
826
827 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
828 // of the intersection of the tracklet with the pixel layer.
829 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
830 else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
831 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
832 Double_t vtxy[2]={vtx[0],vtx[1]};
833 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
834 // this method gives you two interceptions
835 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
836 // set phiAbs in the range [0,2pi]
837 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
838 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
839 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
840 // taking the closest one to phiVtx
841 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
842 } else phiAbs=phiVtx;
843 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
844
845 // now you need to locate the chip within the idet detector,
846 // starting from the local coordinates in such a detector
847
848 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
849 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
850 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
851
852 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
853 return kTRUE;
854}
855//______________________________________________________________________________
856Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
a3b31967 857//
858// Return the average radius of a layer from Geometry
859//
275a301c 860 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
861 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
862
863 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
864 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
865 Double_t r=TMath::Sqrt(x*x + y*y);
866
867 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
868 r += TMath::Sqrt(x*x + y*y);
869 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
870 r += TMath::Sqrt(x*x + y*y);
871 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
872 r += TMath::Sqrt(x*x + y*y);
873 r*=0.25;
874 return r;
875}
876//______________________________________________________________________________
877Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
878 Float_t &xloc, Float_t &zloc) {
879 // this method transform Global Cilindrical coordinates into local (i.e. module)
880 // cartesian coordinates
881 //
882 //Compute Cartesian Global Coordinate
883 Double_t xyzGlob[3],xyzLoc[3];
884 xyzGlob[2]=z;
885 xyzGlob[0]=r*TMath::Cos(phi);
886 xyzGlob[1]=r*TMath::Sin(phi);
887
888 xloc=0.;
889 zloc=0.;
890
891 if(idet<0) return kFALSE;
892
893 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
894 Int_t lad = Int_t(idet/ndet) + 1;
895 Int_t det = idet - (lad-1)*ndet + 1;
896
897 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
898
899 xloc = (Float_t)xyzLoc[0];
900 zloc = (Float_t)xyzLoc[2];
901
902return kTRUE;
903}
904//______________________________________________________________________________
905Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
906 //--------------------------------------------------------------------
a3b31967 907 // This function finds the detector crossed by the track
908 // Input: layer in range [0,1]
909 // phi in ALICE absolute reference system
910 // z " " " " "
275a301c 911 //--------------------------------------------------------------------
912 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
913 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
914 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
915 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
916
917 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
918 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
919 Double_t phiOffset=TMath::ATan2(y,x);
920 Double_t zOffset=z2;
921
922 Double_t dphi;
923 if (zOffset<0) // old geometry
924 dphi = -(phi-phiOffset);
925 else // new geometry
926 dphi = phi-phiOffset;
927
928 if (dphi < 0) dphi += 2*TMath::Pi();
929 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
930 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
931 if (np>=nladders) np-=nladders;
932 if (np<0) np+=nladders;
933
934 Double_t dz=zOffset-z;
935 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
936 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
937 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
938 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
939
940 return np*ndetectors + nz;
941}
942//____________________________________________________________
943Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
944// this method find the intersection in xy between a tracklet (straight line) and
945// a circonference (r=R), using polar coordinates.
946/*
947Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
948 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
949 - R: radius of the circle
950Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
951
952Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
953r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
954In the same system, the equation of a semi-line is: phi=phiVtx;
955Hence you get one interception only: P=(r,phiVtx)
a3b31967 956Finally you want P in the ABSOLUTE ALICE reference system.
275a301c 957*/
958Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
959Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
960Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
961Double_t cC=rO*rO-R*R;
962Double_t dDelta=bB*bB-4*cC;
963if(dDelta<0) return kFALSE;
964Double_t r1,r2;
965r1=(-bB-TMath::Sqrt(dDelta))/2;
966r2=(-bB+TMath::Sqrt(dDelta))/2;
967if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
968Double_t r=TMath::Max(r1,r2); // take the positive
969Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
970Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
971pvtx[0]=r*TMath::Cos(phiVtx);
972pvtx[1]=r*TMath::Sin(phiVtx);
973pP[0]=vtx[0]+pvtx[0];
974pP[1]=vtx[1]+pvtx[1];
975phi=TMath::ATan2(pP[1],pP[0]);
976return kTRUE;
977}
978//___________________________________________________________
979Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
a3b31967 980//
981// simple method to reduce all angles (in rad)
982// in range [0,2pi[
983//
984//
275a301c 985while(angle >=2*TMath::Pi() || angle<0) {
986 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
987 if(angle < 0) angle+=2*TMath::Pi();
988}
989return kTRUE;
990}
991//___________________________________________________________
992Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
a3b31967 993//
994// This method check if a particle is primary; i.e.
995// it comes from the main vertex and it is a "stable" particle, according to
996// AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
997// a stable particle: it has no effect on this analysis).
998// This method can be called only for MC events, where Kinematics is available.
58e8dc31 999// if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
a3b31967 1000// returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
1001// The latter (see below) try to verify if a primary particle is also "detectable".
1002//
275a301c 1003if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
1004if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
1005if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
1006// return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
1007 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
1008 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
1009 TParticle* part = stack->Particle(ipart);
1010 TParticle* part0 = stack->Particle(0); // first primary
1011 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
1012 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
1013 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
1014
1015 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
1016 TParticlePDG* pdgPart = part->GetPDG();
1017 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
1018
1019 Double_t distx = part->Vx() - part0->Vx();
1020 Double_t disty = part->Vy() - part0->Vy();
1021 Double_t distz = part->Vz() - part0->Vz();
1022 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
1023
1024 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
1025 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
1026 return kFALSE; }// primary if within 500 microns from true Vertex
1027
a3b31967 1028 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
275a301c 1029 return kTRUE;
1030}
1031//_____________________________________________________________________________________________
1032Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
a3b31967 1033//
1034// This private method can be applied on MC particles (if stack is available),
1035// provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
1036//
1037// It define "detectable" a primary particle according to the following criteria:
1038//
1039// - if no decay products can be found in the stack (note that this does not
1040// means it is stable, since a particle is stored in stack if it has at least 1 hit in a
1041// sensitive detector)
1042// - if it has at least one decay daughter produced outside or just on the outer pixel layer
1043// - if the last decay particle is an electron (or a muon) which is not produced in-between
1044// the two pixel layers (this is likely to be a kink).
275a301c 1045if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
1046if(!stack) {AliError("null pointer to MC stack"); return 0;}
1047if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
1048
1049TParticle* part = stack->Particle(ipart);
1050//TParticle* part0 = stack->Particle(0); // first primary
1051
1052 Int_t nret=0;
1053 TParticle* dau = 0;
1054 Int_t nDau = 0;
a3b31967 1055 Int_t pdgDau;
1056 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
1057 // its real fate ! But you have to take it !
1058 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
275a301c 1059 Int_t lastDau = part->GetLastDaughter();
1060 nDau = lastDau - firstDau + 1;
a3b31967 1061 Double_t distMax=0.;
1062 Int_t jmax=0;
1063 for(Int_t j=firstDau; j<=lastDau; j++) {
1064 dau = stack->Particle(j);
1065 Double_t distx = dau->Vx();
1066 Double_t disty = dau->Vy();
1067 //Double_t distz = dau->Vz();
1068 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
1069 if(distR<distMax) continue; // considere only the daughter produced at largest radius
1070 distMax=distR;
1071 jmax=j;
1072 }
1073 dau = stack->Particle(jmax);
1074 pdgDau=dau->GetPdgCode();
1075 if (pdgDau == 11 || pdgDau == 13 ) {
1076 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
1077 else nret =0; // delta-ray emission in material (keep it)
1078 }
1079 else {// not ele or muon
1080 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
275a301c 1081 }
a3b31967 1082return nret;
275a301c 1083}
1084//_________________________________________________________________
1085void AliITSTrackleterSPDEff::InitPredictionMC() {
a3b31967 1086//
1087// this method allocate memory for the MC related informations
1088// all the counters are set to 0
1089//
1090//
275a301c 1091if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1092fPredictionPrimary = new Int_t[1200];
1093fPredictionSecondary = new Int_t[1200];
1094fClusterPrimary = new Int_t[1200];
1095fClusterSecondary = new Int_t[1200];
a3b31967 1096fSuccessPP = new Int_t[1200];
1097fSuccessTT = new Int_t[1200];
1098fSuccessS = new Int_t[1200];
1099fSuccessP = new Int_t[1200];
1100fFailureS = new Int_t[1200];
1101fFailureP = new Int_t[1200];
1102fRecons = new Int_t[1200];
1103fNonRecons = new Int_t[1200];
275a301c 1104for(Int_t i=0; i<1200; i++) {
1105 fPredictionPrimary[i]=0;
1106 fPredictionSecondary[i]=0;
1107 fPredictionSecondary[i]=0;
1108 fClusterSecondary[i]=0;
a3b31967 1109 fSuccessPP[i]=0;
1110 fSuccessTT[i]=0;
1111 fSuccessS[i]=0;
1112 fSuccessP[i]=0;
1113 fFailureS[i]=0;
1114 fFailureP[i]=0;
1115 fRecons[i]=0;
1116 fNonRecons[i]=0;
275a301c 1117}
1118return;
1119}
c6a05d92 1120//_________________________________________________________________
1121void AliITSTrackleterSPDEff::DeletePredictionMC() {
1122//
1123// this method deallocate memory for the MC related informations
1124// all the counters are set to 0
1125//
1126//
1127if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1128if(fPredictionPrimary) {
1129 delete fPredictionPrimary; fPredictionPrimary=0;
1130}
1131if(fPredictionSecondary) {
1132 delete fPredictionSecondary; fPredictionSecondary=0;
1133}
1134if(fClusterPrimary) {
1135 delete fClusterPrimary; fClusterPrimary=0;
1136}
1137if(fClusterSecondary) {
1138 delete fClusterSecondary; fClusterSecondary=0;
1139}
1140if(fSuccessPP) {
1141 delete fSuccessPP; fSuccessPP=0;
1142}
1143if(fSuccessTT) {
1144 delete fSuccessTT; fSuccessTT=0;
1145}
1146if(fSuccessS) {
1147 delete fSuccessS; fSuccessS=0;
1148}
1149if(fSuccessP) {
1150 delete fSuccessP; fSuccessP=0;
1151}
1152if(fFailureS) {
1153 delete fFailureS; fFailureS=0;
1154}
1155if(fFailureP) {
1156 delete fFailureP; fFailureP=0;
1157}
1158if(fRecons) {
1159 delete fRecons; fRecons=0;
1160}
1161if(fNonRecons) {
1162 delete fNonRecons; fNonRecons=0;
1163}
1164return;
1165}
275a301c 1166//______________________________________________________________________
1167Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
84161aec 1168//
1169// This method return the Data menmber fPredictionPrimary [1200].
1170// You can call it only for MC events.
1171// fPredictionPrimary[key] contains the number of tracklet predictions on the
1172// given chip key built using a cluster on the other layer produced (at least)
1173// from a primary particle.
1174// Key refers to the chip crossed by the prediction
1175//
1176//
275a301c 1177if (!fMC) {CallWarningMC(); return 0;}
1178if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1179return fPredictionPrimary[(Int_t)key];
1180}
1181//______________________________________________________________________
1182Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
84161aec 1183//
1184// This method return the Data menmber fPredictionSecondary [1200].
1185// You can call it only for MC events.
1186// fPredictionSecondary[key] contains the number of tracklet predictions on the
1187// given chip key built using a cluster on the other layer produced (only)
1188// from a secondary particle
1189// Key refers to the chip crossed by the prediction
1190//
1191//
275a301c 1192if (!fMC) {CallWarningMC(); return 0;}
1193if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1194return fPredictionSecondary[(Int_t)key];
1195}
1196//______________________________________________________________________
1197Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
84161aec 1198//
1199// This method return the Data menmber fClusterPrimary [1200].
1200// You can call it only for MC events.
1201// fClusterPrimary[key] contains the number of tracklet predictions
1202// built using a cluster on that layer produced (only)
1203// from a primary particle
1204// Key refers to the chip used to build the prediction
1205//
1206//
275a301c 1207if (!fMC) {CallWarningMC(); return 0;}
1208if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1209return fClusterPrimary[(Int_t)key];
1210}
1211//______________________________________________________________________
1212Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
84161aec 1213//
1214// This method return the Data menmber fClusterSecondary [1200].
1215// You can call it only for MC events.
1216// fClusterSecondary[key] contains the number of tracklet predictions
1217// built using a cluster on that layer produced (only)
1218// from a secondary particle
1219// Key refers to the chip used to build the prediction
1220//
275a301c 1221if (!fMC) {CallWarningMC(); return 0;}
1222if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1223return fClusterSecondary[(Int_t)key];
1224}
1225//______________________________________________________________________
a3b31967 1226Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1227//
1228// This method return the Data menmber fSuccessPP [1200].
1229// You can call it only for MC events.
1230// fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1231// with a cluster on the other layer) built by using the same primary particle
1232// the unique chip key refers to the chip which get updated its efficiency
1233//
1234if (!fMC) {CallWarningMC(); return 0;}
1235if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1236return fSuccessPP[(Int_t)key];
1237}
1238//______________________________________________________________________
1239Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1240//
1241// This method return the Data menmber fSuccessTT [1200].
1242// You can call it only for MC events.
1243// fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1244// with a cluster on the other layer) built by using the same particle (whatever)
1245// the unique chip key refers to the chip which get updated its efficiency
1246//
1247if (!fMC) {CallWarningMC(); return 0;}
1248if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1249return fSuccessTT[(Int_t)key];
1250}
1251//______________________________________________________________________
1252Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1253//
1254// This method return the Data menmber fSuccessS [1200].
1255// You can call it only for MC events.
1256// fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1257// with a cluster on the other layer) built by using a secondary particle
1258// the unique chip key refers to the chip which get updated its efficiency
1259//
1260if (!fMC) {CallWarningMC(); return 0;}
1261if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1262return fSuccessS[(Int_t)key];
1263}
1264//______________________________________________________________________
1265Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1266//
1267// This method return the Data menmber fSuccessP [1200].
1268// You can call it only for MC events.
1269// fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1270// with a cluster on the other layer) built by using a primary particle
1271// the unique chip key refers to the chip which get updated its efficiency
1272//
1273if (!fMC) {CallWarningMC(); return 0;}
1274if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1275return fSuccessP[(Int_t)key];
1276}
1277//______________________________________________________________________
1278Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1279//
1280// This method return the Data menmber fFailureS [1200].
1281// You can call it only for MC events.
1282// fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1283// with a cluster on the other layer) built by using a secondary particle
1284// the unique chip key refers to the chip which get updated its efficiency
1285//
1286if (!fMC) {CallWarningMC(); return 0;}
1287if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1288return fFailureS[(Int_t)key];
1289}
1290//______________________________________________________________________
1291Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1292//
1293// This method return the Data menmber fFailureP [1200].
1294// You can call it only for MC events.
1295// fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1296// with a cluster on the other layer) built by using a primary particle
1297// the unique chip key refers to the chip which get updated its efficiency
1298//
1299if (!fMC) {CallWarningMC(); return 0;}
1300if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1301return fFailureP[(Int_t)key];
1302}
1303//_____________________________________________________________________
1304Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1305//
1306// This method return the Data menmber fRecons [1200].
1307// You can call it only for MC events.
1308// fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1309// has an hit in the detector)
1310// the unique chip key refers to the chip where fall the prediction
1311//
1312if (!fMC) {CallWarningMC(); return 0;}
1313if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1314return fRecons[(Int_t)key];
1315}
1316//_____________________________________________________________________
1317Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1318//
1319// This method return the Data menmber fNonRecons [1200].
1320// You can call it only for MC events.
1321// fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1322// has not any hit in the detector)
1323// the unique chip key refers to the chip where fall the prediction
1324//
1325if (!fMC) {CallWarningMC(); return 0;}
1326if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1327return fNonRecons[(Int_t)key];
1328}
1329//______________________________________________________________________
275a301c 1330void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1331 // Print out some class data values in Ascii Form to output stream
1332 // Inputs:
1333 // ostream *os Output stream where Ascii data is to be writen
1334 // Outputs:
1335 // none.
1336 // Return:
1337 // none.
a3b31967 1338 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow
1339 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
0fce916f 1340 << " " << fUpdateOncePerEventPlaneEff << " " << fReflectClusterAroundZAxisForLayer0
1341 << " " << fReflectClusterAroundZAxisForLayer1;
275a301c 1342 *os << " " << fMC;
1343 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1344 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1345 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1346 << " " << fUseOnlyStableParticle ;
1347 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1348 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1349 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1350 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
a3b31967 1351 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1352 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1353 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1354 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1355 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1356 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1357 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1358 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
275a301c 1359 return;
1360}
1361//______________________________________________________________________
1362void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1363 // Read in some class data values in Ascii Form to output stream
1364 // Inputs:
1365 // istream *is Input stream where Ascii data is to be read in from
1366 // Outputs:
1367 // none.
1368 // Return:
1369 // none.
1370
c6a05d92 1371 Bool_t tmp= fMC;
a3b31967 1372 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow
1373 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
0fce916f 1374 >> fUpdateOncePerEventPlaneEff >> fReflectClusterAroundZAxisForLayer0
1375 >> fReflectClusterAroundZAxisForLayer1;
c6a05d92 1376 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
275a301c 1377 *is >> fMC;
c6a05d92 1378 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1379 else {
1380 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1381 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1382 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1383 >> fUseOnlyStableParticle;
1384 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1385 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1386 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1387 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1388 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1389 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1390 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1391 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1392 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1393 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1394 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1395 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1396 }
275a301c 1397 return;
1398}
1399//______________________________________________________________________
1400ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1401 // Standard output streaming function
1402 // Inputs:
1403 // ostream &os output steam
1404 // AliITSTrackleterSPDEff &s class to be streamed.
1405 // Output:
1406 // none.
1407 // Return:
1408 // ostream &os The stream pointer
1409
1410 s.PrintAscii(&os);
1411 return os;
1412}
1413//______________________________________________________________________
1414istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1415 // Standard inputput streaming function
1416 // Inputs:
1417 // istream &is input steam
1418 // AliITSTrackleterSPDEff &s class to be streamed.
1419 // Output:
1420 // none.
1421 // Return:
1422 // ostream &os The stream pointer
1423
1424 //printf("prova %d \n", (Int_t)s.GetMC());
1425 s.ReadAscii(&is);
1426 return is;
1427}
1428//______________________________________________________________________
1429void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
84161aec 1430//
58e8dc31 1431// This Method write into an either asci or root file
84161aec 1432// the used cuts and the statistics of the MC related quantities
1433// The method SetMC() has to be called before
1434// Input TString filename: name of file for output (it deletes already existing
1435// file)
1436// Output: none
1437//
1438//
c6a05d92 1439 //if(!fMC) {CallWarningMC(); return;}
1440 if (!filename.Contains(".root")) {
1441 ofstream out(filename.Data(),ios::out | ios::binary);
1442 out << *this;
1443 out.close();
1444 return;
1445 }
1446 else {
1447 TFile* mcfile = TFile::Open(filename, "RECREATE");
1448 TH1F* cuts = new TH1F("cuts", "list of cuts", 10, 0, 10); // TH1I containing cuts
1449 cuts->SetBinContent(1,fPhiWindowL1);
1450 cuts->SetBinContent(2,fZetaWindowL1);
1451 cuts->SetBinContent(3,fPhiWindow);
1452 cuts->SetBinContent(4,fZetaWindow);
1453 cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
1454 cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
1455 cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
1456 cuts->SetBinContent(8,fReflectClusterAroundZAxisForLayer0);
1457 cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer1);
1458 cuts->SetBinContent(10,fMC);
1459 cuts->Write();
1460 delete cuts;
1461 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
1462 else {
1463 TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
1464 mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
1465 mc0->SetBinContent(2,fUseOnlySecondaryForPred);
1466 mc0->SetBinContent(3,fUseOnlySameParticle);
1467 mc0->SetBinContent(4,fUseOnlyDifferentParticle);
1468 mc0->SetBinContent(5,fUseOnlyStableParticle);
1469 mc0->Write();
1470 delete mc0;
1471 TH1I *mc1;
1472 mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
1473 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
1474 mc1->Write();
1475 mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
1476 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
1477 mc1->Write();
1478 mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
1479 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
1480 mc1->Write();
1481 mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
1482 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
1483 mc1->Write();
1484 mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
1485 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
1486 mc1->Write();
1487 mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
1488 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
1489 mc1->Write();
1490 mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
1491 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
1492 mc1->Write();
1493 mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
1494 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
1495 mc1->Write();
1496 mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
1497 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
1498 mc1->Write();
1499 mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
1500 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
1501 mc1->Write();
1502 mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
1503 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
1504 mc1->Write();
1505 mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
1506 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
1507 mc1->Write();
1508 delete mc1;
1509 }
1510 mcfile->Close();
1511 }
275a301c 1512return;
1513}
1514//____________________________________________________________________
1515void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
84161aec 1516//
1517// This Method read from an asci file (do not know why binary does not work)
1518// the cuts to be used and the statistics of the MC related quantities
1519// Input TString filename: name of input file for output
1520// The method SetMC() has to be called before
1521// Output: none
1522//
1523//
c6a05d92 1524 //if(!fMC) {CallWarningMC(); return;}
275a301c 1525 if( gSystem->AccessPathName( filename.Data() ) ) {
1526 AliError( Form( "file (%s) not found", filename.Data() ) );
1527 return;
1528 }
1529
c6a05d92 1530 if (!filename.Contains(".root")) {
1531 ifstream in(filename.Data(),ios::in | ios::binary);
1532 in >> *this;
1533 in.close();
1534 return;
1535 }
1536 else {
1537 Bool_t tmp= fMC;
1538 TFile *mcfile = TFile::Open(filename);
1539 TH1F *cuts = (TH1F*)mcfile->Get("cuts");
1540 fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
1541 fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
1542 fPhiWindow=(Float_t)cuts->GetBinContent(3);
1543 fZetaWindow=(Float_t)cuts->GetBinContent(4);
1544 fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
1545 fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
1546 fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
1547 fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(8);
1548 fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(9);
1549 fMC=(Bool_t)cuts->GetBinContent(10);
1550 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1551 else { // only if file with MC predictions
1552 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1553 TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
1554 fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
1555 fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
1556 fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
1557 fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
1558 fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
1559 TH1I *mc1;
1560 mc1 =(TH1I*)mcfile->Get("mc1");
1561 for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1562 mc1 =(TH1I*)mcfile->Get("mc2");
1563 for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1564 mc1 =(TH1I*)mcfile->Get("mc3");
1565 for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1566 mc1 =(TH1I*)mcfile->Get("mc4");
1567 for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1568 mc1 =(TH1I*)mcfile->Get("mc5");
1569 for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1570 mc1 =(TH1I*)mcfile->Get("mc6");
1571 for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
1572 mc1 =(TH1I*)mcfile->Get("mc7");
1573 for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1574 mc1 =(TH1I*)mcfile->Get("mc8");
1575 for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1576 mc1 =(TH1I*)mcfile->Get("mc9");
1577 for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1578 mc1 =(TH1I*)mcfile->Get("mc10");
1579 for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1580 mc1 =(TH1I*)mcfile->Get("mc11");
1581 for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1582 mc1 =(TH1I*)mcfile->Get("mc12");
1583 for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1584 }
1585 mcfile->Close();
1586 }
275a301c 1587 return;
1588}
1589//____________________________________________________________________
1590Bool_t AliITSTrackleterSPDEff::SaveHists() {
84161aec 1591 // This (private) method save the histograms on the output file
275a301c 1592 // (only if fHistOn is TRUE).
84161aec 1593 // Also the histograms from the base class are saved through the
1594 // AliITSMultReconstructor::SaveHists() call
275a301c 1595
1596 if (!GetHistOn()) return kFALSE;
1597
58e8dc31 1598// AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1599 fhClustersDPhiAll->Write();
1600 fhClustersDThetaAll->Write();
1601 fhClustersDZetaAll->Write();
1602 fhDPhiVsDThetaAll->Write();
1603 fhDPhiVsDZetaAll->Write();
1604
1605 fhClustersDPhiAcc->Write();
1606 fhClustersDThetaAcc->Write();
1607 fhClustersDZetaAcc->Write();
1608 fhDPhiVsDThetaAcc->Write();
1609 fhDPhiVsDZetaAcc->Write();
1610
1611 fhetaTracklets->Write();
1612 fhphiTracklets->Write();
1613 fhetaClustersLay1->Write();
1614 fhphiClustersLay1->Write();
275a301c 1615
1616 fhClustersDPhiInterpAll->Write();
1617 fhClustersDThetaInterpAll->Write();
1618 fhClustersDZetaInterpAll->Write();
1619 fhDPhiVsDThetaInterpAll->Write();
1620 fhDPhiVsDZetaInterpAll->Write();
1621
1622 fhClustersDPhiInterpAcc->Write();
1623 fhClustersDThetaInterpAcc->Write();
1624 fhClustersDZetaInterpAcc->Write();
1625 fhDPhiVsDThetaInterpAcc->Write();
1626 fhDPhiVsDZetaInterpAcc->Write();
1627
1628 fhetaClustersLay2->Write();
1629 fhphiClustersLay2->Write();
1630 return kTRUE;
1631}
1632//__________________________________________________________
1633Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1634 //
1635 // Saves the histograms into a tree and saves the trees into a file
84161aec 1636 // Also the histograms from the base class are saved
275a301c 1637 //
1638 if (!GetHistOn()) return kFALSE;
58e8dc31 1639 if (filename.Data()=="") {
275a301c 1640 AliWarning("WriteHistosToFile: null output filename!");
1641 return kFALSE;
1642 }
1643 TFile *hFile=new TFile(filename.Data(),option,
1644 "The File containing the histos for SPD efficiency studies with tracklets");
1645 if(!SaveHists()) return kFALSE;
1646 hFile->Write();
1647 hFile->Close();
1648 return kTRUE;
1649}
1650//____________________________________________________________
1651void AliITSTrackleterSPDEff::BookHistos() {
84161aec 1652//
1653// This method books addtitional histograms
1654// w.r.t. those of the base class.
1655// In particular, the differences of cluster coordinate between the two SPD
1656// layers are computed in the interpolation phase
1657//
275a301c 1658 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
58e8dc31 1659//
1660 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
1661 fhClustersDPhiAcc->SetDirectory(0);
1662 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
1663 fhClustersDThetaAcc->SetDirectory(0);
1664 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
1665 fhClustersDZetaAcc->SetDirectory(0);
1666
1667 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
1668 fhDPhiVsDZetaAcc->SetDirectory(0);
1669 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
1670 fhDPhiVsDThetaAcc->SetDirectory(0);
1671
1672 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
1673 fhClustersDPhiAll->SetDirectory(0);
1674 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
1675 fhClustersDThetaAll->SetDirectory(0);
1676 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
1677 fhClustersDZetaAll->SetDirectory(0);
1678
1679 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
1680 fhDPhiVsDZetaAll->SetDirectory(0);
1681 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
1682 fhDPhiVsDThetaAll->SetDirectory(0);
1683
1684 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
1685 fhetaTracklets->SetDirectory(0);
1686 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
1687 fhphiTracklets->SetDirectory(0);
1688 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
1689 fhetaClustersLay1->SetDirectory(0);
1690 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
1691 fhphiClustersLay1->SetDirectory(0);
1692//
275a301c 1693 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1694 fhClustersDPhiInterpAcc->SetDirectory(0);
1695 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1696 fhClustersDThetaInterpAcc->SetDirectory(0);
1697 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1698 fhClustersDZetaInterpAcc->SetDirectory(0);
1699
1700 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1701 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1702 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1703 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1704
1705 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1706 fhClustersDPhiInterpAll->SetDirectory(0);
1707 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1708 fhClustersDThetaInterpAll->SetDirectory(0);
1709 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1710 fhClustersDZetaInterpAll->SetDirectory(0);
1711
1712 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1713 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1714 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1715 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1716
1717 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1718 fhetaClustersLay2->SetDirectory(0);
1719 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1720 fhphiClustersLay2->SetDirectory(0);
1721 return;
1722}
1723//____________________________________________________________
1724void AliITSTrackleterSPDEff::DeleteHistos() {
84161aec 1725//
1726// Private method to delete Histograms from memory
1727// it is called. e.g., by the destructor.
58e8dc31 1728//
1729// form AliITSMultReconstructor
1730 if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
1731 if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
1732 if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
1733 if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
1734 if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
1735 if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
1736 if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
1737 if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
1738 if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
1739 if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
1740 if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
1741 if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
1742 if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
1743 if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
84161aec 1744//
275a301c 1745 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1746 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1747 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1748 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1749 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1750 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1751 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1752 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1753 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1754 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1755 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1756 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1757}
1758//_______________________________________________________________
a3b31967 1759Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1760 Float_t* vtx, AliStack *stack, TTree *ref) {
1761// This (private) method can be used only for MC events, where both AliStack and the TrackReference
1762// are available.
1763// It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1764// Input:
1765// - Int_t layer (either 0 or 1): layer which you want to chech if the tracklete can be
1766// reconstructed at
1767// - Int_t iC : cluster index used to build the tracklet prediction
1768// if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1769// - Float_t* vtx: actual event vertex
1770// - stack: pointer to Stack
1771// - ref: pointer to TTRee of TrackReference
1772Bool_t ret=kFALSE; // returned value
1773Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1774if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1775if(!stack) {AliError("null pointer to MC stack"); return ret;}
1776if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1777if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1778if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1779
1780AliTrackReference *tref=0x0;
1781Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1782Int_t nentries = (Int_t)ref->GetEntries();
1783TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1784TBranch *br = ref->GetBranch("TrackReferences");
1785br->SetAddress(&tcaRef);
1786for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1787 br->GetEntry(itrack);
1788 Int_t nref=tcaRef->GetEntriesFast();
1789 if(nref>0) { //it is enough to look at the first one
1790 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1791 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1792 }
1793}
1794if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1795br->GetEntry(imatch); // redundant, nevertheless ...
1796Int_t nref=tcaRef->GetEntriesFast();
1797for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1798 tref=(AliTrackReference*)tcaRef->At(iref);
1799 if(tref->R()>10) continue; // not SPD ref
1800 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1801 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1802
1803// compute the proper quantities for this tref, as was done for fClustersLay1/2
1804 Float_t x = tref->X() - vtx[0];
1805 Float_t y = tref->Y() - vtx[1];
1806 Float_t z = tref->Z() - vtx[2];
1807
1808 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1809
1810 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1811 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1812 trefLayExtr[2] = z; // Store z
1813
1814 if(layer==1) { // try to see if it is reconstructable at the outer layer
1815// find the difference in angles
1816 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1817 // take into account boundary condition
1818 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1819
1820 // find the difference in z (between linear projection from layer 1
1821 // and the actual point: Dzeta= z1/r1*r2 -z2)
1822 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1823 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1824
1825 // make "elliptical" cut in Phi and Zeta!
1826 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
1827 dZeta*dZeta/fZetaWindow/fZetaWindow);
1828 if (d<1) {ret=kTRUE; break;}
1829 }
1830 if(layer==0) { // try to see if it is reconstructable at the inner layer
1831
1832 // find the difference in angles
1833 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1834 // take into account boundary condition
1835 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1836
1837 // find the difference in z (between linear projection from layer 2
1838 // and the actual point: Dzeta= z2/r2*r1 -z1)
1839 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1840 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1841
1842 // make "elliptical" cut in Phi and Zeta!
1843 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1844 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1845 if (d<1) {ret=kTRUE; break;};
1846 }
1847}
1848delete tcaRef;
1849return ret;
1850}
0fce916f 1851//_________________________________________________________________________
1852void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
1853//
1854// this method apply a rotation by 180 degree around the Z (beam) axis to all
1855// the RecPoints in a given layer to be used to build tracklets.
1856// **************** VERY IMPORTANT:: ***************
1857// It must be called just after LoadClusterArrays, since afterwards the datamember
1858// fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
1859// instead of Cartesian
1860//
1861if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1862AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1863if(ilayer==0) {
1864 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1865 fClustersLay1[iC1][0]*=-1;
1866 fClustersLay1[iC1][1]*=-1;
1867 }
1868}
1869if(ilayer==1) {
1870 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1871 fClustersLay2[iC2][0]*=-1;
1872 fClustersLay2[iC2][1]*=-1;
1873 }
1874}
1875return;
1876}
58e8dc31 1877//____________________________________________________________________________
1878Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *){
1879// This method is used to find the tracklets.
1880// It is called from AliReconstruction
1881// The vertex is supposed to be associated to the Tracker (i.e. to this) already
1882// The cluster is supposed to be associated to the Tracker already
1883// In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
1884//
1885 Int_t rc=1;
1886 AliRunLoader* runLoader = AliRunLoader::Instance();
1887 if (!runLoader) {
1888 Error("Clusters2Tracks", "no run loader found");
1889 return rc;
1890 }
1891 AliStack *pStack=0x0; TTree *tRefTree=0x0;
1892 if(GetMC()) {
1893 runLoader->LoadKinematics("read");
1894 runLoader->LoadTrackRefs("read");
1895 pStack= runLoader->Stack();
1896 tRefTree= runLoader->TreeTR();
1897 }
1898 Reconstruct(pStack,tRefTree);
1899 return 0;
1900}
1901//____________________________________________________________________________
1902Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
1903//
1904// It is called from AliReconstruction
1905//
1906//
1907//
1908//
1909 Int_t rc=0;
1910 if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
1911 if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
1912 return rc;
1913}
1914//____________________________________________________________________
1915void
1916AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
1917 // This method
1918 // - gets the clusters from the cluster tree
1919 // - convert them into global coordinates
1920 // - store them in the internal arrays
1921 // - count the number of cluster-fired chips
1922
1923 //AliDebug(1,"Loading clusters and cluster-fired chips ...");
1924
1925 fNClustersLay1 = 0;
1926 fNClustersLay2 = 0;
1927
1928 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
1929 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
1930
1931 itsClusterBranch->SetAddress(&itsClusters);
1932
1933 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
1934 Float_t cluGlo[3]={0.,0.,0.};
1935
1936 // loop over the its subdetectors
1937 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
1938
1939 if (!itsClusterTree->GetEvent(iIts))
1940 continue;
1941
1942 Int_t nClusters = itsClusters->GetEntriesFast();
1943
1944 // number of clusters in each chip of the current module
1945 Int_t layer = 0;
1946
1947 // loop over clusters
1948 while(nClusters--) {
1949 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
1950
1951 layer = cluster->GetLayer();
1952 if (layer>1) continue;
1953
1954 cluster->GetGlobalXYZ(cluGlo);
1955 Float_t x = cluGlo[0];
1956 Float_t y = cluGlo[1];
1957 Float_t z = cluGlo[2];
1958
1959 if (layer==0) {
1960 fClustersLay1[fNClustersLay1][0] = x;
1961 fClustersLay1[fNClustersLay1][1] = y;
1962 fClustersLay1[fNClustersLay1][2] = z;
1963
1964 for (Int_t i=0; i<3; i++)
1965 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
1966 fNClustersLay1++;
1967 }
1968 if (layer==1) {
1969 fClustersLay2[fNClustersLay2][0] = x;
1970 fClustersLay2[fNClustersLay2][1] = y;
1971 fClustersLay2[fNClustersLay2][2] = z;
1972
1973 for (Int_t i=0; i<3; i++)
1974 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
1975 fNClustersLay2++;
1976 }
1977
1978 }// end of cluster loop
1979
1980 } // end of its "subdetector" loop
1981 if (itsClusters) {
1982 itsClusters->Delete();
1983 delete itsClusters;
1984 itsClusters = 0;
1985 }
1986 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
1987}
1988