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