1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
15 //____________________________________________________________________
17 // AliITSTrackleterSPDEff - find SPD chips efficiencies by using tracklets.
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.
22 // The trackleter prediction is built using the vertex and 1 cluster.
25 // Author : Giuseppe Eugenio Bruno, based on the skeleton of Reconstruct method provided by Tiziano Virgili
26 // email: giuseppe.bruno@ba.infn.it
28 //____________________________________________________________________
34 #include <TParticle.h>
36 #include <Riostream.h>
37 #include <TClonesArray.h>
39 #include "AliITSMultReconstructor.h"
40 #include "AliITSTrackleterSPDEff.h"
41 #include "AliITSgeomTGeo.h"
43 #include "AliITSPlaneEffSPD.h"
45 #include "AliTrackReference.h"
47 //____________________________________________________________________
48 ClassImp(AliITSTrackleterSPDEff)
51 //____________________________________________________________________
52 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
53 AliITSMultReconstructor(),
60 fOnlyOneTrackletPerC1(0),
61 fUpdateOncePerEventPlaneEff(0),
62 fChipUpdatedInEvent(0),
64 fReflectClusterAroundZAxisForLayer0(kFALSE),
65 fReflectClusterAroundZAxisForLayer1(kFALSE),
67 fUseOnlyPrimaryForPred(0),
68 fUseOnlySecondaryForPred(0),
69 fUseOnlySameParticle(0),
70 fUseOnlyDifferentParticle(0),
71 fUseOnlyStableParticle(0),
72 fPredictionPrimary(0),
73 fPredictionSecondary(0),
84 fhClustersDPhiInterpAcc(0),
85 fhClustersDThetaInterpAcc(0),
86 fhClustersDZetaInterpAcc(0),
87 fhClustersDPhiInterpAll(0),
88 fhClustersDThetaInterpAll(0),
89 fhClustersDZetaInterpAll(0),
90 fhDPhiVsDThetaInterpAll(0),
91 fhDPhiVsDThetaInterpAcc(0),
92 fhDPhiVsDZetaInterpAll(0),
93 fhDPhiVsDZetaInterpAcc(0),
97 // default constructor
101 SetOnlyOneTrackletPerC1();
103 fAssociationFlag1 = new Bool_t[300000];
104 fChipPredOnLay2 = new UInt_t[300000];
105 fChipPredOnLay1 = new UInt_t[300000];
106 fChipUpdatedInEvent = new Bool_t[1200];
108 for(Int_t i=0; i<300000; i++) {
109 fAssociationFlag1[i] = kFALSE;
111 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
113 if (GetHistOn()) BookHistos();
115 fPlaneEffSPD = new AliITSPlaneEffSPD();
117 //______________________________________________________________________
118 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) : AliITSMultReconstructor(mr),
119 fAssociationFlag1(mr.fAssociationFlag1),
120 fChipPredOnLay2(mr.fChipPredOnLay2),
121 fChipPredOnLay1(mr.fChipPredOnLay1),
122 fNTracklets1(mr.fNTracklets1),
123 fPhiWindowL1(mr.fPhiWindowL1),
124 fZetaWindowL1(mr.fZetaWindowL1),
125 fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
126 fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
127 fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
128 fPlaneEffSPD(mr.fPlaneEffSPD),
129 fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
130 fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
132 fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
133 fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
134 fUseOnlySameParticle(mr.fUseOnlySameParticle),
135 fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
136 fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
137 fPredictionPrimary(mr.fPredictionPrimary),
138 fPredictionSecondary(mr.fPredictionSecondary),
139 fClusterPrimary(mr.fClusterPrimary),
140 fClusterSecondary(mr.fClusterSecondary),
141 fSuccessPP(mr.fSuccessPP),
142 fSuccessTT(mr.fSuccessTT),
143 fSuccessS(mr.fSuccessS),
144 fSuccessP(mr.fSuccessP),
145 fFailureS(mr.fFailureS),
146 fFailureP(mr.fFailureP),
148 fNonRecons(mr.fNonRecons),
149 fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
150 fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
151 fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
152 fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
153 fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
154 fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
155 fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
156 fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
157 fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
158 fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
159 fhetaClustersLay2(mr.fhetaClustersLay2),
160 fhphiClustersLay2(mr.fhphiClustersLay2)
165 //______________________________________________________________________
166 AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
167 // Assignment operator
168 this->~AliITSTrackleterSPDEff();
169 new(this) AliITSTrackleterSPDEff(mr);
172 //______________________________________________________________________
173 AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
179 delete [] fAssociationFlag1;
181 delete [] fChipPredOnLay2;
182 delete [] fChipPredOnLay1;
184 delete [] fChipUpdatedInEvent;
186 delete [] fPredictionPrimary;
187 delete [] fPredictionSecondary;
188 delete [] fClusterPrimary;
189 delete [] fClusterSecondary;
190 delete [] fSuccessPP;
191 delete [] fSuccessTT;
197 delete [] fNonRecons;
202 //____________________________________________________________________
204 AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t*, AliStack *pStack, TTree *tRef) {
206 // - calls LoadClusterArray that finds the position of the clusters
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
219 // loading the clusters
220 LoadClusterArrays(clusterTree);
221 // to study residual background (i.e. contribution from TT' to measured efficiency)
222 if(fReflectClusterAroundZAxisForLayer0) ReflectClusterAroundZAxisForLayer(0);
223 if(fReflectClusterAroundZAxisForLayer1) ReflectClusterAroundZAxisForLayer(1);
225 if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
226 if(fMC && !tRef) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
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;
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;
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]));
241 //###########################################################
242 // Loop on layer 1 : finding theta, phi and z
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];
249 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
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
255 // find the Radius and the chip corresponding to the extrapolation point
257 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
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
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;
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]);
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];
281 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
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
287 // find the Radius and the chip corresponding to the extrapolation point
289 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
291 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
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;
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]);
308 //###########################################################
310 // First part : Extrapolation to Layer 2
313 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
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;
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!
326 // in any case, if MC has been required, store statistics of primaries and secondaries
327 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
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)); }
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]++;
343 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
344 } else { // this cluster is from a secondary particle
345 fClusterSecondary[key]++;
347 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
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]]++;
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]]++;
365 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
367 // The following excludes double associations
368 if (!fAssociationFlag[iC2]) {
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;
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];
382 fhClustersDPhiAll->Fill(dPhi);
383 fhClustersDThetaAll->Fill(dTheta);
384 fhClustersDZetaAll->Fill(dZeta);
385 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
386 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
389 // make "elliptical" cut in Phi and Zeta!
390 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
391 dZeta*dZeta/fZetaWindow/fZetaWindow);
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));
401 iC2WithBestDist = iC2;
404 } // end of loop over clusters in layer 2
406 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
409 fhClustersDPhiAcc->Fill(dPhimin);
410 fhClustersDThetaAcc->Fill(dThetamin);
411 fhClustersDZetaAcc->Fill(dZetamin);
412 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
413 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
416 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
417 // flag the association
419 // store the tracklet
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];
433 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
445 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
449 fTracklets[fNTracklets][3] = -2;
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]);
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]);
464 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
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];
476 if(primary) fSuccessPP[key]++;
478 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
479 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
482 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
484 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
485 fChipUpdatedInEvent[key]=kTRUE;
487 if(primary) fSuccessP[key]++;
488 if(secondary) fSuccessS[key]++;
492 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
493 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
495 if(primary) fSuccessP[key]++;
496 if(secondary) fSuccessS[key]++;
503 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
504 else if (fChipPredOnLay2[iC1]<1200) {
505 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
506 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
508 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
509 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
512 } // end of loop over clusters in layer 1
514 fNTracklets1=fNTracklets;
516 //###################################################################
518 // Second part : Interpolation to Layer 1
521 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
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;
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!
534 // in any case, if MC has been required, store statistics of primaries and secondaries
535 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
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)); }
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]++;
551 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
552 } else { // this cluster is from a secondary particle
553 fClusterSecondary[key]++;
555 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
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]]++;
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]]++;
573 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
575 // The following excludes double associations
576 if (!fAssociationFlag1[iC1]) {
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;
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];
592 fhClustersDPhiInterpAll->Fill(dPhi);
593 fhClustersDThetaInterpAll->Fill(dTheta);
594 fhClustersDZetaInterpAll->Fill(dZeta);
595 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
596 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
598 // make "elliptical" cut in Phi and Zeta!
599 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
600 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
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));
610 iC1WithBestDist = iC1;
613 } // end of loop over clusters in layer 1
615 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
618 fhClustersDPhiInterpAcc->Fill(dPhimin);
619 fhClustersDThetaInterpAcc->Fill(dThetamin);
620 fhClustersDZetaInterpAcc->Fill(dZetamin);
621 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
622 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
625 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
626 // flag the association
628 // store the tracklet
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];
642 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
654 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
658 fTracklets[fNTracklets][3] = -2;
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]);
665 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
668 nfClu1+=(Int_t)found; // this for debugging purpose
669 ntClu1++; // to check efficiency of the method FindChip
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];
675 if (label2 < 3) { // same label
677 if(primary) fSuccessPP[key]++;
679 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
680 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
683 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
685 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
686 fChipUpdatedInEvent[key]=kTRUE;
688 if(primary) fSuccessP[key]++;
689 if(secondary) fSuccessS[key]++;
692 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
693 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
695 if(primary) fSuccessP[key]++;
696 if(secondary) fSuccessS[key]++;
703 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
704 else if (fChipPredOnLay1[iC2]<1200) {
705 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
706 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
708 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
709 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
712 } // end of loop over clusters in layer 2
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));
720 //____________________________________________________________________
721 Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
722 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
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
729 // Output: Unique key to locate a chip
730 // return: kTRUE if succesfull
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));
736 // set phiVtx in the range [0,2pi]
737 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
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
757 // now you need to locate the chip within the idet detector,
758 // starting from the local coordinates in such a detector
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;
764 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
767 //______________________________________________________________________________
768 Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
770 // Return the average radius of a layer from Geometry
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 !
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);
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);
788 //______________________________________________________________________________
789 Bool_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
794 //Compute Cartesian Global Coordinate
795 Double_t xyzGlob[3],xyzLoc[3];
797 xyzGlob[0]=r*TMath::Cos(phi);
798 xyzGlob[1]=r*TMath::Sin(phi);
803 if(idet<0) return kFALSE;
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;
809 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
811 xloc = (Float_t)xyzLoc[0];
812 zloc = (Float_t)xyzLoc[2];
816 //______________________________________________________________________________
817 Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
818 //--------------------------------------------------------------------
819 // This function finds the detector crossed by the track
820 // Input: layer in range [0,1]
821 // phi in ALICE absolute reference system
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);
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);
835 if (zOffset<0) // old geometry
836 dphi = -(phi-phiOffset);
838 dphi = phi-phiOffset;
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;
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;}
852 return np*ndetectors + nz;
854 //____________________________________________________________
855 Bool_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.
859 Input: - 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
862 Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
864 Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
865 r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
866 In the same system, the equation of a semi-line is: phi=phiVtx;
867 Hence you get one interception only: P=(r,phiVtx)
868 Finally you want P in the ABSOLUTE ALICE reference system.
870 Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
871 Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
872 Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
873 Double_t cC=rO*rO-R*R;
874 Double_t dDelta=bB*bB-4*cC;
875 if(dDelta<0) return kFALSE;
877 r1=(-bB-TMath::Sqrt(dDelta))/2;
878 r2=(-bB+TMath::Sqrt(dDelta))/2;
879 if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
880 Double_t r=TMath::Max(r1,r2); // take the positive
881 Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
882 Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
883 pvtx[0]=r*TMath::Cos(phiVtx);
884 pvtx[1]=r*TMath::Sin(phiVtx);
885 pP[0]=vtx[0]+pvtx[0];
886 pP[1]=vtx[1]+pvtx[1];
887 phi=TMath::ATan2(pP[1],pP[0]);
890 //___________________________________________________________
891 Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
893 // simple method to reduce all angles (in rad)
897 while(angle >=2*TMath::Pi() || angle<0) {
898 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
899 if(angle < 0) angle+=2*TMath::Pi();
903 //___________________________________________________________
904 Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
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".
915 if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
916 if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
917 if(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() ));
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;}
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);
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
940 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
943 //_____________________________________________________________________________________________
944 Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
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).
949 // It define "detectable" a primary particle according to the following criteria:
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).
957 if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
958 if(!stack) {AliError("null pointer to MC stack"); return 0;}
959 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
961 TParticle* part = stack->Particle(ipart);
962 //TParticle* part0 = stack->Particle(0); // first primary
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
971 Int_t lastDau = part->GetLastDaughter();
972 nDau = lastDau - firstDau + 1;
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
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)
991 else {// not ele or muon
992 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
996 //_________________________________________________________________
997 void AliITSTrackleterSPDEff::InitPredictionMC() {
999 // this method allocate memory for the MC related informations
1000 // all the counters are set to 0
1003 if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1004 fPredictionPrimary = new Int_t[1200];
1005 fPredictionSecondary = new Int_t[1200];
1006 fClusterPrimary = new Int_t[1200];
1007 fClusterSecondary = new Int_t[1200];
1008 fSuccessPP = new Int_t[1200];
1009 fSuccessTT = new Int_t[1200];
1010 fSuccessS = new Int_t[1200];
1011 fSuccessP = new Int_t[1200];
1012 fFailureS = new Int_t[1200];
1013 fFailureP = new Int_t[1200];
1014 fRecons = new Int_t[1200];
1015 fNonRecons = new Int_t[1200];
1016 for(Int_t i=0; i<1200; i++) {
1017 fPredictionPrimary[i]=0;
1018 fPredictionSecondary[i]=0;
1019 fPredictionSecondary[i]=0;
1020 fClusterSecondary[i]=0;
1032 //_________________________________________________________________
1033 void AliITSTrackleterSPDEff::DeletePredictionMC() {
1035 // this method deallocate memory for the MC related informations
1036 // all the counters are set to 0
1039 if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1040 if(fPredictionPrimary) {
1041 delete fPredictionPrimary; fPredictionPrimary=0;
1043 if(fPredictionSecondary) {
1044 delete fPredictionSecondary; fPredictionSecondary=0;
1046 if(fClusterPrimary) {
1047 delete fClusterPrimary; fClusterPrimary=0;
1049 if(fClusterSecondary) {
1050 delete fClusterSecondary; fClusterSecondary=0;
1053 delete fSuccessPP; fSuccessPP=0;
1056 delete fSuccessTT; fSuccessTT=0;
1059 delete fSuccessS; fSuccessS=0;
1062 delete fSuccessP; fSuccessP=0;
1065 delete fFailureS; fFailureS=0;
1068 delete fFailureP; fFailureP=0;
1071 delete fRecons; fRecons=0;
1074 delete fNonRecons; fNonRecons=0;
1078 //______________________________________________________________________
1079 Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
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
1089 if (!fMC) {CallWarningMC(); return 0;}
1090 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1091 return fPredictionPrimary[(Int_t)key];
1093 //______________________________________________________________________
1094 Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
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
1104 if (!fMC) {CallWarningMC(); return 0;}
1105 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1106 return fPredictionSecondary[(Int_t)key];
1108 //______________________________________________________________________
1109 Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
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
1119 if (!fMC) {CallWarningMC(); return 0;}
1120 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1121 return fClusterPrimary[(Int_t)key];
1123 //______________________________________________________________________
1124 Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
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
1133 if (!fMC) {CallWarningMC(); return 0;}
1134 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1135 return fClusterSecondary[(Int_t)key];
1137 //______________________________________________________________________
1138 Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
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
1146 if (!fMC) {CallWarningMC(); return 0;}
1147 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1148 return fSuccessPP[(Int_t)key];
1150 //______________________________________________________________________
1151 Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
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
1159 if (!fMC) {CallWarningMC(); return 0;}
1160 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1161 return fSuccessTT[(Int_t)key];
1163 //______________________________________________________________________
1164 Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
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
1172 if (!fMC) {CallWarningMC(); return 0;}
1173 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1174 return fSuccessS[(Int_t)key];
1176 //______________________________________________________________________
1177 Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
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
1185 if (!fMC) {CallWarningMC(); return 0;}
1186 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1187 return fSuccessP[(Int_t)key];
1189 //______________________________________________________________________
1190 Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
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
1198 if (!fMC) {CallWarningMC(); return 0;}
1199 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1200 return fFailureS[(Int_t)key];
1202 //______________________________________________________________________
1203 Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
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
1211 if (!fMC) {CallWarningMC(); return 0;}
1212 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1213 return fFailureP[(Int_t)key];
1215 //_____________________________________________________________________
1216 Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
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
1224 if (!fMC) {CallWarningMC(); return 0;}
1225 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1226 return fRecons[(Int_t)key];
1228 //_____________________________________________________________________
1229 Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
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
1237 if (!fMC) {CallWarningMC(); return 0;}
1238 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1239 return fNonRecons[(Int_t)key];
1241 //______________________________________________________________________
1242 void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1243 // Print out some class data values in Ascii Form to output stream
1245 // ostream *os Output stream where Ascii data is to be writen
1250 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow
1251 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1252 << " " << fUpdateOncePerEventPlaneEff << " " << fReflectClusterAroundZAxisForLayer0
1253 << " " << fReflectClusterAroundZAxisForLayer1;
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) ;
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) ;
1273 //______________________________________________________________________
1274 void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1275 // Read in some class data values in Ascii Form to output stream
1277 // istream *is Input stream where Ascii data is to be read in from
1284 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow
1285 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1286 >> fUpdateOncePerEventPlaneEff >> fReflectClusterAroundZAxisForLayer0
1287 >> fReflectClusterAroundZAxisForLayer1;
1288 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
1290 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
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] ;
1311 //______________________________________________________________________
1312 ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1313 // Standard output streaming function
1315 // ostream &os output steam
1316 // AliITSTrackleterSPDEff &s class to be streamed.
1320 // ostream &os The stream pointer
1325 //______________________________________________________________________
1326 istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1327 // Standard inputput streaming function
1329 // istream &is input steam
1330 // AliITSTrackleterSPDEff &s class to be streamed.
1334 // ostream &os The stream pointer
1336 //printf("prova %d \n", (Int_t)s.GetMC());
1340 //______________________________________________________________________
1341 void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
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
1351 //if(!fMC) {CallWarningMC(); return;}
1352 if (!filename.Contains(".root")) {
1353 ofstream out(filename.Data(),ios::out | ios::binary);
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);
1373 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
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);
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)) ;
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)) ;
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)) ;
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)) ;
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)) ;
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)) ;
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)) ;
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)) ;
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)) ;
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)) ;
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)) ;
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)) ;
1426 //____________________________________________________________________
1427 void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
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
1436 //if(!fMC) {CallWarningMC(); return;}
1437 if( gSystem->AccessPathName( filename.Data() ) ) {
1438 AliError( Form( "file (%s) not found", filename.Data() ) );
1442 if (!filename.Contains(".root")) {
1443 ifstream in(filename.Data(),ios::in | ios::binary);
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);
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) ;
1501 //____________________________________________________________________
1502 Bool_t AliITSTrackleterSPDEff::SaveHists() {
1503 // This (private) method save the histograms on the output file
1504 // (only if fHistOn is TRUE).
1505 // Also the histograms from the base class are saved through the
1506 // AliITSMultReconstructor::SaveHists() call
1508 if (!GetHistOn()) return kFALSE;
1510 AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1512 fhClustersDPhiInterpAll->Write();
1513 fhClustersDThetaInterpAll->Write();
1514 fhClustersDZetaInterpAll->Write();
1515 fhDPhiVsDThetaInterpAll->Write();
1516 fhDPhiVsDZetaInterpAll->Write();
1518 fhClustersDPhiInterpAcc->Write();
1519 fhClustersDThetaInterpAcc->Write();
1520 fhClustersDZetaInterpAcc->Write();
1521 fhDPhiVsDThetaInterpAcc->Write();
1522 fhDPhiVsDZetaInterpAcc->Write();
1524 fhetaClustersLay2->Write();
1525 fhphiClustersLay2->Write();
1528 //__________________________________________________________
1529 Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1531 // Saves the histograms into a tree and saves the trees into a file
1532 // Also the histograms from the base class are saved
1534 if (!GetHistOn()) return kFALSE;
1535 if (filename.Data()=="") {
1536 AliWarning("WriteHistosToFile: null output filename!");
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;
1546 //____________________________________________________________
1547 void AliITSTrackleterSPDEff::BookHistos() {
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
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);
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);
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);
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);
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);
1585 //____________________________________________________________
1586 void AliITSTrackleterSPDEff::DeleteHistos() {
1588 // Private method to delete Histograms from memory
1589 // it is called. e.g., by the destructor.
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;}
1604 //_______________________________________________________________
1605 Bool_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
1609 // It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1611 // - Int_t layer (either 0 or 1): layer which you want to chech if the tracklete can be
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
1618 Bool_t ret=kFALSE; // returned value
1619 Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1620 if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1621 if(!stack) {AliError("null pointer to MC stack"); return ret;}
1622 if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1623 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1624 if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1626 AliTrackReference *tref=0x0;
1627 Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1628 Int_t nentries = (Int_t)ref->GetEntries();
1629 TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1630 TBranch *br = ref->GetBranch("TrackReferences");
1631 br->SetAddress(&tcaRef);
1632 for(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;}
1640 if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1641 br->GetEntry(imatch); // redundant, nevertheless ...
1642 Int_t nref=tcaRef->GetEntriesFast();
1643 for(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
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];
1654 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
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
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;
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];
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;}
1676 if(layer==0) { // try to see if it is reconstructable at the inner layer
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;
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];
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;};
1697 //_________________________________________________________________________
1698 void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
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
1707 if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1708 AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1710 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1711 fClustersLay1[iC1][0]*=-1;
1712 fClustersLay1[iC1][1]*=-1;
1716 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1717 fClustersLay2[iC2][0]*=-1;
1718 fClustersLay2[iC2][1]*=-1;