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),
96 fReflectClusterAroundZAxisForLayer0(kFALSE),
97 fReflectClusterAroundZAxisForLayer1(kFALSE),
99 fUseOnlyPrimaryForPred(1),
100 fUseOnlySecondaryForPred(0),
101 fUseOnlySameParticle(0),
102 fUseOnlyDifferentParticle(0),
103 fUseOnlyStableParticle(0),
104 fPredictionPrimary(0),
105 fPredictionSecondary(0),
107 fClusterSecondary(0),
116 fhClustersDPhiInterpAcc(0),
117 fhClustersDThetaInterpAcc(0),
118 fhClustersDZetaInterpAcc(0),
119 fhClustersDPhiInterpAll(0),
120 fhClustersDThetaInterpAll(0),
121 fhClustersDZetaInterpAll(0),
122 fhDPhiVsDThetaInterpAll(0),
123 fhDPhiVsDThetaInterpAcc(0),
124 fhDPhiVsDZetaInterpAll(0),
125 fhDPhiVsDZetaInterpAcc(0),
126 fhetaClustersLay2(0),
127 fhphiClustersLay2(0),
129 fhClustersInModuleLay1(0),
130 fhClustersInModuleLay2(0)
132 // default constructor
133 // from AliITSMultReconstructor
136 SetOnlyOneTrackletPerC2();
137 fClustersLay1 = new Float_t*[300000];
138 fClustersLay2 = new Float_t*[300000];
139 fTracklets = new Float_t*[300000];
140 fAssociationFlag = new Bool_t[300000];
144 SetOnlyOneTrackletPerC1();
146 fAssociationFlag1 = new Bool_t[300000];
147 fChipPredOnLay2 = new UInt_t[300000];
148 fChipPredOnLay1 = new UInt_t[300000];
149 fChipUpdatedInEvent = new Bool_t[1200];
151 for(Int_t i=0; i<300000; i++) {
152 // from AliITSMultReconstructor
153 fClustersLay1[i] = new Float_t[6];
154 fClustersLay2[i] = new Float_t[6];
155 fTracklets[i] = new Float_t[5];
156 fAssociationFlag[i] = kFALSE;
158 fAssociationFlag1[i] = kFALSE;
160 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
162 if (GetHistOn()) BookHistos();
164 fPlaneEffSPD = new AliITSPlaneEffSPD();
166 //______________________________________________________________________
167 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
169 // from AliITSMultReconstructor
170 fClustersLay1(mr.fClustersLay1),
171 fClustersLay2(mr.fClustersLay2),
172 fTracklets(mr.fTracklets),
173 fAssociationFlag(mr.fAssociationFlag),
174 fNClustersLay1(mr.fNClustersLay1),
175 fNClustersLay2(mr.fNClustersLay2),
176 fNTracklets(mr.fNTracklets),
177 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
178 fPhiWindowL2(mr.fPhiWindowL2),
179 fZetaWindowL2(mr.fZetaWindowL2),
180 fPhiOverlapCut(mr.fPhiOverlapCut),
181 fZetaOverlapCut(mr.fZetaOverlapCut),
183 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
184 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
185 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
186 fhClustersDPhiAll(mr.fhClustersDPhiAll),
187 fhClustersDThetaAll(mr.fhClustersDThetaAll),
188 fhClustersDZetaAll(mr.fhClustersDZetaAll),
189 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
190 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
191 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
192 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
193 fhetaTracklets(mr.fhetaTracklets),
194 fhphiTracklets(mr.fhphiTracklets),
195 fhetaClustersLay1(mr.fhetaClustersLay1),
196 fhphiClustersLay1(mr.fhphiClustersLay1),
198 fAssociationFlag1(mr.fAssociationFlag1),
199 fChipPredOnLay2(mr.fChipPredOnLay2),
200 fChipPredOnLay1(mr.fChipPredOnLay1),
201 fNTracklets1(mr.fNTracklets1),
202 fPhiWindowL1(mr.fPhiWindowL1),
203 fZetaWindowL1(mr.fZetaWindowL1),
204 fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
205 fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
206 fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
207 fPlaneEffSPD(mr.fPlaneEffSPD),
208 fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
209 fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
211 fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
212 fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
213 fUseOnlySameParticle(mr.fUseOnlySameParticle),
214 fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
215 fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
216 fPredictionPrimary(mr.fPredictionPrimary),
217 fPredictionSecondary(mr.fPredictionSecondary),
218 fClusterPrimary(mr.fClusterPrimary),
219 fClusterSecondary(mr.fClusterSecondary),
220 fSuccessPP(mr.fSuccessPP),
221 fSuccessTT(mr.fSuccessTT),
222 fSuccessS(mr.fSuccessS),
223 fSuccessP(mr.fSuccessP),
224 fFailureS(mr.fFailureS),
225 fFailureP(mr.fFailureP),
227 fNonRecons(mr.fNonRecons),
228 fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
229 fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
230 fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
231 fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
232 fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
233 fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
234 fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
235 fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
236 fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
237 fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
238 fhetaClustersLay2(mr.fhetaClustersLay2),
239 fhphiClustersLay2(mr.fhphiClustersLay2),
240 fhClustersInChip(mr.fhClustersInChip),
241 fhClustersInModuleLay1(mr.fhClustersInModuleLay1),
242 fhClustersInModuleLay2(mr.fhClustersInModuleLay2)
247 //______________________________________________________________________
248 AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
249 // Assignment operator
250 this->~AliITSTrackleterSPDEff();
251 new(this) AliITSTrackleterSPDEff(mr);
254 //______________________________________________________________________
255 AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
257 // from AliITSMultReconstructor
259 for(Int_t i=0; i<300000; i++) {
260 delete [] fClustersLay1[i];
261 delete [] fClustersLay2[i];
262 delete [] fTracklets[i];
264 delete [] fClustersLay1;
265 delete [] fClustersLay2;
266 delete [] fTracklets;
267 delete [] fAssociationFlag;
272 delete [] fAssociationFlag1;
274 delete [] fChipPredOnLay2;
275 delete [] fChipPredOnLay1;
277 delete [] fChipUpdatedInEvent;
279 delete [] fPredictionPrimary;
280 delete [] fPredictionSecondary;
281 delete [] fClusterPrimary;
282 delete [] fClusterSecondary;
283 delete [] fSuccessPP;
284 delete [] fSuccessTT;
290 delete [] fNonRecons;
295 //____________________________________________________________________
297 //AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t*, AliStack *pStack, TTree *tRef) {
298 AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef) {
300 // - you have to take care of the following, before of using Reconstruct
301 // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
302 // and convert the cluster coordinates to theta, phi (seen from the
303 // interaction vertex).
304 // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
305 // - Find the extrapolation/interpolation point.
306 // - Find the chip corresponding to that
307 // - Check if there is a cluster near that point
311 // retrieve the vertex position
313 vtx[0]=(Float_t)GetX();
314 vtx[1]=(Float_t)GetY();
315 vtx[2]=(Float_t)GetZ();
316 // to study residual background (i.e. contribution from TT' to measured efficiency)
317 if(fReflectClusterAroundZAxisForLayer0) ReflectClusterAroundZAxisForLayer(0);
318 if(fReflectClusterAroundZAxisForLayer1) ReflectClusterAroundZAxisForLayer(1);
320 if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
321 if(fMC && !tRef) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
323 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
324 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
325 Int_t nfClu1=0; Int_t ntClu1=0;
326 Int_t nfClu2=0; Int_t ntClu2=0;
328 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
329 // for this new event)
330 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
332 // find the tracklets
333 AliDebug(1,"Looking for tracklets... ");
334 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
336 //###########################################################
337 // Loop on layer 1 : finding theta, phi and z
339 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
340 Float_t x = fClustersLay1[iC1][0] - vtx[0];
341 Float_t y = fClustersLay1[iC1][1] - vtx[1];
342 Float_t z = fClustersLay1[iC1][2] - vtx[2];
344 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
346 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
347 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
348 fClustersLay1[iC1][2] = z; // Store z
350 // find the Radius and the chip corresponding to the extrapolation point
352 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
354 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
355 key=999999; // also some other actions should be taken if not Found
357 nfTraPred2+=(Int_t)found; // this for debugging purpose
358 ntTraPred2++; // to check efficiency of the method FindChip
359 fChipPredOnLay2[iC1] = key;
360 fAssociationFlag1[iC1] = kFALSE;
363 Float_t eta=fClustersLay1[iC1][0];
364 eta= TMath::Tan(eta/2.);
365 eta=-TMath::Log(eta);
366 fhetaClustersLay1->Fill(eta);
367 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
368 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
371 // Loop on layer 2 : finding theta, phi and r
372 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
373 Float_t x = fClustersLay2[iC2][0] - vtx[0];
374 Float_t y = fClustersLay2[iC2][1] - vtx[1];
375 Float_t z = fClustersLay2[iC2][2] - vtx[2];
377 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
379 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
380 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
381 fClustersLay2[iC2][2] = z; // Store z
383 // find the Radius and the chip corresponding to the extrapolation point
385 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
387 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
390 nfTraPred1+=(Int_t)found; // this for debugging purpose
391 ntTraPred1++; // to check efficiency of the method FindChip
392 fChipPredOnLay1[iC2] = key;
393 fAssociationFlag[iC2] = kFALSE;
396 Float_t eta=fClustersLay2[iC2][0];
397 eta= TMath::Tan(eta/2.);
398 eta=-TMath::Log(eta);
399 fhetaClustersLay2->Fill(eta);
400 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
401 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
405 //###########################################################
407 // First part : Extrapolation to Layer 2
410 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
412 // here the control to check whether the efficiency of the chip traversed by this tracklet
413 // prediction has already been updated in this event using another tracklet prediction
414 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
416 // reset of variables for multiple candidates
417 Int_t iC2WithBestDist = 0; // reset
418 Float_t distmin = 100.; // just to put a huge number!
419 Float_t dPhimin = 0.; // Used for histograms only!
420 Float_t dThetamin = 0.; // Used for histograms only!
421 Float_t dZetamin = 0.; // Used for histograms only!
423 // in any case, if MC has been required, store statistics of primaries and secondaries
424 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
426 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
427 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
428 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
429 // do it always as a function of the chip number used to built the prediction
430 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
431 if (!found) {AliWarning(
432 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
434 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
435 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
436 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
437 { // this cluster is from a primary particle
438 fClusterPrimary[key]++;
440 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
441 } else { // this cluster is from a secondary particle
442 fClusterSecondary[key]++;
444 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
447 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
448 // (in case the prediction is reliable)
449 if( fChipPredOnLay2[iC1]<1200) {
450 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
451 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
452 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
453 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
454 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
455 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
456 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
457 else fNonRecons[fChipPredOnLay2[iC1]]++;
462 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
464 // The following excludes double associations
465 if (!fAssociationFlag[iC2]) {
467 // find the difference in angles
468 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
469 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
470 // take into account boundary condition
471 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
473 // find the difference in z (between linear projection from layer 1
474 // and the actual point: Dzeta= z1/r1*r2 -z2)
475 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
476 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
479 fhClustersDPhiAll->Fill(dPhi);
480 fhClustersDThetaAll->Fill(dTheta);
481 fhClustersDZetaAll->Fill(dZeta);
482 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
483 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
486 // make "elliptical" cut in Phi and Zeta!
487 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
488 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
492 //look for the minimum distance: the minimum is in iC2WithBestDist
493 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
494 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
498 iC2WithBestDist = iC2;
501 } // end of loop over clusters in layer 2
503 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
506 fhClustersDPhiAcc->Fill(dPhimin);
507 fhClustersDThetaAcc->Fill(dThetamin);
508 fhClustersDZetaAcc->Fill(dZetamin);
509 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
510 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
513 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
514 // flag the association
516 // store the tracklet
518 // use the theta from the clusters in the first layer
519 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
520 // use the phi from the clusters in the first layer
521 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
522 // Store the difference between phi1 and phi2
523 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
530 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
542 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
546 fTracklets[fNTracklets][3] = -2;
550 Float_t eta=fTracklets[fNTracklets][0];
551 eta= TMath::Tan(eta/2.);
552 eta=-TMath::Log(eta);
553 fhetaTracklets->Fill(eta);
554 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
557 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
558 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
561 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
564 nfClu2+=(Int_t)found; // this for debugging purpose
565 ntClu2++; // to check efficiency of the method FindChip
566 if(key<1200) { // the Chip has been found
567 if(fMC) { // this part only for MC
568 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
569 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
570 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
573 if(primary) fSuccessPP[key]++;
575 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
576 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
579 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
581 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
582 fChipUpdatedInEvent[key]=kTRUE;
584 if(primary) fSuccessP[key]++;
585 if(secondary) fSuccessS[key]++;
589 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
590 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
592 if(primary) fSuccessP[key]++;
593 if(secondary) fSuccessS[key]++;
600 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
601 else if (fChipPredOnLay2[iC1]<1200) {
602 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
603 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
605 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
606 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
609 } // end of loop over clusters in layer 1
611 fNTracklets1=fNTracklets;
613 //###################################################################
615 // Second part : Interpolation to Layer 1
618 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
620 // here the control to check whether the efficiency of the chip traversed by this tracklet
621 // prediction has already been updated in this event using another tracklet prediction
622 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
624 // reset of variables for multiple candidates
625 Int_t iC1WithBestDist = 0; // reset
626 Float_t distmin = 100.; // just to put a huge number!
627 Float_t dPhimin = 0.; // Used for histograms only!
628 Float_t dThetamin = 0.; // Used for histograms only!
629 Float_t dZetamin = 0.; // Used for histograms only!
631 // in any case, if MC has been required, store statistics of primaries and secondaries
632 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
634 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
635 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
636 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
637 // do it always as a function of the chip number used to built the prediction
638 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
639 if (!found) {AliWarning(
640 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
642 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
643 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
644 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
645 { // this cluster is from a primary particle
646 fClusterPrimary[key]++;
648 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
649 } else { // this cluster is from a secondary particle
650 fClusterSecondary[key]++;
652 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
655 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
656 // (in case the prediction is reliable)
657 if( fChipPredOnLay1[iC2]<1200) {
658 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
659 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
660 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
661 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
662 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
663 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
664 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
665 else fNonRecons[fChipPredOnLay1[iC2]]++;
670 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
672 // The following excludes double associations
673 if (!fAssociationFlag1[iC1]) {
675 // find the difference in angles
676 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
677 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
678 // take into account boundary condition
679 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
682 // find the difference in z (between linear projection from layer 2
683 // and the actual point: Dzeta= z2/r2*r1 -z1)
684 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
685 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
689 fhClustersDPhiInterpAll->Fill(dPhi);
690 fhClustersDThetaInterpAll->Fill(dTheta);
691 fhClustersDZetaInterpAll->Fill(dZeta);
692 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
693 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
695 // make "elliptical" cut in Phi and Zeta!
696 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
697 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
701 //look for the minimum distance: the minimum is in iC1WithBestDist
702 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
703 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
707 iC1WithBestDist = iC1;
710 } // end of loop over clusters in layer 1
712 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
715 fhClustersDPhiInterpAcc->Fill(dPhimin);
716 fhClustersDThetaInterpAcc->Fill(dThetamin);
717 fhClustersDZetaInterpAcc->Fill(dZetamin);
718 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
719 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
722 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
723 // flag the association
725 // store the tracklet
727 // use the theta from the clusters in the first layer
728 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
729 // use the phi from the clusters in the first layer
730 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
731 // Store the difference between phi1 and phi2
732 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
739 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
751 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
755 fTracklets[fNTracklets][3] = -2;
758 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
759 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
762 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
765 nfClu1+=(Int_t)found; // this for debugging purpose
766 ntClu1++; // to check efficiency of the method FindChip
768 if(fMC) { // this part only for MC
769 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
770 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
771 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
772 if (label2 < 3) { // same label
774 if(primary) fSuccessPP[key]++;
776 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
777 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
780 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
782 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
783 fChipUpdatedInEvent[key]=kTRUE;
785 if(primary) fSuccessP[key]++;
786 if(secondary) fSuccessS[key]++;
789 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
790 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
792 if(primary) fSuccessP[key]++;
793 if(secondary) fSuccessS[key]++;
800 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
801 else if (fChipPredOnLay1[iC2]<1200) {
802 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
803 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
805 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
806 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
809 } // end of loop over clusters in layer 2
811 AliDebug(1,Form("%d tracklets found", fNTracklets));
812 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
813 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
814 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
815 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
817 //____________________________________________________________________
818 Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
819 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
821 // Input: a) layer number in the range [0,1]
822 // b) vtx[3]: actual vertex
823 // c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
824 // d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
826 // Output: Unique key to locate a chip
827 // return: kTRUE if succesfull
829 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
830 Double_t r=GetRLayer(layer);
831 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
833 // set phiVtx in the range [0,2pi]
834 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
836 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
837 // of the intersection of the tracklet with the pixel layer.
838 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
839 else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
840 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
841 Double_t vtxy[2]={vtx[0],vtx[1]};
842 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
843 // this method gives you two interceptions
844 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
845 // set phiAbs in the range [0,2pi]
846 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
847 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
848 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
849 // taking the closest one to phiVtx
850 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
851 } else phiAbs=phiVtx;
852 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
854 // now you need to locate the chip within the idet detector,
855 // starting from the local coordinates in such a detector
857 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
858 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
859 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
861 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
864 //______________________________________________________________________________
865 Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
867 // Return the average radius of a layer from Geometry
869 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
870 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
872 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
873 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
874 Double_t r=TMath::Sqrt(x*x + y*y);
876 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
877 r += TMath::Sqrt(x*x + y*y);
878 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
879 r += TMath::Sqrt(x*x + y*y);
880 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
881 r += TMath::Sqrt(x*x + y*y);
885 //______________________________________________________________________________
886 Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
887 Float_t &xloc, Float_t &zloc) {
888 // this method transform Global Cilindrical coordinates into local (i.e. module)
889 // cartesian coordinates
891 //Compute Cartesian Global Coordinate
892 Double_t xyzGlob[3],xyzLoc[3];
894 xyzGlob[0]=r*TMath::Cos(phi);
895 xyzGlob[1]=r*TMath::Sin(phi);
900 if(idet<0) return kFALSE;
902 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
903 Int_t lad = Int_t(idet/ndet) + 1;
904 Int_t det = idet - (lad-1)*ndet + 1;
906 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
908 xloc = (Float_t)xyzLoc[0];
909 zloc = (Float_t)xyzLoc[2];
913 //______________________________________________________________________________
914 Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
915 //--------------------------------------------------------------------
916 // This function finds the detector crossed by the track
917 // Input: layer in range [0,1]
918 // phi in ALICE absolute reference system
920 //--------------------------------------------------------------------
921 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
922 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
923 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
924 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
926 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
927 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
928 Double_t phiOffset=TMath::ATan2(y,x);
932 if (zOffset<0) // old geometry
933 dphi = -(phi-phiOffset);
935 dphi = phi-phiOffset;
937 if (dphi < 0) dphi += 2*TMath::Pi();
938 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
939 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
940 if (np>=nladders) np-=nladders;
941 if (np<0) np+=nladders;
943 Double_t dz=zOffset-z;
944 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
945 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
946 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
947 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
949 return np*ndetectors + nz;
951 //____________________________________________________________
952 Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
953 // this method find the intersection in xy between a tracklet (straight line) and
954 // a circonference (r=R), using polar coordinates.
956 Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
957 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
958 - R: radius of the circle
959 Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
961 Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
962 r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
963 In the same system, the equation of a semi-line is: phi=phiVtx;
964 Hence you get one interception only: P=(r,phiVtx)
965 Finally you want P in the ABSOLUTE ALICE reference system.
967 Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
968 Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
969 Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
970 Double_t cC=rO*rO-R*R;
971 Double_t dDelta=bB*bB-4*cC;
972 if(dDelta<0) return kFALSE;
974 r1=(-bB-TMath::Sqrt(dDelta))/2;
975 r2=(-bB+TMath::Sqrt(dDelta))/2;
976 if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
977 Double_t r=TMath::Max(r1,r2); // take the positive
978 Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
979 Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
980 pvtx[0]=r*TMath::Cos(phiVtx);
981 pvtx[1]=r*TMath::Sin(phiVtx);
982 pP[0]=vtx[0]+pvtx[0];
983 pP[1]=vtx[1]+pvtx[1];
984 phi=TMath::ATan2(pP[1],pP[0]);
987 //___________________________________________________________
988 Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
990 // simple method to reduce all angles (in rad)
994 while(angle >=2*TMath::Pi() || angle<0) {
995 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
996 if(angle < 0) angle+=2*TMath::Pi();
1000 //___________________________________________________________
1001 Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
1003 // This method check if a particle is primary; i.e.
1004 // it comes from the main vertex and it is a "stable" particle, according to
1005 // AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
1006 // a stable particle: it has no effect on this analysis).
1007 // This method can be called only for MC events, where Kinematics is available.
1008 // if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
1009 // returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
1010 // The latter (see below) try to verify if a primary particle is also "detectable".
1012 if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
1013 if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
1014 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
1015 // return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
1016 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
1017 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
1018 TParticle* part = stack->Particle(ipart);
1019 TParticle* part0 = stack->Particle(0); // first primary
1020 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
1021 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
1022 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
1024 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
1025 TParticlePDG* pdgPart = part->GetPDG();
1026 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
1028 Double_t distx = part->Vx() - part0->Vx();
1029 Double_t disty = part->Vy() - part0->Vy();
1030 Double_t distz = part->Vz() - part0->Vz();
1031 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
1033 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
1034 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
1035 return kFALSE; }// primary if within 500 microns from true Vertex
1037 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
1040 //_____________________________________________________________________________________________
1041 Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
1043 // This private method can be applied on MC particles (if stack is available),
1044 // provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
1046 // It define "detectable" a primary particle according to the following criteria:
1048 // - if no decay products can be found in the stack (note that this does not
1049 // means it is stable, since a particle is stored in stack if it has at least 1 hit in a
1050 // sensitive detector)
1051 // - if it has at least one decay daughter produced outside or just on the outer pixel layer
1052 // - if the last decay particle is an electron (or a muon) which is not produced in-between
1053 // the two pixel layers (this is likely to be a kink).
1054 if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
1055 if(!stack) {AliError("null pointer to MC stack"); return 0;}
1056 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
1058 TParticle* part = stack->Particle(ipart);
1059 //TParticle* part0 = stack->Particle(0); // first primary
1065 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
1066 // its real fate ! But you have to take it !
1067 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
1068 Int_t lastDau = part->GetLastDaughter();
1069 nDau = lastDau - firstDau + 1;
1070 Double_t distMax=0.;
1072 for(Int_t j=firstDau; j<=lastDau; j++) {
1073 dau = stack->Particle(j);
1074 Double_t distx = dau->Vx();
1075 Double_t disty = dau->Vy();
1076 //Double_t distz = dau->Vz();
1077 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
1078 if(distR<distMax) continue; // considere only the daughter produced at largest radius
1082 dau = stack->Particle(jmax);
1083 pdgDau=dau->GetPdgCode();
1084 if (pdgDau == 11 || pdgDau == 13 ) {
1085 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
1086 else nret =0; // delta-ray emission in material (keep it)
1088 else {// not ele or muon
1089 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
1093 //_________________________________________________________________
1094 void AliITSTrackleterSPDEff::InitPredictionMC() {
1096 // this method allocate memory for the MC related informations
1097 // all the counters are set to 0
1100 if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1101 fPredictionPrimary = new Int_t[1200];
1102 fPredictionSecondary = new Int_t[1200];
1103 fClusterPrimary = new Int_t[1200];
1104 fClusterSecondary = new Int_t[1200];
1105 fSuccessPP = new Int_t[1200];
1106 fSuccessTT = new Int_t[1200];
1107 fSuccessS = new Int_t[1200];
1108 fSuccessP = new Int_t[1200];
1109 fFailureS = new Int_t[1200];
1110 fFailureP = new Int_t[1200];
1111 fRecons = new Int_t[1200];
1112 fNonRecons = new Int_t[1200];
1113 for(Int_t i=0; i<1200; i++) {
1114 fPredictionPrimary[i]=0;
1115 fPredictionSecondary[i]=0;
1116 fPredictionSecondary[i]=0;
1117 fClusterSecondary[i]=0;
1129 //_________________________________________________________________
1130 void AliITSTrackleterSPDEff::DeletePredictionMC() {
1132 // this method deallocate memory for the MC related informations
1133 // all the counters are set to 0
1136 if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1137 if(fPredictionPrimary) {
1138 delete fPredictionPrimary; fPredictionPrimary=0;
1140 if(fPredictionSecondary) {
1141 delete fPredictionSecondary; fPredictionSecondary=0;
1143 if(fClusterPrimary) {
1144 delete fClusterPrimary; fClusterPrimary=0;
1146 if(fClusterSecondary) {
1147 delete fClusterSecondary; fClusterSecondary=0;
1150 delete fSuccessPP; fSuccessPP=0;
1153 delete fSuccessTT; fSuccessTT=0;
1156 delete fSuccessS; fSuccessS=0;
1159 delete fSuccessP; fSuccessP=0;
1162 delete fFailureS; fFailureS=0;
1165 delete fFailureP; fFailureP=0;
1168 delete fRecons; fRecons=0;
1171 delete fNonRecons; fNonRecons=0;
1175 //______________________________________________________________________
1176 Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
1178 // This method return the Data menmber fPredictionPrimary [1200].
1179 // You can call it only for MC events.
1180 // fPredictionPrimary[key] contains the number of tracklet predictions on the
1181 // given chip key built using a cluster on the other layer produced (at least)
1182 // from a primary particle.
1183 // Key refers to the chip crossed by the prediction
1186 if (!fMC) {CallWarningMC(); return 0;}
1187 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1188 return fPredictionPrimary[(Int_t)key];
1190 //______________________________________________________________________
1191 Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
1193 // This method return the Data menmber fPredictionSecondary [1200].
1194 // You can call it only for MC events.
1195 // fPredictionSecondary[key] contains the number of tracklet predictions on the
1196 // given chip key built using a cluster on the other layer produced (only)
1197 // from a secondary particle
1198 // Key refers to the chip crossed by the prediction
1201 if (!fMC) {CallWarningMC(); return 0;}
1202 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1203 return fPredictionSecondary[(Int_t)key];
1205 //______________________________________________________________________
1206 Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
1208 // This method return the Data menmber fClusterPrimary [1200].
1209 // You can call it only for MC events.
1210 // fClusterPrimary[key] contains the number of tracklet predictions
1211 // built using a cluster on that layer produced (only)
1212 // from a primary particle
1213 // Key refers to the chip used to build the prediction
1216 if (!fMC) {CallWarningMC(); return 0;}
1217 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1218 return fClusterPrimary[(Int_t)key];
1220 //______________________________________________________________________
1221 Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
1223 // This method return the Data menmber fClusterSecondary [1200].
1224 // You can call it only for MC events.
1225 // fClusterSecondary[key] contains the number of tracklet predictions
1226 // built using a cluster on that layer produced (only)
1227 // from a secondary particle
1228 // Key refers to the chip used to build the prediction
1230 if (!fMC) {CallWarningMC(); return 0;}
1231 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1232 return fClusterSecondary[(Int_t)key];
1234 //______________________________________________________________________
1235 Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1237 // This method return the Data menmber fSuccessPP [1200].
1238 // You can call it only for MC events.
1239 // fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1240 // with a cluster on the other layer) built by using the same primary particle
1241 // the unique chip key refers to the chip which get updated its efficiency
1243 if (!fMC) {CallWarningMC(); return 0;}
1244 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1245 return fSuccessPP[(Int_t)key];
1247 //______________________________________________________________________
1248 Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1250 // This method return the Data menmber fSuccessTT [1200].
1251 // You can call it only for MC events.
1252 // fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1253 // with a cluster on the other layer) built by using the same particle (whatever)
1254 // the unique chip key refers to the chip which get updated its efficiency
1256 if (!fMC) {CallWarningMC(); return 0;}
1257 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1258 return fSuccessTT[(Int_t)key];
1260 //______________________________________________________________________
1261 Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1263 // This method return the Data menmber fSuccessS [1200].
1264 // You can call it only for MC events.
1265 // fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1266 // with a cluster on the other layer) built by using a secondary particle
1267 // the unique chip key refers to the chip which get updated its efficiency
1269 if (!fMC) {CallWarningMC(); return 0;}
1270 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1271 return fSuccessS[(Int_t)key];
1273 //______________________________________________________________________
1274 Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1276 // This method return the Data menmber fSuccessP [1200].
1277 // You can call it only for MC events.
1278 // fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1279 // with a cluster on the other layer) built by using a primary particle
1280 // the unique chip key refers to the chip which get updated its efficiency
1282 if (!fMC) {CallWarningMC(); return 0;}
1283 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1284 return fSuccessP[(Int_t)key];
1286 //______________________________________________________________________
1287 Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1289 // This method return the Data menmber fFailureS [1200].
1290 // You can call it only for MC events.
1291 // fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1292 // with a cluster on the other layer) built by using a secondary particle
1293 // the unique chip key refers to the chip which get updated its efficiency
1295 if (!fMC) {CallWarningMC(); return 0;}
1296 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1297 return fFailureS[(Int_t)key];
1299 //______________________________________________________________________
1300 Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1302 // This method return the Data menmber fFailureP [1200].
1303 // You can call it only for MC events.
1304 // fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1305 // with a cluster on the other layer) built by using a primary particle
1306 // the unique chip key refers to the chip which get updated its efficiency
1308 if (!fMC) {CallWarningMC(); return 0;}
1309 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1310 return fFailureP[(Int_t)key];
1312 //_____________________________________________________________________
1313 Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1315 // This method return the Data menmber fRecons [1200].
1316 // You can call it only for MC events.
1317 // fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1318 // has an hit in the detector)
1319 // the unique chip key refers to the chip where fall the prediction
1321 if (!fMC) {CallWarningMC(); return 0;}
1322 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1323 return fRecons[(Int_t)key];
1325 //_____________________________________________________________________
1326 Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1328 // This method return the Data menmber fNonRecons [1200].
1329 // You can call it only for MC events.
1330 // fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1331 // has not any hit in the detector)
1332 // the unique chip key refers to the chip where fall the prediction
1334 if (!fMC) {CallWarningMC(); return 0;}
1335 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1336 return fNonRecons[(Int_t)key];
1338 //______________________________________________________________________
1339 void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1340 // Print out some class data values in Ascii Form to output stream
1342 // ostream *os Output stream where Ascii data is to be writen
1347 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindowL2 <<" "<< fZetaWindowL2
1348 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1349 << " " << fUpdateOncePerEventPlaneEff << " " << fReflectClusterAroundZAxisForLayer0
1350 << " " << fReflectClusterAroundZAxisForLayer1;
1352 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1353 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1354 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1355 << " " << fUseOnlyStableParticle ;
1356 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1357 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1358 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1359 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
1360 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1361 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1362 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1363 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1364 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1365 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1366 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1367 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
1370 //______________________________________________________________________
1371 void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1372 // Read in some class data values in Ascii Form to output stream
1374 // istream *is Input stream where Ascii data is to be read in from
1381 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindowL2 >> fZetaWindowL2
1382 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1383 >> fUpdateOncePerEventPlaneEff >> fReflectClusterAroundZAxisForLayer0
1384 >> fReflectClusterAroundZAxisForLayer1;
1385 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
1387 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1389 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1390 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1391 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1392 >> fUseOnlyStableParticle;
1393 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1394 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1395 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1396 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1397 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1398 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1399 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1400 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1401 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1402 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1403 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1404 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1408 //______________________________________________________________________
1409 ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1410 // Standard output streaming function
1412 // ostream &os output steam
1413 // AliITSTrackleterSPDEff &s class to be streamed.
1417 // ostream &os The stream pointer
1422 //______________________________________________________________________
1423 istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1424 // Standard inputput streaming function
1426 // istream &is input steam
1427 // AliITSTrackleterSPDEff &s class to be streamed.
1431 // ostream &os The stream pointer
1433 //printf("prova %d \n", (Int_t)s.GetMC());
1437 //______________________________________________________________________
1438 void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
1440 // This Method write into an either asci or root file
1441 // the used cuts and the statistics of the MC related quantities
1442 // The method SetMC() has to be called before
1443 // Input TString filename: name of file for output (it deletes already existing
1448 //if(!fMC) {CallWarningMC(); return;}
1449 if (!filename.Contains(".root")) {
1450 ofstream out(filename.Data(),ios::out | ios::binary);
1456 TFile* mcfile = TFile::Open(filename, "RECREATE");
1457 TH1F* cuts = new TH1F("cuts", "list of cuts", 10, 0, 10); // TH1I containing cuts
1458 cuts->SetBinContent(1,fPhiWindowL1);
1459 cuts->SetBinContent(2,fZetaWindowL1);
1460 cuts->SetBinContent(3,fPhiWindowL2);
1461 cuts->SetBinContent(4,fZetaWindowL2);
1462 cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
1463 cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
1464 cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
1465 cuts->SetBinContent(8,fReflectClusterAroundZAxisForLayer0);
1466 cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer1);
1467 cuts->SetBinContent(10,fMC);
1470 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
1472 TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
1473 mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
1474 mc0->SetBinContent(2,fUseOnlySecondaryForPred);
1475 mc0->SetBinContent(3,fUseOnlySameParticle);
1476 mc0->SetBinContent(4,fUseOnlyDifferentParticle);
1477 mc0->SetBinContent(5,fUseOnlyStableParticle);
1481 mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
1482 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
1484 mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
1485 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
1487 mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
1488 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
1490 mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
1491 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
1493 mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
1494 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
1496 mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
1497 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
1499 mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
1500 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
1502 mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
1503 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
1505 mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
1506 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
1508 mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
1509 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
1511 mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
1512 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
1514 mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
1515 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
1523 //____________________________________________________________________
1524 void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
1526 // This Method read from an asci file (do not know why binary does not work)
1527 // the cuts to be used and the statistics of the MC related quantities
1528 // Input TString filename: name of input file for output
1529 // The method SetMC() has to be called before
1533 //if(!fMC) {CallWarningMC(); return;}
1534 if( gSystem->AccessPathName( filename.Data() ) ) {
1535 AliError( Form( "file (%s) not found", filename.Data() ) );
1539 if (!filename.Contains(".root")) {
1540 ifstream in(filename.Data(),ios::in | ios::binary);
1547 TFile *mcfile = TFile::Open(filename);
1548 TH1F *cuts = (TH1F*)mcfile->Get("cuts");
1549 fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
1550 fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
1551 fPhiWindowL2=(Float_t)cuts->GetBinContent(3);
1552 fZetaWindowL2=(Float_t)cuts->GetBinContent(4);
1553 fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
1554 fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
1555 fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
1556 fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(8);
1557 fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(9);
1558 fMC=(Bool_t)cuts->GetBinContent(10);
1559 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1560 else { // only if file with MC predictions
1561 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1562 TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
1563 fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
1564 fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
1565 fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
1566 fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
1567 fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
1569 mc1 =(TH1I*)mcfile->Get("mc1");
1570 for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1571 mc1 =(TH1I*)mcfile->Get("mc2");
1572 for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1573 mc1 =(TH1I*)mcfile->Get("mc3");
1574 for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1575 mc1 =(TH1I*)mcfile->Get("mc4");
1576 for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1577 mc1 =(TH1I*)mcfile->Get("mc5");
1578 for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1579 mc1 =(TH1I*)mcfile->Get("mc6");
1580 for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
1581 mc1 =(TH1I*)mcfile->Get("mc7");
1582 for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1583 mc1 =(TH1I*)mcfile->Get("mc8");
1584 for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1585 mc1 =(TH1I*)mcfile->Get("mc9");
1586 for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1587 mc1 =(TH1I*)mcfile->Get("mc10");
1588 for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1589 mc1 =(TH1I*)mcfile->Get("mc11");
1590 for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1591 mc1 =(TH1I*)mcfile->Get("mc12");
1592 for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1598 //____________________________________________________________________
1599 Bool_t AliITSTrackleterSPDEff::SaveHists() {
1600 // This (private) method save the histograms on the output file
1601 // (only if fHistOn is TRUE).
1602 // Also the histograms from the base class are saved through the
1603 // AliITSMultReconstructor::SaveHists() call
1605 if (!GetHistOn()) return kFALSE;
1607 // AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1608 fhClustersDPhiAll->Write();
1609 fhClustersDThetaAll->Write();
1610 fhClustersDZetaAll->Write();
1611 fhDPhiVsDThetaAll->Write();
1612 fhDPhiVsDZetaAll->Write();
1614 fhClustersDPhiAcc->Write();
1615 fhClustersDThetaAcc->Write();
1616 fhClustersDZetaAcc->Write();
1617 fhDPhiVsDThetaAcc->Write();
1618 fhDPhiVsDZetaAcc->Write();
1620 fhetaTracklets->Write();
1621 fhphiTracklets->Write();
1622 fhetaClustersLay1->Write();
1623 fhphiClustersLay1->Write();
1625 fhClustersDPhiInterpAll->Write();
1626 fhClustersDThetaInterpAll->Write();
1627 fhClustersDZetaInterpAll->Write();
1628 fhDPhiVsDThetaInterpAll->Write();
1629 fhDPhiVsDZetaInterpAll->Write();
1631 fhClustersDPhiInterpAcc->Write();
1632 fhClustersDThetaInterpAcc->Write();
1633 fhClustersDZetaInterpAcc->Write();
1634 fhDPhiVsDThetaInterpAcc->Write();
1635 fhDPhiVsDZetaInterpAcc->Write();
1637 fhetaClustersLay2->Write();
1638 fhphiClustersLay2->Write();
1639 fhClustersInChip->Write();
1640 for (Int_t nhist=0;nhist<80;nhist++){
1641 fhClustersInModuleLay1[nhist]->Write();
1643 for (Int_t nhist=0;nhist<160;nhist++){
1644 fhClustersInModuleLay2[nhist]->Write();
1648 //__________________________________________________________
1649 Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1651 // Saves the histograms into a tree and saves the trees into a file
1652 // Also the histograms from the base class are saved
1654 if (!GetHistOn()) return kFALSE;
1655 if (!strcmp(filename.Data(),"")) {
1656 AliWarning("WriteHistosToFile: null output filename!");
1659 TFile *hFile=new TFile(filename.Data(),option,
1660 "The File containing the histos for SPD efficiency studies with tracklets");
1661 if(!SaveHists()) return kFALSE;
1666 //____________________________________________________________
1667 void AliITSTrackleterSPDEff::BookHistos() {
1669 // This method books addtitional histograms
1670 // w.r.t. those of the base class.
1671 // In particular, the differences of cluster coordinate between the two SPD
1672 // layers are computed in the interpolation phase
1674 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1676 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
1677 fhClustersDPhiAcc->SetDirectory(0);
1678 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
1679 fhClustersDThetaAcc->SetDirectory(0);
1680 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
1681 fhClustersDZetaAcc->SetDirectory(0);
1683 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
1684 fhDPhiVsDZetaAcc->SetDirectory(0);
1685 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
1686 fhDPhiVsDThetaAcc->SetDirectory(0);
1688 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
1689 fhClustersDPhiAll->SetDirectory(0);
1690 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
1691 fhClustersDThetaAll->SetDirectory(0);
1692 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
1693 fhClustersDZetaAll->SetDirectory(0);
1695 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
1696 fhDPhiVsDZetaAll->SetDirectory(0);
1697 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
1698 fhDPhiVsDThetaAll->SetDirectory(0);
1700 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
1701 fhetaTracklets->SetDirectory(0);
1702 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
1703 fhphiTracklets->SetDirectory(0);
1704 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
1705 fhetaClustersLay1->SetDirectory(0);
1706 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
1707 fhphiClustersLay1->SetDirectory(0);
1709 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1710 fhClustersDPhiInterpAcc->SetDirectory(0);
1711 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1712 fhClustersDThetaInterpAcc->SetDirectory(0);
1713 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1714 fhClustersDZetaInterpAcc->SetDirectory(0);
1716 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1717 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1718 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1719 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1721 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1722 fhClustersDPhiInterpAll->SetDirectory(0);
1723 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1724 fhClustersDThetaInterpAll->SetDirectory(0);
1725 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1726 fhClustersDZetaInterpAll->SetDirectory(0);
1728 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1729 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1730 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1731 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1733 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1734 fhetaClustersLay2->SetDirectory(0);
1735 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1736 fhphiClustersLay2->SetDirectory(0);
1737 fhClustersInChip = new TH1F("fhClustersInChip", "ClustersPerChip", 1200, -0.5, 1199.5);
1738 fhClustersInChip->SetDirectory(0);
1739 // each chip is divided 8(z) x 4(y), i.e. in 32 squares, each containing 4 columns and 64 rows.
1740 Float_t bz[160]; const Float_t kconv = 1.0E-04; // converts microns to cm.
1741 for(Int_t i=0;i<160;i++) bz[i] = 425.0; // most are 425 microns except below
1742 bz[ 31] = bz[ 32] = 625.0; // first chip boundry
1743 bz[ 63] = bz[ 64] = 625.0; // first chip boundry
1744 bz[ 95] = bz[ 96] = 625.0; // first chip boundry
1745 bz[127] = bz[128] = 625.0; // first chip boundry
1746 Double_t xbins[41]; // each bin in x (Z loc coordinate) includes 4 columns
1748 Float_t xmn,xmx,zmn,zmx;
1749 if(!fPlaneEffSPD->GetBlockBoundaries(0,xmn,xmx,zmn,zmx)) AliWarning("Could not book histo properly");
1750 xbins[0]=(Double_t)zmn;
1751 for(Int_t i=0;i<40;i++) {
1752 xbins[i+1]=xbins[i] + (bz[4*i]+bz[4*i+1]+bz[4*i+2]+bz[4*i+3])*kconv;
1754 TString histname="ClustersLay1_mod_",aux;
1755 fhClustersInModuleLay1 =new TH2F*[80];
1756 for (Int_t nhist=0;nhist<80;nhist++){
1760 fhClustersInModuleLay1[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1761 fhClustersInModuleLay1[nhist]->SetName(aux.Data());
1762 fhClustersInModuleLay1[nhist]->SetTitle(aux.Data());
1763 fhClustersInModuleLay1[nhist]->SetDirectory(0);
1765 histname="ClustersLay2_mod_";
1766 fhClustersInModuleLay2 =new TH2F*[160];
1767 for (Int_t nhist=0;nhist<160;nhist++){
1770 fhClustersInModuleLay2[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1771 fhClustersInModuleLay2[nhist]->SetName(aux.Data());
1772 fhClustersInModuleLay2[nhist]->SetTitle(aux.Data());
1773 fhClustersInModuleLay2[nhist]->SetDirectory(0);
1778 //____________________________________________________________
1779 void AliITSTrackleterSPDEff::DeleteHistos() {
1781 // Private method to delete Histograms from memory
1782 // it is called. e.g., by the destructor.
1784 // form AliITSMultReconstructor
1785 if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
1786 if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
1787 if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
1788 if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
1789 if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
1790 if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
1791 if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
1792 if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
1793 if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
1794 if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
1795 if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
1796 if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
1797 if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
1798 if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
1800 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1801 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1802 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1803 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1804 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1805 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1806 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1807 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1808 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1809 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1810 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1811 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1812 if(fhClustersInChip) {delete fhClustersInChip; fhClustersInChip=0;}
1813 if(fhClustersInModuleLay1) {
1814 for (Int_t i=0; i<80; i++ ) delete fhClustersInModuleLay1[i];
1815 delete [] fhClustersInModuleLay1; fhClustersInModuleLay1=0;
1817 if(fhClustersInModuleLay2) {
1818 for (Int_t i=0; i<160; i++ ) delete fhClustersInModuleLay2[i];
1819 delete [] fhClustersInModuleLay2; fhClustersInModuleLay2=0;
1822 //_______________________________________________________________
1823 Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1824 Float_t* vtx, AliStack *stack, TTree *ref) {
1825 // This (private) method can be used only for MC events, where both AliStack and the TrackReference
1827 // It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1829 // - Int_t layer (either 0 or 1): layer which you want to chech if the tracklete can be
1831 // - Int_t iC : cluster index used to build the tracklet prediction
1832 // if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1833 // - Float_t* vtx: actual event vertex
1834 // - stack: pointer to Stack
1835 // - ref: pointer to TTRee of TrackReference
1836 Bool_t ret=kFALSE; // returned value
1837 Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1838 if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1839 if(!stack) {AliError("null pointer to MC stack"); return ret;}
1840 if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1841 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1842 if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1844 AliTrackReference *tref=0x0;
1845 Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1846 Int_t nentries = (Int_t)ref->GetEntries();
1847 TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1848 TBranch *br = ref->GetBranch("TrackReferences");
1849 br->SetAddress(&tcaRef);
1850 for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1851 br->GetEntry(itrack);
1852 Int_t nref=tcaRef->GetEntriesFast();
1853 if(nref>0) { //it is enough to look at the first one
1854 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1855 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1858 if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1859 br->GetEntry(imatch); // redundant, nevertheless ...
1860 Int_t nref=tcaRef->GetEntriesFast();
1861 for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1862 tref=(AliTrackReference*)tcaRef->At(iref);
1863 if(tref->R()>10) continue; // not SPD ref
1864 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1865 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1867 // compute the proper quantities for this tref, as was done for fClustersLay1/2
1868 Float_t x = tref->X() - vtx[0];
1869 Float_t y = tref->Y() - vtx[1];
1870 Float_t z = tref->Z() - vtx[2];
1872 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1874 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1875 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1876 trefLayExtr[2] = z; // Store z
1878 if(layer==1) { // try to see if it is reconstructable at the outer layer
1879 // find the difference in angles
1880 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1881 // take into account boundary condition
1882 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1884 // find the difference in z (between linear projection from layer 1
1885 // and the actual point: Dzeta= z1/r1*r2 -z2)
1886 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1887 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1889 // make "elliptical" cut in Phi and Zeta!
1890 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
1891 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
1892 if (d<1) {ret=kTRUE; break;}
1894 if(layer==0) { // try to see if it is reconstructable at the inner layer
1896 // find the difference in angles
1897 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1898 // take into account boundary condition
1899 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1901 // find the difference in z (between linear projection from layer 2
1902 // and the actual point: Dzeta= z2/r2*r1 -z1)
1903 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1904 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1906 // make "elliptical" cut in Phi and Zeta!
1907 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1908 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1909 if (d<1) {ret=kTRUE; break;};
1915 //_________________________________________________________________________
1916 void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
1918 // this method apply a rotation by 180 degree around the Z (beam) axis to all
1919 // the RecPoints in a given layer to be used to build tracklets.
1920 // **************** VERY IMPORTANT:: ***************
1921 // It must be called just after LoadClusterArrays, since afterwards the datamember
1922 // fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
1923 // instead of Cartesian
1925 if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1926 AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1928 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1929 fClustersLay1[iC1][0]*=-1;
1930 fClustersLay1[iC1][1]*=-1;
1934 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1935 fClustersLay2[iC2][0]*=-1;
1936 fClustersLay2[iC2][1]*=-1;
1941 //____________________________________________________________________________
1942 Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *){
1943 // This method is used to find the tracklets.
1944 // It is called from AliReconstruction
1945 // The vertex is supposed to be associated to the Tracker (i.e. to this) already
1946 // The cluster is supposed to be associated to the Tracker already
1947 // In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
1950 AliRunLoader* runLoader = AliRunLoader::Instance();
1952 Error("Clusters2Tracks", "no run loader found");
1955 AliStack *pStack=0x0; TTree *tRefTree=0x0;
1957 runLoader->LoadKinematics("read");
1958 runLoader->LoadTrackRefs("read");
1959 pStack= runLoader->Stack();
1960 tRefTree= runLoader->TreeTR();
1962 Reconstruct(pStack,tRefTree);
1965 //____________________________________________________________________________
1966 Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
1968 // It is called from AliReconstruction
1974 if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
1975 if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
1978 //____________________________________________________________________
1980 AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
1982 // - gets the clusters from the cluster tree
1983 // - convert them into global coordinates
1984 // - store them in the internal arrays
1985 // - count the number of cluster-fired chips
1987 //AliDebug(1,"Loading clusters and cluster-fired chips ...");
1992 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
1993 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
1995 itsClusterBranch->SetAddress(&itsClusters);
1997 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
1998 Float_t cluGlo[3]={0.,0.,0.};
2000 // loop over the its subdetectors
2001 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
2003 if (!itsClusterTree->GetEvent(iIts))
2006 Int_t nClusters = itsClusters->GetEntriesFast();
2008 // number of clusters in each chip of the current module
2011 // loop over clusters
2012 while(nClusters--) {
2013 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
2015 layer = cluster->GetLayer();
2016 if (layer>1) continue;
2018 cluster->GetGlobalXYZ(cluGlo);
2019 Float_t x = cluGlo[0];
2020 Float_t y = cluGlo[1];
2021 Float_t z = cluGlo[2];
2024 fClustersLay1[fNClustersLay1][0] = x;
2025 fClustersLay1[fNClustersLay1][1] = y;
2026 fClustersLay1[fNClustersLay1][2] = z;
2028 for (Int_t i=0; i<3; i++)
2029 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
2032 Int_t det=cluster->GetDetectorIndex();
2033 if(det<0 || det>79) {AliError("Cluster with det. index out of boundaries"); return;}
2034 fhClustersInModuleLay1[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2038 fClustersLay2[fNClustersLay2][0] = x;
2039 fClustersLay2[fNClustersLay2][1] = y;
2040 fClustersLay2[fNClustersLay2][2] = z;
2042 for (Int_t i=0; i<3; i++)
2043 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
2046 Int_t det=cluster->GetDetectorIndex();
2047 if(det<0 || det>159) {AliError("Cluster with det. index out of boundaries"); return;}
2048 fhClustersInModuleLay2[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2052 }// end of cluster loop
2054 } // end of its "subdetector" loop
2056 itsClusters->Delete();
2060 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));