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