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