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