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