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