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