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>
35 #include <TParticlePDG.h>
37 #include <Riostream.h>
38 #include <TClonesArray.h>
40 #include "AliTracker.h"
41 #include "AliITSTrackleterSPDEff.h"
42 #include "AliITSgeomTGeo.h"
44 #include "AliITSPlaneEffSPD.h"
46 #include "AliTrackReference.h"
47 #include "AliRunLoader.h"
48 #include "AliITSReconstructor.h"
49 #include "AliITSRecPoint.h"
50 //____________________________________________________________________
51 ClassImp(AliITSTrackleterSPDEff)
54 //____________________________________________________________________
55 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
65 fOnlyOneTrackletPerC2(0),
72 fhClustersDThetaAcc(0),
73 fhClustersDZetaAcc(0),
75 fhClustersDThetaAll(0),
76 fhClustersDZetaAll(0),
92 fOnlyOneTrackletPerC1(0),
93 fUpdateOncePerEventPlaneEff(0),
94 fChipUpdatedInEvent(0),
97 fReflectClusterAroundZAxisForLayer0(kFALSE),
98 fReflectClusterAroundZAxisForLayer1(kFALSE),
99 fLightBkgStudyInParallel(kFALSE),
101 fUseOnlyPrimaryForPred(0),
102 fUseOnlySecondaryForPred(0),
103 fUseOnlySameParticle(0),
104 fUseOnlyDifferentParticle(0),
105 fUseOnlyStableParticle(1),
106 fPredictionPrimary(0),
107 fPredictionSecondary(0),
109 fClusterSecondary(0),
118 fhClustersDPhiInterpAcc(0),
119 fhClustersDThetaInterpAcc(0),
120 fhClustersDZetaInterpAcc(0),
121 fhClustersDPhiInterpAll(0),
122 fhClustersDThetaInterpAll(0),
123 fhClustersDZetaInterpAll(0),
124 fhDPhiVsDThetaInterpAll(0),
125 fhDPhiVsDThetaInterpAcc(0),
126 fhDPhiVsDZetaInterpAll(0),
127 fhDPhiVsDZetaInterpAcc(0),
128 fhetaClustersLay2(0),
129 fhphiClustersLay2(0),
131 fhClustersInModuleLay1(0),
132 fhClustersInModuleLay2(0)
134 // default constructor
135 // from AliITSMultReconstructor
138 SetOnlyOneTrackletPerC2();
139 fClustersLay1 = new Float_t*[300000];
140 fClustersLay2 = new Float_t*[300000];
141 fTracklets = new Float_t*[300000];
142 fAssociationFlag = new Bool_t[300000];
146 SetOnlyOneTrackletPerC1();
148 fAssociationFlag1 = new Bool_t[300000];
149 fChipPredOnLay2 = new UInt_t[300000];
150 fChipPredOnLay1 = new UInt_t[300000];
151 fChipUpdatedInEvent = new Bool_t[1200];
153 for(Int_t i=0; i<300000; i++) {
154 // from AliITSMultReconstructor
155 fClustersLay1[i] = new Float_t[6];
156 fClustersLay2[i] = new Float_t[6];
157 fTracklets[i] = new Float_t[5];
158 fAssociationFlag[i] = kFALSE;
160 fAssociationFlag1[i] = kFALSE;
162 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
164 if (GetHistOn()) BookHistos();
166 fPlaneEffSPD = new AliITSPlaneEffSPD();
167 SetLightBkgStudyInParallel();
169 //______________________________________________________________________
170 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
172 // from AliITSMultReconstructor
173 fClustersLay1(mr.fClustersLay1),
174 fClustersLay2(mr.fClustersLay2),
175 fTracklets(mr.fTracklets),
176 fAssociationFlag(mr.fAssociationFlag),
177 fNClustersLay1(mr.fNClustersLay1),
178 fNClustersLay2(mr.fNClustersLay2),
179 fNTracklets(mr.fNTracklets),
180 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
181 fPhiWindowL2(mr.fPhiWindowL2),
182 fZetaWindowL2(mr.fZetaWindowL2),
183 fPhiOverlapCut(mr.fPhiOverlapCut),
184 fZetaOverlapCut(mr.fZetaOverlapCut),
186 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
187 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
188 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
189 fhClustersDPhiAll(mr.fhClustersDPhiAll),
190 fhClustersDThetaAll(mr.fhClustersDThetaAll),
191 fhClustersDZetaAll(mr.fhClustersDZetaAll),
192 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
193 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
194 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
195 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
196 fhetaTracklets(mr.fhetaTracklets),
197 fhphiTracklets(mr.fhphiTracklets),
198 fhetaClustersLay1(mr.fhetaClustersLay1),
199 fhphiClustersLay1(mr.fhphiClustersLay1),
201 fAssociationFlag1(mr.fAssociationFlag1),
202 fChipPredOnLay2(mr.fChipPredOnLay2),
203 fChipPredOnLay1(mr.fChipPredOnLay1),
204 fNTracklets1(mr.fNTracklets1),
205 fPhiWindowL1(mr.fPhiWindowL1),
206 fZetaWindowL1(mr.fZetaWindowL1),
207 fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
208 fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
209 fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
210 fPlaneEffSPD(mr.fPlaneEffSPD),
211 fPlaneEffBkg(mr.fPlaneEffBkg),
212 fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
213 fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
214 fLightBkgStudyInParallel(mr.fLightBkgStudyInParallel),
216 fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
217 fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
218 fUseOnlySameParticle(mr.fUseOnlySameParticle),
219 fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
220 fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
221 fPredictionPrimary(mr.fPredictionPrimary),
222 fPredictionSecondary(mr.fPredictionSecondary),
223 fClusterPrimary(mr.fClusterPrimary),
224 fClusterSecondary(mr.fClusterSecondary),
225 fSuccessPP(mr.fSuccessPP),
226 fSuccessTT(mr.fSuccessTT),
227 fSuccessS(mr.fSuccessS),
228 fSuccessP(mr.fSuccessP),
229 fFailureS(mr.fFailureS),
230 fFailureP(mr.fFailureP),
232 fNonRecons(mr.fNonRecons),
233 fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
234 fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
235 fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
236 fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
237 fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
238 fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
239 fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
240 fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
241 fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
242 fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
243 fhetaClustersLay2(mr.fhetaClustersLay2),
244 fhphiClustersLay2(mr.fhphiClustersLay2),
245 fhClustersInChip(mr.fhClustersInChip),
246 fhClustersInModuleLay1(mr.fhClustersInModuleLay1),
247 fhClustersInModuleLay2(mr.fhClustersInModuleLay2)
252 //______________________________________________________________________
253 AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
254 // Assignment operator
255 this->~AliITSTrackleterSPDEff();
256 new(this) AliITSTrackleterSPDEff(mr);
259 //______________________________________________________________________
260 AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
262 // from AliITSMultReconstructor
264 for(Int_t i=0; i<300000; i++) {
265 delete [] fClustersLay1[i];
266 delete [] fClustersLay2[i];
267 delete [] fTracklets[i];
269 delete [] fClustersLay1;
270 delete [] fClustersLay2;
271 delete [] fTracklets;
272 delete [] fAssociationFlag;
277 delete [] fAssociationFlag1;
279 delete [] fChipPredOnLay2;
280 delete [] fChipPredOnLay1;
282 delete [] fChipUpdatedInEvent;
284 delete [] fPredictionPrimary;
285 delete [] fPredictionSecondary;
286 delete [] fClusterPrimary;
287 delete [] fClusterSecondary;
288 delete [] fSuccessPP;
289 delete [] fSuccessTT;
295 delete [] fNonRecons;
306 //____________________________________________________________________
308 AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef, Bool_t lbkg) {
310 // - you have to take care of the following, before of using Reconstruct
311 // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
312 // and convert the cluster coordinates to theta, phi (seen from the
313 // interaction vertex).
314 // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
315 // - Find the extrapolation/interpolation point.
316 // - Find the chip corresponding to that
317 // - Check if there is a cluster near that point
320 if(lbkg && !GetLightBkgStudyInParallel()) {
321 AliError("You asked for lightBackground in the Reconstruction without proper call to SetLightBkgStudyInParallel(1)");
324 AliITSPlaneEffSPD *pe;
331 // retrieve the vertex position
333 vtx[0]=(Float_t)GetX();
334 vtx[1]=(Float_t)GetY();
335 vtx[2]=(Float_t)GetZ();
336 // to study residual background (i.e. contribution from TT' to measured efficiency)
337 if(fReflectClusterAroundZAxisForLayer0 && !lbkg) ReflectClusterAroundZAxisForLayer(0);
338 if(fReflectClusterAroundZAxisForLayer1 && !lbkg) ReflectClusterAroundZAxisForLayer(1);
340 if(fMC && !pStack && !lbkg) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
341 if(fMC && !tRef && !lbkg) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
343 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
344 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
345 Int_t nfClu1=0; Int_t ntClu1=0;
346 Int_t nfClu2=0; Int_t ntClu2=0;
348 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
349 // for this new event)
350 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
352 // find the tracklets
353 AliDebug(1,"Looking for tracklets... ");
354 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
356 //###########################################################
357 // Loop on layer 1 : finding theta, phi and z
359 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
360 Float_t x = fClustersLay1[iC1][0] - vtx[0];
361 Float_t y = fClustersLay1[iC1][1] - vtx[1];
362 Float_t z = fClustersLay1[iC1][2] - vtx[2];
364 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
366 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
367 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
368 fClustersLay1[iC1][2] = z; // Store z
370 // find the Radius and the chip corresponding to the extrapolation point
372 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
374 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
375 key=999999; // also some other actions should be taken if not Found
377 nfTraPred2+=(Int_t)found; // this for debugging purpose
378 ntTraPred2++; // to check efficiency of the method FindChip
379 fChipPredOnLay2[iC1] = key;
380 fAssociationFlag1[iC1] = kFALSE;
382 if (fHistOn && !lbkg) {
383 Float_t eta=fClustersLay1[iC1][0];
384 eta= TMath::Tan(eta/2.);
385 eta=-TMath::Log(eta);
386 fhetaClustersLay1->Fill(eta);
387 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
388 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
391 // Loop on layer 2 : finding theta, phi and r
392 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
393 Float_t x = fClustersLay2[iC2][0] - vtx[0];
394 Float_t y = fClustersLay2[iC2][1] - vtx[1];
395 Float_t z = fClustersLay2[iC2][2] - vtx[2];
397 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
399 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
400 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
401 fClustersLay2[iC2][2] = z; // Store z
403 // find the Radius and the chip corresponding to the extrapolation point
405 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
407 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
410 nfTraPred1+=(Int_t)found; // this for debugging purpose
411 ntTraPred1++; // to check efficiency of the method FindChip
412 fChipPredOnLay1[iC2] = key;
413 fAssociationFlag[iC2] = kFALSE;
415 if (fHistOn && !lbkg) {
416 Float_t eta=fClustersLay2[iC2][0];
417 eta= TMath::Tan(eta/2.);
418 eta=-TMath::Log(eta);
419 fhetaClustersLay2->Fill(eta);
420 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
421 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
425 //###########################################################
427 // First part : Extrapolation to Layer 2
430 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
432 // here the control to check whether the efficiency of the chip traversed by this tracklet
433 // prediction has already been updated in this event using another tracklet prediction
434 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
436 // reset of variables for multiple candidates
437 Int_t iC2WithBestDist = 0; // reset
438 Float_t distmin = 100.; // just to put a huge number!
439 Float_t dPhimin = 0.; // Used for histograms only!
440 Float_t dThetamin = 0.; // Used for histograms only!
441 Float_t dZetamin = 0.; // Used for histograms only!
443 // in any case, if MC has been required, store statistics of primaries and secondaries
444 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
446 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
447 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
448 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
449 // do it always as a function of the chip number used to built the prediction
450 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
451 if (!found) {AliWarning(
452 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
454 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
455 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
456 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
457 { // this cluster is from a primary particle
458 fClusterPrimary[key]++;
460 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
461 } else { // this cluster is from a secondary particle
462 fClusterSecondary[key]++;
464 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
467 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
468 // (in case the prediction is reliable)
469 if( fChipPredOnLay2[iC1]<1200) {
470 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
471 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
472 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
473 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
474 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
475 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
476 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
477 else fNonRecons[fChipPredOnLay2[iC1]]++;
482 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
484 // The following excludes double associations
485 if (!fAssociationFlag[iC2]) {
487 // find the difference in angles
488 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
489 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
490 // take into account boundary condition
491 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
493 // find the difference in z (between linear projection from layer 1
494 // and the actual point: Dzeta= z1/r1*r2 -z2)
495 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
496 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
498 if (fHistOn && !lbkg) {
499 fhClustersDPhiAll->Fill(dPhi);
500 fhClustersDThetaAll->Fill(dTheta);
501 fhClustersDZetaAll->Fill(dZeta);
502 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
503 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
506 // make "elliptical" cut in Phi and Zeta!
507 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
508 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
512 //look for the minimum distance: the minimum is in iC2WithBestDist
513 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
514 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
518 iC2WithBestDist = iC2;
521 } // end of loop over clusters in layer 2
523 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
525 if (fHistOn && !lbkg) {
526 fhClustersDPhiAcc->Fill(dPhimin);
527 fhClustersDThetaAcc->Fill(dThetamin);
528 fhClustersDZetaAcc->Fill(dZetamin);
529 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
530 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
533 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
534 // flag the association
536 // store the tracklet
538 // use the theta from the clusters in the first layer
539 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
540 // use the phi from the clusters in the first layer
541 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
542 // Store the difference between phi1 and phi2
543 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
550 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
562 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
566 fTracklets[fNTracklets][3] = -2;
569 if (fHistOn && !lbkg) {
570 Float_t eta=fTracklets[fNTracklets][0];
571 eta= TMath::Tan(eta/2.);
572 eta=-TMath::Log(eta);
573 fhetaTracklets->Fill(eta);
574 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
577 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
578 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
581 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
584 nfClu2+=(Int_t)found; // this for debugging purpose
585 ntClu2++; // to check efficiency of the method FindChip
586 if(key<1200) { // the Chip has been found
587 if(fMC && !lbkg) { // this part only for MC
588 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
589 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
590 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
593 if(primary) fSuccessPP[key]++;
595 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
596 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
599 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
601 pe->UpDatePlaneEff(kTRUE,key); // success
602 fChipUpdatedInEvent[key]=kTRUE;
604 if(primary) fSuccessP[key]++;
605 if(secondary) fSuccessS[key]++;
609 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
610 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
612 if(primary) fSuccessP[key]++;
613 if(secondary) fSuccessS[key]++;
620 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
621 else if (fChipPredOnLay2[iC1]<1200) {
622 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
623 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
625 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
626 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
629 } // end of loop over clusters in layer 1
631 fNTracklets1=fNTracklets;
633 //###################################################################
635 // Second part : Interpolation to Layer 1
638 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
640 // here the control to check whether the efficiency of the chip traversed by this tracklet
641 // prediction has already been updated in this event using another tracklet prediction
642 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
644 // reset of variables for multiple candidates
645 Int_t iC1WithBestDist = 0; // reset
646 Float_t distmin = 100.; // just to put a huge number!
647 Float_t dPhimin = 0.; // Used for histograms only!
648 Float_t dThetamin = 0.; // Used for histograms only!
649 Float_t dZetamin = 0.; // Used for histograms only!
651 // in any case, if MC has been required, store statistics of primaries and secondaries
652 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
654 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
655 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
656 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
657 // do it always as a function of the chip number used to built the prediction
658 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
659 if (!found) {AliWarning(
660 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
662 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
663 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
664 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
665 { // this cluster is from a primary particle
666 fClusterPrimary[key]++;
668 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
669 } else { // this cluster is from a secondary particle
670 fClusterSecondary[key]++;
672 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
675 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
676 // (in case the prediction is reliable)
677 if( fChipPredOnLay1[iC2]<1200) {
678 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
679 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
680 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
681 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
682 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
683 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
684 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
685 else fNonRecons[fChipPredOnLay1[iC2]]++;
690 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
692 // The following excludes double associations
693 if (!fAssociationFlag1[iC1]) {
695 // find the difference in angles
696 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
697 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
698 // take into account boundary condition
699 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
702 // find the difference in z (between linear projection from layer 2
703 // and the actual point: Dzeta= z2/r2*r1 -z1)
704 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
705 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
708 if (fHistOn && !lbkg) {
709 fhClustersDPhiInterpAll->Fill(dPhi);
710 fhClustersDThetaInterpAll->Fill(dTheta);
711 fhClustersDZetaInterpAll->Fill(dZeta);
712 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
713 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
715 // make "elliptical" cut in Phi and Zeta!
716 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
717 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
721 //look for the minimum distance: the minimum is in iC1WithBestDist
722 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
723 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
727 iC1WithBestDist = iC1;
730 } // end of loop over clusters in layer 1
732 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
734 if (fHistOn && !lbkg) {
735 fhClustersDPhiInterpAcc->Fill(dPhimin);
736 fhClustersDThetaInterpAcc->Fill(dThetamin);
737 fhClustersDZetaInterpAcc->Fill(dZetamin);
738 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
739 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
742 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
743 // flag the association
745 // store the tracklet
747 // use the theta from the clusters in the first layer
748 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
749 // use the phi from the clusters in the first layer
750 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
751 // Store the difference between phi1 and phi2
752 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
759 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
771 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
775 fTracklets[fNTracklets][3] = -2;
778 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
779 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
782 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
785 nfClu1+=(Int_t)found; // this for debugging purpose
786 ntClu1++; // to check efficiency of the method FindChip
788 if(fMC && !lbkg) { // this part only for MC
789 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
790 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
791 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
792 if (label2 < 3) { // same label
794 if(primary) fSuccessPP[key]++;
796 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
797 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
800 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
802 pe->UpDatePlaneEff(kTRUE,key); // success
803 fChipUpdatedInEvent[key]=kTRUE;
805 if(primary) fSuccessP[key]++;
806 if(secondary) fSuccessS[key]++;
809 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
810 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
812 if(primary) fSuccessP[key]++;
813 if(secondary) fSuccessS[key]++;
820 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
821 else if (fChipPredOnLay1[iC2]<1200) {
822 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
823 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
825 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
826 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
829 } // end of loop over clusters in layer 2
831 AliDebug(1,Form("%d tracklets found", fNTracklets));
832 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
833 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
834 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
835 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
837 //____________________________________________________________________
838 Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
839 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
841 // Input: a) layer number in the range [0,1]
842 // b) vtx[3]: actual vertex
843 // c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
844 // d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
846 // Output: Unique key to locate a chip
847 // return: kTRUE if succesfull
849 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
850 Double_t r=GetRLayer(layer);
851 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
853 // set phiVtx in the range [0,2pi]
854 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
856 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
857 // of the intersection of the tracklet with the pixel layer.
858 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
859 else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
860 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
861 Double_t vtxy[2]={vtx[0],vtx[1]};
862 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
863 // this method gives you two interceptions
864 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
865 // set phiAbs in the range [0,2pi]
866 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
867 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
868 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
869 // taking the closest one to phiVtx
870 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
871 } else phiAbs=phiVtx;
872 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
874 // now you need to locate the chip within the idet detector,
875 // starting from the local coordinates in such a detector
877 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
878 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
879 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
881 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
884 //______________________________________________________________________________
885 Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
887 // Return the average radius of a layer from Geometry
889 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
890 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
892 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
893 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
894 Double_t r=TMath::Sqrt(x*x + y*y);
896 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
897 r += TMath::Sqrt(x*x + y*y);
898 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
899 r += TMath::Sqrt(x*x + y*y);
900 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
901 r += TMath::Sqrt(x*x + y*y);
905 //______________________________________________________________________________
906 Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
907 Float_t &xloc, Float_t &zloc) {
908 // this method transform Global Cilindrical coordinates into local (i.e. module)
909 // cartesian coordinates
911 //Compute Cartesian Global Coordinate
912 Double_t xyzGlob[3],xyzLoc[3];
914 xyzGlob[0]=r*TMath::Cos(phi);
915 xyzGlob[1]=r*TMath::Sin(phi);
920 if(idet<0) return kFALSE;
922 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
923 Int_t lad = Int_t(idet/ndet) + 1;
924 Int_t det = idet - (lad-1)*ndet + 1;
926 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
928 xloc = (Float_t)xyzLoc[0];
929 zloc = (Float_t)xyzLoc[2];
933 //______________________________________________________________________________
934 Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
935 //--------------------------------------------------------------------
936 // This function finds the detector crossed by the track
937 // Input: layer in range [0,1]
938 // phi in ALICE absolute reference system
940 //--------------------------------------------------------------------
941 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
942 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
943 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
944 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
946 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
947 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
948 Double_t phiOffset=TMath::ATan2(y,x);
952 if (zOffset<0) // old geometry
953 dphi = -(phi-phiOffset);
955 dphi = phi-phiOffset;
957 if (dphi < 0) dphi += 2*TMath::Pi();
958 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
959 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
960 if (np>=nladders) np-=nladders;
961 if (np<0) np+=nladders;
963 Double_t dz=zOffset-z;
964 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
965 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
966 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
967 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
969 return np*ndetectors + nz;
971 //____________________________________________________________
972 Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
973 // this method find the intersection in xy between a tracklet (straight line) and
974 // a circonference (r=R), using polar coordinates.
976 Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
977 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
978 - R: radius of the circle
979 Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
981 Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
982 r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
983 In the same system, the equation of a semi-line is: phi=phiVtx;
984 Hence you get one interception only: P=(r,phiVtx)
985 Finally you want P in the ABSOLUTE ALICE reference system.
987 Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
988 Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
989 Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
990 Double_t cC=rO*rO-R*R;
991 Double_t dDelta=bB*bB-4*cC;
992 if(dDelta<0) return kFALSE;
994 r1=(-bB-TMath::Sqrt(dDelta))/2;
995 r2=(-bB+TMath::Sqrt(dDelta))/2;
996 if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
997 Double_t r=TMath::Max(r1,r2); // take the positive
998 Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
999 Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
1000 pvtx[0]=r*TMath::Cos(phiVtx);
1001 pvtx[1]=r*TMath::Sin(phiVtx);
1002 pP[0]=vtx[0]+pvtx[0];
1003 pP[1]=vtx[1]+pvtx[1];
1004 phi=TMath::ATan2(pP[1],pP[0]);
1007 //___________________________________________________________
1008 Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
1010 // simple method to reduce all angles (in rad)
1014 while(angle >=2*TMath::Pi() || angle<0) {
1015 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
1016 if(angle < 0) angle+=2*TMath::Pi();
1020 //___________________________________________________________
1021 Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
1023 // This method check if a particle is primary; i.e.
1024 // it comes from the main vertex and it is a "stable" particle, according to
1025 // AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
1026 // a stable particle: it has no effect on this analysis).
1027 // This method can be called only for MC events, where Kinematics is available.
1028 // if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
1029 // returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
1030 // The latter (see below) try to verify if a primary particle is also "detectable".
1032 if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
1033 if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
1034 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
1035 // return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
1036 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
1037 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
1038 TParticle* part = stack->Particle(ipart);
1039 TParticle* part0 = stack->Particle(0); // first primary
1040 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
1041 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
1042 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
1044 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
1045 TParticlePDG* pdgPart = part->GetPDG();
1046 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
1048 Double_t distx = part->Vx() - part0->Vx();
1049 Double_t disty = part->Vy() - part0->Vy();
1050 Double_t distz = part->Vz() - part0->Vz();
1051 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
1053 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
1054 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
1055 return kFALSE; }// primary if within 500 microns from true Vertex
1057 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
1060 //_____________________________________________________________________________________________
1061 Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
1063 // This private method can be applied on MC particles (if stack is available),
1064 // provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
1066 // It define "detectable" a primary particle according to the following criteria:
1068 // - if no decay products can be found in the stack (note that this does not
1069 // means it is stable, since a particle is stored in stack if it has at least 1 hit in a
1070 // sensitive detector)
1071 // - if it has at least one decay daughter produced outside or just on the outer pixel layer
1072 // - if the last decay particle is an electron (or a muon) which is not produced in-between
1073 // the two pixel layers (this is likely to be a kink).
1074 if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
1075 if(!stack) {AliError("null pointer to MC stack"); return 0;}
1076 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
1078 TParticle* part = stack->Particle(ipart);
1079 //TParticle* part0 = stack->Particle(0); // first primary
1085 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
1086 // its real fate ! But you have to take it !
1087 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
1088 Int_t lastDau = part->GetLastDaughter();
1089 nDau = lastDau - firstDau + 1;
1090 Double_t distMax=0.;
1092 for(Int_t j=firstDau; j<=lastDau; j++) {
1093 dau = stack->Particle(j);
1094 Double_t distx = dau->Vx();
1095 Double_t disty = dau->Vy();
1096 //Double_t distz = dau->Vz();
1097 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
1098 if(distR<distMax) continue; // considere only the daughter produced at largest radius
1102 dau = stack->Particle(jmax);
1103 pdgDau=dau->GetPdgCode();
1104 if (pdgDau == 11 || pdgDau == 13 ) {
1105 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
1106 else nret =0; // delta-ray emission in material (keep it)
1108 else {// not ele or muon
1109 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
1113 //_________________________________________________________________
1114 void AliITSTrackleterSPDEff::InitPredictionMC() {
1116 // this method allocate memory for the MC related informations
1117 // all the counters are set to 0
1120 if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1121 fPredictionPrimary = new Int_t[1200];
1122 fPredictionSecondary = new Int_t[1200];
1123 fClusterPrimary = new Int_t[1200];
1124 fClusterSecondary = new Int_t[1200];
1125 fSuccessPP = new Int_t[1200];
1126 fSuccessTT = new Int_t[1200];
1127 fSuccessS = new Int_t[1200];
1128 fSuccessP = new Int_t[1200];
1129 fFailureS = new Int_t[1200];
1130 fFailureP = new Int_t[1200];
1131 fRecons = new Int_t[1200];
1132 fNonRecons = new Int_t[1200];
1133 for(Int_t i=0; i<1200; i++) {
1134 fPredictionPrimary[i]=0;
1135 fPredictionSecondary[i]=0;
1136 fPredictionSecondary[i]=0;
1137 fClusterSecondary[i]=0;
1149 //_________________________________________________________________
1150 void AliITSTrackleterSPDEff::DeletePredictionMC() {
1152 // this method deallocate memory for the MC related informations
1153 // all the counters are set to 0
1156 if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1157 if(fPredictionPrimary) {
1158 delete fPredictionPrimary; fPredictionPrimary=0;
1160 if(fPredictionSecondary) {
1161 delete fPredictionSecondary; fPredictionSecondary=0;
1163 if(fClusterPrimary) {
1164 delete fClusterPrimary; fClusterPrimary=0;
1166 if(fClusterSecondary) {
1167 delete fClusterSecondary; fClusterSecondary=0;
1170 delete fSuccessPP; fSuccessPP=0;
1173 delete fSuccessTT; fSuccessTT=0;
1176 delete fSuccessS; fSuccessS=0;
1179 delete fSuccessP; fSuccessP=0;
1182 delete fFailureS; fFailureS=0;
1185 delete fFailureP; fFailureP=0;
1188 delete fRecons; fRecons=0;
1191 delete fNonRecons; fNonRecons=0;
1195 //______________________________________________________________________
1196 Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
1198 // This method return the Data menmber fPredictionPrimary [1200].
1199 // You can call it only for MC events.
1200 // fPredictionPrimary[key] contains the number of tracklet predictions on the
1201 // given chip key built using a cluster on the other layer produced (at least)
1202 // from a primary particle.
1203 // Key refers to the chip crossed by the prediction
1206 if (!fMC) {CallWarningMC(); return 0;}
1207 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1208 return fPredictionPrimary[(Int_t)key];
1210 //______________________________________________________________________
1211 Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
1213 // This method return the Data menmber fPredictionSecondary [1200].
1214 // You can call it only for MC events.
1215 // fPredictionSecondary[key] contains the number of tracklet predictions on the
1216 // given chip key built using a cluster on the other layer produced (only)
1217 // from a secondary particle
1218 // Key refers to the chip crossed by the prediction
1221 if (!fMC) {CallWarningMC(); return 0;}
1222 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1223 return fPredictionSecondary[(Int_t)key];
1225 //______________________________________________________________________
1226 Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
1228 // This method return the Data menmber fClusterPrimary [1200].
1229 // You can call it only for MC events.
1230 // fClusterPrimary[key] contains the number of tracklet predictions
1231 // built using a cluster on that layer produced (only)
1232 // from a primary particle
1233 // Key refers to the chip used to build the prediction
1236 if (!fMC) {CallWarningMC(); return 0;}
1237 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1238 return fClusterPrimary[(Int_t)key];
1240 //______________________________________________________________________
1241 Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
1243 // This method return the Data menmber fClusterSecondary [1200].
1244 // You can call it only for MC events.
1245 // fClusterSecondary[key] contains the number of tracklet predictions
1246 // built using a cluster on that layer produced (only)
1247 // from a secondary particle
1248 // Key refers to the chip used to build the prediction
1250 if (!fMC) {CallWarningMC(); return 0;}
1251 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1252 return fClusterSecondary[(Int_t)key];
1254 //______________________________________________________________________
1255 Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1257 // This method return the Data menmber fSuccessPP [1200].
1258 // You can call it only for MC events.
1259 // fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1260 // with a cluster on the other layer) built by using the same primary particle
1261 // the unique chip key refers to the chip which get updated its efficiency
1263 if (!fMC) {CallWarningMC(); return 0;}
1264 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1265 return fSuccessPP[(Int_t)key];
1267 //______________________________________________________________________
1268 Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1270 // This method return the Data menmber fSuccessTT [1200].
1271 // You can call it only for MC events.
1272 // fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1273 // with a cluster on the other layer) built by using the same particle (whatever)
1274 // the unique chip key refers to the chip which get updated its efficiency
1276 if (!fMC) {CallWarningMC(); return 0;}
1277 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1278 return fSuccessTT[(Int_t)key];
1280 //______________________________________________________________________
1281 Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1283 // This method return the Data menmber fSuccessS [1200].
1284 // You can call it only for MC events.
1285 // fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1286 // with a cluster on the other layer) built by using a secondary particle
1287 // the unique chip key refers to the chip which get updated its efficiency
1289 if (!fMC) {CallWarningMC(); return 0;}
1290 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1291 return fSuccessS[(Int_t)key];
1293 //______________________________________________________________________
1294 Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1296 // This method return the Data menmber fSuccessP [1200].
1297 // You can call it only for MC events.
1298 // fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1299 // with a cluster on the other layer) built by using a primary particle
1300 // the unique chip key refers to the chip which get updated its efficiency
1302 if (!fMC) {CallWarningMC(); return 0;}
1303 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1304 return fSuccessP[(Int_t)key];
1306 //______________________________________________________________________
1307 Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1309 // This method return the Data menmber fFailureS [1200].
1310 // You can call it only for MC events.
1311 // fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1312 // with a cluster on the other layer) built by using a secondary particle
1313 // the unique chip key refers to the chip which get updated its efficiency
1315 if (!fMC) {CallWarningMC(); return 0;}
1316 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1317 return fFailureS[(Int_t)key];
1319 //______________________________________________________________________
1320 Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1322 // This method return the Data menmber fFailureP [1200].
1323 // You can call it only for MC events.
1324 // fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1325 // with a cluster on the other layer) built by using a primary particle
1326 // the unique chip key refers to the chip which get updated its efficiency
1328 if (!fMC) {CallWarningMC(); return 0;}
1329 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1330 return fFailureP[(Int_t)key];
1332 //_____________________________________________________________________
1333 Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1335 // This method return the Data menmber fRecons [1200].
1336 // You can call it only for MC events.
1337 // fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1338 // has an hit in the detector)
1339 // the unique chip key refers to the chip where fall the prediction
1341 if (!fMC) {CallWarningMC(); return 0;}
1342 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1343 return fRecons[(Int_t)key];
1345 //_____________________________________________________________________
1346 Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1348 // This method return the Data menmber fNonRecons [1200].
1349 // You can call it only for MC events.
1350 // fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1351 // has not any hit in the detector)
1352 // the unique chip key refers to the chip where fall the prediction
1354 if (!fMC) {CallWarningMC(); return 0;}
1355 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1356 return fNonRecons[(Int_t)key];
1358 //______________________________________________________________________
1359 void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1360 // Print out some class data values in Ascii Form to output stream
1362 // ostream *os Output stream where Ascii data is to be writen
1367 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindowL2 <<" "<< fZetaWindowL2
1368 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1369 << " " << fUpdateOncePerEventPlaneEff << " " << fReflectClusterAroundZAxisForLayer0
1370 << " " << fReflectClusterAroundZAxisForLayer1;
1372 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1373 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1374 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1375 << " " << fUseOnlyStableParticle ;
1376 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1377 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1378 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1379 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
1380 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1381 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1382 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1383 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1384 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1385 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1386 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1387 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
1390 //______________________________________________________________________
1391 void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1392 // Read in some class data values in Ascii Form to output stream
1394 // istream *is Input stream where Ascii data is to be read in from
1401 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindowL2 >> fZetaWindowL2
1402 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1403 >> fUpdateOncePerEventPlaneEff >> fReflectClusterAroundZAxisForLayer0
1404 >> fReflectClusterAroundZAxisForLayer1;
1405 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
1407 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1409 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1410 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1411 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1412 >> fUseOnlyStableParticle;
1413 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1414 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1415 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1416 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1417 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1418 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1419 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1420 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1421 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1422 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1423 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1424 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1428 //______________________________________________________________________
1429 ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1430 // Standard output streaming function
1432 // ostream &os output steam
1433 // AliITSTrackleterSPDEff &s class to be streamed.
1437 // ostream &os The stream pointer
1442 //______________________________________________________________________
1443 istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1444 // Standard inputput streaming function
1446 // istream &is input steam
1447 // AliITSTrackleterSPDEff &s class to be streamed.
1451 // ostream &os The stream pointer
1453 //printf("prova %d \n", (Int_t)s.GetMC());
1457 //______________________________________________________________________
1458 void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
1460 // This Method write into an either asci or root file
1461 // the used cuts and the statistics of the MC related quantities
1462 // The method SetMC() has to be called before
1463 // Input TString filename: name of file for output (it deletes already existing
1468 //if(!fMC) {CallWarningMC(); return;}
1469 if (!filename.Contains(".root")) {
1470 ofstream out(filename.Data(),ios::out | ios::binary);
1476 TFile* mcfile = TFile::Open(filename, "RECREATE");
1477 TH1F* cuts = new TH1F("cuts", "list of cuts", 10, 0, 10); // TH1I containing cuts
1478 cuts->SetBinContent(1,fPhiWindowL1);
1479 cuts->SetBinContent(2,fZetaWindowL1);
1480 cuts->SetBinContent(3,fPhiWindowL2);
1481 cuts->SetBinContent(4,fZetaWindowL2);
1482 cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
1483 cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
1484 cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
1485 cuts->SetBinContent(8,fReflectClusterAroundZAxisForLayer0);
1486 cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer1);
1487 cuts->SetBinContent(10,fMC);
1490 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
1492 TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
1493 mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
1494 mc0->SetBinContent(2,fUseOnlySecondaryForPred);
1495 mc0->SetBinContent(3,fUseOnlySameParticle);
1496 mc0->SetBinContent(4,fUseOnlyDifferentParticle);
1497 mc0->SetBinContent(5,fUseOnlyStableParticle);
1501 mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
1502 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
1505 mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
1506 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
1509 mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
1510 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
1513 mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
1514 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
1517 mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
1518 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
1521 mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
1522 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
1525 mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
1526 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
1529 mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
1530 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
1533 mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
1534 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
1537 mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
1538 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
1541 mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
1542 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
1545 mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
1546 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
1554 //____________________________________________________________________
1555 void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
1557 // This Method read from an asci file (do not know why binary does not work)
1558 // the cuts to be used and the statistics of the MC related quantities
1559 // Input TString filename: name of input file for output
1560 // The method SetMC() has to be called before
1564 //if(!fMC) {CallWarningMC(); return;}
1565 if( gSystem->AccessPathName( filename.Data() ) ) {
1566 AliError( Form( "file (%s) not found", filename.Data() ) );
1570 if (!filename.Contains(".root")) {
1571 ifstream in(filename.Data(),ios::in | ios::binary);
1578 TFile *mcfile = TFile::Open(filename);
1579 TH1F *cuts = (TH1F*)mcfile->Get("cuts");
1580 fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
1581 fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
1582 fPhiWindowL2=(Float_t)cuts->GetBinContent(3);
1583 fZetaWindowL2=(Float_t)cuts->GetBinContent(4);
1584 fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
1585 fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
1586 fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
1587 fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(8);
1588 fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(9);
1589 fMC=(Bool_t)cuts->GetBinContent(10);
1590 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1591 else { // only if file with MC predictions
1592 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1593 TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
1594 fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
1595 fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
1596 fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
1597 fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
1598 fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
1600 mc1 =(TH1I*)mcfile->Get("mc1");
1601 for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1602 mc1 =(TH1I*)mcfile->Get("mc2");
1603 for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1604 mc1 =(TH1I*)mcfile->Get("mc3");
1605 for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1606 mc1 =(TH1I*)mcfile->Get("mc4");
1607 for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1608 mc1 =(TH1I*)mcfile->Get("mc5");
1609 for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1610 mc1 =(TH1I*)mcfile->Get("mc6");
1611 for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
1612 mc1 =(TH1I*)mcfile->Get("mc7");
1613 for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1614 mc1 =(TH1I*)mcfile->Get("mc8");
1615 for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1616 mc1 =(TH1I*)mcfile->Get("mc9");
1617 for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1618 mc1 =(TH1I*)mcfile->Get("mc10");
1619 for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1620 mc1 =(TH1I*)mcfile->Get("mc11");
1621 for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1622 mc1 =(TH1I*)mcfile->Get("mc12");
1623 for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1629 //____________________________________________________________________
1630 Bool_t AliITSTrackleterSPDEff::SaveHists() {
1631 // This (private) method save the histograms on the output file
1632 // (only if fHistOn is TRUE).
1633 // Also the histograms from the base class are saved through the
1634 // AliITSMultReconstructor::SaveHists() call
1636 if (!GetHistOn()) return kFALSE;
1638 // AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1639 fhClustersDPhiAll->Write();
1640 fhClustersDThetaAll->Write();
1641 fhClustersDZetaAll->Write();
1642 fhDPhiVsDThetaAll->Write();
1643 fhDPhiVsDZetaAll->Write();
1645 fhClustersDPhiAcc->Write();
1646 fhClustersDThetaAcc->Write();
1647 fhClustersDZetaAcc->Write();
1648 fhDPhiVsDThetaAcc->Write();
1649 fhDPhiVsDZetaAcc->Write();
1651 fhetaTracklets->Write();
1652 fhphiTracklets->Write();
1653 fhetaClustersLay1->Write();
1654 fhphiClustersLay1->Write();
1656 fhClustersDPhiInterpAll->Write();
1657 fhClustersDThetaInterpAll->Write();
1658 fhClustersDZetaInterpAll->Write();
1659 fhDPhiVsDThetaInterpAll->Write();
1660 fhDPhiVsDZetaInterpAll->Write();
1662 fhClustersDPhiInterpAcc->Write();
1663 fhClustersDThetaInterpAcc->Write();
1664 fhClustersDZetaInterpAcc->Write();
1665 fhDPhiVsDThetaInterpAcc->Write();
1666 fhDPhiVsDZetaInterpAcc->Write();
1668 fhetaClustersLay2->Write();
1669 fhphiClustersLay2->Write();
1670 fhClustersInChip->Write();
1671 for (Int_t nhist=0;nhist<80;nhist++){
1672 fhClustersInModuleLay1[nhist]->Write();
1674 for (Int_t nhist=0;nhist<160;nhist++){
1675 fhClustersInModuleLay2[nhist]->Write();
1679 //__________________________________________________________
1680 Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1682 // Saves the histograms into a tree and saves the trees into a file
1683 // Also the histograms from the base class are saved
1685 if (!GetHistOn()) return kFALSE;
1686 if (!strcmp(filename.Data(),"")) {
1687 AliWarning("WriteHistosToFile: null output filename!");
1690 TFile *hFile=new TFile(filename.Data(),option,
1691 "The File containing the histos for SPD efficiency studies with tracklets");
1692 if(!SaveHists()) return kFALSE;
1697 //____________________________________________________________
1698 void AliITSTrackleterSPDEff::BookHistos() {
1700 // This method books addtitional histograms
1701 // w.r.t. those of the base class.
1702 // In particular, the differences of cluster coordinate between the two SPD
1703 // layers are computed in the interpolation phase
1705 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1707 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
1708 fhClustersDPhiAcc->SetDirectory(0);
1709 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
1710 fhClustersDThetaAcc->SetDirectory(0);
1711 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
1712 fhClustersDZetaAcc->SetDirectory(0);
1714 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
1715 fhDPhiVsDZetaAcc->SetDirectory(0);
1716 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
1717 fhDPhiVsDThetaAcc->SetDirectory(0);
1719 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
1720 fhClustersDPhiAll->SetDirectory(0);
1721 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
1722 fhClustersDThetaAll->SetDirectory(0);
1723 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
1724 fhClustersDZetaAll->SetDirectory(0);
1726 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
1727 fhDPhiVsDZetaAll->SetDirectory(0);
1728 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
1729 fhDPhiVsDThetaAll->SetDirectory(0);
1731 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
1732 fhetaTracklets->SetDirectory(0);
1733 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
1734 fhphiTracklets->SetDirectory(0);
1735 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
1736 fhetaClustersLay1->SetDirectory(0);
1737 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
1738 fhphiClustersLay1->SetDirectory(0);
1740 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1741 fhClustersDPhiInterpAcc->SetDirectory(0);
1742 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1743 fhClustersDThetaInterpAcc->SetDirectory(0);
1744 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1745 fhClustersDZetaInterpAcc->SetDirectory(0);
1747 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1748 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1749 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1750 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1752 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1753 fhClustersDPhiInterpAll->SetDirectory(0);
1754 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1755 fhClustersDThetaInterpAll->SetDirectory(0);
1756 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1757 fhClustersDZetaInterpAll->SetDirectory(0);
1759 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1760 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1761 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1762 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1764 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1765 fhetaClustersLay2->SetDirectory(0);
1766 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1767 fhphiClustersLay2->SetDirectory(0);
1768 fhClustersInChip = new TH1F("fhClustersInChip", "ClustersPerChip", 1200, -0.5, 1199.5);
1769 fhClustersInChip->SetDirectory(0);
1770 // each chip is divided 8(z) x 4(y), i.e. in 32 squares, each containing 4 columns and 64 rows.
1771 Float_t bz[160]; const Float_t kconv = 1.0E-04; // converts microns to cm.
1772 for(Int_t i=0;i<160;i++) bz[i] = 425.0; // most are 425 microns except below
1773 bz[ 31] = bz[ 32] = 625.0; // first chip boundry
1774 bz[ 63] = bz[ 64] = 625.0; // first chip boundry
1775 bz[ 95] = bz[ 96] = 625.0; // first chip boundry
1776 bz[127] = bz[128] = 625.0; // first chip boundry
1777 Double_t xbins[41]; // each bin in x (Z loc coordinate) includes 4 columns
1779 Float_t xmn,xmx,zmn,zmx;
1780 if(!fPlaneEffSPD->GetBlockBoundaries(0,xmn,xmx,zmn,zmx)) AliWarning("Could not book histo properly");
1781 xbins[0]=(Double_t)zmn;
1782 for(Int_t i=0;i<40;i++) {
1783 xbins[i+1]=xbins[i] + (bz[4*i]+bz[4*i+1]+bz[4*i+2]+bz[4*i+3])*kconv;
1785 TString histname="ClustersLay1_mod_",aux;
1786 fhClustersInModuleLay1 =new TH2F*[80];
1787 for (Int_t nhist=0;nhist<80;nhist++){
1791 fhClustersInModuleLay1[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1792 fhClustersInModuleLay1[nhist]->SetName(aux.Data());
1793 fhClustersInModuleLay1[nhist]->SetTitle(aux.Data());
1794 fhClustersInModuleLay1[nhist]->SetDirectory(0);
1796 histname="ClustersLay2_mod_";
1797 fhClustersInModuleLay2 =new TH2F*[160];
1798 for (Int_t nhist=0;nhist<160;nhist++){
1801 fhClustersInModuleLay2[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1802 fhClustersInModuleLay2[nhist]->SetName(aux.Data());
1803 fhClustersInModuleLay2[nhist]->SetTitle(aux.Data());
1804 fhClustersInModuleLay2[nhist]->SetDirectory(0);
1809 //____________________________________________________________
1810 void AliITSTrackleterSPDEff::DeleteHistos() {
1812 // Private method to delete Histograms from memory
1813 // it is called. e.g., by the destructor.
1815 // form AliITSMultReconstructor
1816 if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
1817 if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
1818 if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
1819 if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
1820 if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
1821 if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
1822 if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
1823 if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
1824 if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
1825 if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
1826 if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
1827 if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
1828 if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
1829 if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
1831 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1832 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1833 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1834 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1835 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1836 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1837 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1838 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1839 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1840 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1841 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1842 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1843 if(fhClustersInChip) {delete fhClustersInChip; fhClustersInChip=0;}
1844 if(fhClustersInModuleLay1) {
1845 for (Int_t i=0; i<80; i++ ) delete fhClustersInModuleLay1[i];
1846 delete [] fhClustersInModuleLay1; fhClustersInModuleLay1=0;
1848 if(fhClustersInModuleLay2) {
1849 for (Int_t i=0; i<160; i++ ) delete fhClustersInModuleLay2[i];
1850 delete [] fhClustersInModuleLay2; fhClustersInModuleLay2=0;
1853 //_______________________________________________________________
1854 Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1855 Float_t* vtx, AliStack *stack, TTree *ref) {
1856 // This (private) method can be used only for MC events, where both AliStack and the TrackReference
1858 // It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1860 // - Int_t layer (either 0 or 1): layer which you want to chech if the tracklete can be
1862 // - Int_t iC : cluster index used to build the tracklet prediction
1863 // if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1864 // - Float_t* vtx: actual event vertex
1865 // - stack: pointer to Stack
1866 // - ref: pointer to TTRee of TrackReference
1867 Bool_t ret=kFALSE; // returned value
1868 Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1869 if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1870 if(!stack) {AliError("null pointer to MC stack"); return ret;}
1871 if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1872 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1873 if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1875 AliTrackReference *tref=0x0;
1876 Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1877 Int_t nentries = (Int_t)ref->GetEntries();
1878 TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1879 TBranch *br = ref->GetBranch("TrackReferences");
1880 br->SetAddress(&tcaRef);
1881 for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1882 br->GetEntry(itrack);
1883 Int_t nref=tcaRef->GetEntriesFast();
1884 if(nref>0) { //it is enough to look at the first one
1885 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1886 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1889 if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1890 br->GetEntry(imatch); // redundant, nevertheless ...
1891 Int_t nref=tcaRef->GetEntriesFast();
1892 for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1893 tref=(AliTrackReference*)tcaRef->At(iref);
1894 if(tref->R()>10) continue; // not SPD ref
1895 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1896 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1898 // compute the proper quantities for this tref, as was done for fClustersLay1/2
1899 Float_t x = tref->X() - vtx[0];
1900 Float_t y = tref->Y() - vtx[1];
1901 Float_t z = tref->Z() - vtx[2];
1903 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1905 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1906 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1907 trefLayExtr[2] = z; // Store z
1909 if(layer==1) { // try to see if it is reconstructable at the outer layer
1910 // find the difference in angles
1911 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1912 // take into account boundary condition
1913 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1915 // find the difference in z (between linear projection from layer 1
1916 // and the actual point: Dzeta= z1/r1*r2 -z2)
1917 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1918 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1920 // make "elliptical" cut in Phi and Zeta!
1921 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
1922 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
1923 if (d<1) {ret=kTRUE; break;}
1925 if(layer==0) { // try to see if it is reconstructable at the inner layer
1927 // find the difference in angles
1928 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1929 // take into account boundary condition
1930 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1932 // find the difference in z (between linear projection from layer 2
1933 // and the actual point: Dzeta= z2/r2*r1 -z1)
1934 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1935 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1937 // make "elliptical" cut in Phi and Zeta!
1938 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1939 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1940 if (d<1) {ret=kTRUE; break;};
1946 //_________________________________________________________________________
1947 void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
1949 // this method apply a rotation by 180 degree around the Z (beam) axis to all
1950 // the RecPoints in a given layer to be used to build tracklets.
1951 // **************** VERY IMPORTANT:: ***************
1952 // It must be called just after LoadClusterArrays, since afterwards the datamember
1953 // fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
1954 // instead of Cartesian
1956 if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1957 AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1959 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1960 fClustersLay1[iC1][0]*=-1;
1961 fClustersLay1[iC1][1]*=-1;
1965 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1966 fClustersLay2[iC2][0]*=-1;
1967 fClustersLay2[iC2][1]*=-1;
1972 //____________________________________________________________________________
1973 Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *){
1974 // This method is used to find the tracklets.
1975 // It is called from AliReconstruction
1976 // The vertex is supposed to be associated to the Tracker (i.e. to this) already
1977 // The cluster is supposed to be associated to the Tracker already
1978 // In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
1981 AliRunLoader* runLoader = AliRunLoader::Instance();
1983 Error("Clusters2Tracks", "no run loader found");
1986 AliStack *pStack=0x0; TTree *tRefTree=0x0;
1988 runLoader->LoadKinematics("read");
1989 runLoader->LoadTrackRefs("read");
1990 pStack= runLoader->Stack();
1991 tRefTree= runLoader->TreeTR();
1993 Reconstruct(pStack,tRefTree);
1995 if (GetLightBkgStudyInParallel()) {
1996 AliStack *dummy1=0x0; TTree *dummy2=0x0;
1997 ReflectClusterAroundZAxisForLayer(1);
1998 Reconstruct(dummy1,dummy2,kTRUE);
2002 //____________________________________________________________________________
2003 Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
2005 // It is called from AliReconstruction
2011 if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
2012 if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
2013 if(GetLightBkgStudyInParallel()) {
2014 TString name="AliITSPlaneEffSPDtrackletBkg.root";
2015 TFile* pefile = TFile::Open(name, "RECREATE");
2016 rc*=fPlaneEffBkg->Write();
2021 //____________________________________________________________________
2023 AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
2025 // - gets the clusters from the cluster tree
2026 // - convert them into global coordinates
2027 // - store them in the internal arrays
2028 // - count the number of cluster-fired chips
2030 //AliDebug(1,"Loading clusters and cluster-fired chips ...");
2035 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
2036 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
2038 itsClusterBranch->SetAddress(&itsClusters);
2040 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
2041 Float_t cluGlo[3]={0.,0.,0.};
2043 // loop over the its subdetectors
2044 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
2046 if (!itsClusterTree->GetEvent(iIts))
2049 Int_t nClusters = itsClusters->GetEntriesFast();
2051 // number of clusters in each chip of the current module
2054 // loop over clusters
2055 while(nClusters--) {
2056 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
2058 layer = cluster->GetLayer();
2059 if (layer>1) continue;
2061 cluster->GetGlobalXYZ(cluGlo);
2062 Float_t x = cluGlo[0];
2063 Float_t y = cluGlo[1];
2064 Float_t z = cluGlo[2];
2067 fClustersLay1[fNClustersLay1][0] = x;
2068 fClustersLay1[fNClustersLay1][1] = y;
2069 fClustersLay1[fNClustersLay1][2] = z;
2071 for (Int_t i=0; i<3; i++)
2072 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
2075 Int_t det=cluster->GetDetectorIndex();
2076 if(det<0 || det>79) {AliError("Cluster with det. index out of boundaries"); return;}
2077 fhClustersInModuleLay1[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2081 fClustersLay2[fNClustersLay2][0] = x;
2082 fClustersLay2[fNClustersLay2][1] = y;
2083 fClustersLay2[fNClustersLay2][2] = z;
2085 for (Int_t i=0; i<3; i++)
2086 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
2089 Int_t det=cluster->GetDetectorIndex();
2090 if(det<0 || det>159) {AliError("Cluster with det. index out of boundaries"); return;}
2091 fhClustersInModuleLay2[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2095 }// end of cluster loop
2097 } // end of its "subdetector" loop
2099 itsClusters->Delete();
2103 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
2105 //_________________________________________________________________________
2107 AliITSTrackleterSPDEff::SetLightBkgStudyInParallel(Bool_t b) {
2109 // - set Bool_t fLightBackgroundStudyInParallel = b
2110 // a) if you set this kTRUE, then the estimation of the
2111 // SPD efficiency is done as usual for data, but in
2112 // parallel a light (i.e. without control histograms, etc.)
2113 // evaluation of combinatorial background is performed
2114 // with the usual ReflectClusterAroundZAxisForLayer method.
2115 // b) if you set this kFALSE, then you would not have a second
2116 // container for PlaneEfficiency statistics to be used for background
2117 // (fPlaneEffBkg=0). If you want to have a full evaluation of the
2118 // background (with all control histograms and additional data
2119 // members referring to the background) then you have to call the
2120 // method SetReflectClusterAroundZAxisForLayer(kTRUE) esplicitily
2121 fLightBkgStudyInParallel=b;
2122 if(fLightBkgStudyInParallel) {
2123 if(!fPlaneEffBkg) fPlaneEffBkg = new AliITSPlaneEffSPD();
2126 delete fPlaneEffBkg;