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