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