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 #include "AliESDEvent.h"
51 #include "AliESDVertex.h"
52 //____________________________________________________________________
53 ClassImp(AliITSTrackleterSPDEff)
56 //____________________________________________________________________
57 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
67 fOnlyOneTrackletPerC2(0),
74 fhClustersDThetaAcc(0),
75 fhClustersDZetaAcc(0),
77 fhClustersDThetaAll(0),
78 fhClustersDZetaAll(0),
94 fOnlyOneTrackletPerC1(0),
95 fUpdateOncePerEventPlaneEff(0),
97 fChipUpdatedInEvent(0),
100 fReflectClusterAroundZAxisForLayer0(kFALSE),
101 fReflectClusterAroundZAxisForLayer1(kFALSE),
102 fLightBkgStudyInParallel(kFALSE),
104 fUseOnlyPrimaryForPred(0),
105 fUseOnlySecondaryForPred(0),
106 fUseOnlySameParticle(0),
107 fUseOnlyDifferentParticle(0),
108 fUseOnlyStableParticle(1),
109 fPredictionPrimary(0),
110 fPredictionSecondary(0),
112 fClusterSecondary(0),
121 fhClustersDPhiInterpAcc(0),
122 fhClustersDThetaInterpAcc(0),
123 fhClustersDZetaInterpAcc(0),
124 fhClustersDPhiInterpAll(0),
125 fhClustersDThetaInterpAll(0),
126 fhClustersDZetaInterpAll(0),
127 fhDPhiVsDThetaInterpAll(0),
128 fhDPhiVsDThetaInterpAcc(0),
129 fhDPhiVsDZetaInterpAll(0),
130 fhDPhiVsDZetaInterpAcc(0),
131 fhetaClustersLay2(0),
132 fhphiClustersLay2(0),
134 fhClustersInModuleLay1(0),
135 fhClustersInModuleLay2(0)
137 // default constructor
138 // from AliITSMultReconstructor
141 //______________________________________________________________________
142 void AliITSTrackleterSPDEff::Init() {
145 SetOnlyOneTrackletPerC2();
146 fClustersLay1 = new Float_t*[300000];
147 fClustersLay2 = new Float_t*[300000];
148 fTracklets = new Float_t*[300000];
149 fAssociationFlag = new Bool_t[300000];
153 SetOnlyOneTrackletPerC1();
155 fAssociationFlag1 = new Bool_t[300000];
156 fChipPredOnLay2 = new UInt_t[300000];
157 fChipPredOnLay1 = new UInt_t[300000];
158 fChipUpdatedInEvent = new Bool_t[1200];
160 for(Int_t i=0; i<300000; i++) {
161 // from AliITSMultReconstructor
162 fClustersLay1[i] = new Float_t[6];
163 fClustersLay2[i] = new Float_t[6];
164 fTracklets[i] = new Float_t[5];
165 fAssociationFlag[i] = kFALSE;
167 fAssociationFlag1[i] = kFALSE;
169 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
171 if (GetHistOn()) BookHistos();
173 fPlaneEffSPD = new AliITSPlaneEffSPD();
174 SetLightBkgStudyInParallel();
176 //______________________________________________________________________
177 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
179 // from AliITSMultReconstructor
180 fClustersLay1(mr.fClustersLay1),
181 fClustersLay2(mr.fClustersLay2),
182 fTracklets(mr.fTracklets),
183 fAssociationFlag(mr.fAssociationFlag),
184 fNClustersLay1(mr.fNClustersLay1),
185 fNClustersLay2(mr.fNClustersLay2),
186 fNTracklets(mr.fNTracklets),
187 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
188 fPhiWindowL2(mr.fPhiWindowL2),
189 fZetaWindowL2(mr.fZetaWindowL2),
190 fPhiOverlapCut(mr.fPhiOverlapCut),
191 fZetaOverlapCut(mr.fZetaOverlapCut),
193 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
194 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
195 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
196 fhClustersDPhiAll(mr.fhClustersDPhiAll),
197 fhClustersDThetaAll(mr.fhClustersDThetaAll),
198 fhClustersDZetaAll(mr.fhClustersDZetaAll),
199 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
200 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
201 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
202 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
203 fhetaTracklets(mr.fhetaTracklets),
204 fhphiTracklets(mr.fhphiTracklets),
205 fhetaClustersLay1(mr.fhetaClustersLay1),
206 fhphiClustersLay1(mr.fhphiClustersLay1),
208 fAssociationFlag1(mr.fAssociationFlag1),
209 fChipPredOnLay2(mr.fChipPredOnLay2),
210 fChipPredOnLay1(mr.fChipPredOnLay1),
211 fNTracklets1(mr.fNTracklets1),
212 fPhiWindowL1(mr.fPhiWindowL1),
213 fZetaWindowL1(mr.fZetaWindowL1),
214 fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
215 fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
216 fMinContVtx(mr.fMinContVtx),
217 fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
218 fPlaneEffSPD(mr.fPlaneEffSPD),
219 fPlaneEffBkg(mr.fPlaneEffBkg),
220 fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
221 fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
222 fLightBkgStudyInParallel(mr.fLightBkgStudyInParallel),
224 fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
225 fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
226 fUseOnlySameParticle(mr.fUseOnlySameParticle),
227 fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
228 fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
229 fPredictionPrimary(mr.fPredictionPrimary),
230 fPredictionSecondary(mr.fPredictionSecondary),
231 fClusterPrimary(mr.fClusterPrimary),
232 fClusterSecondary(mr.fClusterSecondary),
233 fSuccessPP(mr.fSuccessPP),
234 fSuccessTT(mr.fSuccessTT),
235 fSuccessS(mr.fSuccessS),
236 fSuccessP(mr.fSuccessP),
237 fFailureS(mr.fFailureS),
238 fFailureP(mr.fFailureP),
240 fNonRecons(mr.fNonRecons),
241 fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
242 fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
243 fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
244 fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
245 fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
246 fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
247 fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
248 fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
249 fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
250 fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
251 fhetaClustersLay2(mr.fhetaClustersLay2),
252 fhphiClustersLay2(mr.fhphiClustersLay2),
253 fhClustersInChip(mr.fhClustersInChip),
254 fhClustersInModuleLay1(mr.fhClustersInModuleLay1),
255 fhClustersInModuleLay2(mr.fhClustersInModuleLay2)
260 //______________________________________________________________________
261 AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
262 // Assignment operator
263 this->~AliITSTrackleterSPDEff();
264 new(this) AliITSTrackleterSPDEff(mr);
267 //______________________________________________________________________
268 AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
270 // from AliITSMultReconstructor
272 for(Int_t i=0; i<300000; i++) {
273 delete [] fClustersLay1[i];
274 delete [] fClustersLay2[i];
275 delete [] fTracklets[i];
277 delete [] fClustersLay1;
278 delete [] fClustersLay2;
279 delete [] fTracklets;
280 delete [] fAssociationFlag;
285 delete [] fAssociationFlag1;
287 delete [] fChipPredOnLay2;
288 delete [] fChipPredOnLay1;
290 delete [] fChipUpdatedInEvent;
292 delete [] fPredictionPrimary;
293 delete [] fPredictionSecondary;
294 delete [] fClusterPrimary;
295 delete [] fClusterSecondary;
296 delete [] fSuccessPP;
297 delete [] fSuccessTT;
303 delete [] fNonRecons;
314 //____________________________________________________________________
316 AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef, Bool_t lbkg) {
318 // - you have to take care of the following, before of using Reconstruct
319 // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
320 // and convert the cluster coordinates to theta, phi (seen from the
321 // interaction vertex).
322 // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
323 // - Find the extrapolation/interpolation point.
324 // - Find the chip corresponding to that
325 // - Check if there is a cluster near that point
328 if(lbkg && !GetLightBkgStudyInParallel()) {
329 AliError("You asked for lightBackground in the Reconstruction without proper call to SetLightBkgStudyInParallel(1)");
332 AliITSPlaneEffSPD *pe;
339 // retrieve the vertex position
341 vtx[0]=(Float_t)GetX();
342 vtx[1]=(Float_t)GetY();
343 vtx[2]=(Float_t)GetZ();
344 // to study residual background (i.e. contribution from TT' to measured efficiency)
345 if(fReflectClusterAroundZAxisForLayer0 && !lbkg) ReflectClusterAroundZAxisForLayer(0);
346 if(fReflectClusterAroundZAxisForLayer1 && !lbkg) ReflectClusterAroundZAxisForLayer(1);
348 if(fMC && !pStack && !lbkg) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
349 if(fMC && !tRef && !lbkg) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
351 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
352 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
353 Int_t nfClu1=0; Int_t ntClu1=0;
354 Int_t nfClu2=0; Int_t ntClu2=0;
356 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
357 // for this new event)
358 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
360 // find the tracklets
361 AliDebug(1,"Looking for tracklets... ");
362 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
364 //###########################################################
365 // Loop on layer 1 : finding theta, phi and z
367 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
368 Float_t x = fClustersLay1[iC1][0] - vtx[0];
369 Float_t y = fClustersLay1[iC1][1] - vtx[1];
370 Float_t z = fClustersLay1[iC1][2] - vtx[2];
372 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
374 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
375 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
376 fClustersLay1[iC1][2] = z; // Store z
378 // find the Radius and the chip corresponding to the extrapolation point
380 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
382 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
383 key=999999; // also some other actions should be taken if not Found
385 nfTraPred2+=(Int_t)found; // this for debugging purpose
386 ntTraPred2++; // to check efficiency of the method FindChip
387 fChipPredOnLay2[iC1] = key;
388 fAssociationFlag1[iC1] = kFALSE;
390 if (fHistOn && !lbkg) {
391 Float_t eta=fClustersLay1[iC1][0];
392 eta= TMath::Tan(eta/2.);
393 eta=-TMath::Log(eta);
394 fhetaClustersLay1->Fill(eta);
395 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
396 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
399 // Loop on layer 2 : finding theta, phi and r
400 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
401 Float_t x = fClustersLay2[iC2][0] - vtx[0];
402 Float_t y = fClustersLay2[iC2][1] - vtx[1];
403 Float_t z = fClustersLay2[iC2][2] - vtx[2];
405 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
407 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
408 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
409 fClustersLay2[iC2][2] = z; // Store z
411 // find the Radius and the chip corresponding to the extrapolation point
413 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
415 AliDebug(1,Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
418 nfTraPred1+=(Int_t)found; // this for debugging purpose
419 ntTraPred1++; // to check efficiency of the method FindChip
420 fChipPredOnLay1[iC2] = key;
421 fAssociationFlag[iC2] = kFALSE;
423 if (fHistOn && !lbkg) {
424 Float_t eta=fClustersLay2[iC2][0];
425 eta= TMath::Tan(eta/2.);
426 eta=-TMath::Log(eta);
427 fhetaClustersLay2->Fill(eta);
428 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
429 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
433 //###########################################################
435 // First part : Extrapolation to Layer 2
438 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
440 // here the control to check whether the efficiency of the chip traversed by this tracklet
441 // prediction has already been updated in this event using another tracklet prediction
442 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
444 // reset of variables for multiple candidates
445 Int_t iC2WithBestDist = 0; // reset
446 Float_t distmin = 100.; // just to put a huge number!
447 Float_t dPhimin = 0.; // Used for histograms only!
448 Float_t dThetamin = 0.; // Used for histograms only!
449 Float_t dZetamin = 0.; // Used for histograms only!
451 // in any case, if MC has been required, store statistics of primaries and secondaries
452 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
454 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
455 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
456 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
457 // do it always as a function of the chip number used to built the prediction
458 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
459 if (!found) {AliDebug(1,
460 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
462 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
463 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
464 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
465 { // this cluster is from a primary particle
466 fClusterPrimary[key]++;
468 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
469 } else { // this cluster is from a secondary particle
470 fClusterSecondary[key]++;
472 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
475 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
476 // (in case the prediction is reliable)
477 if( fChipPredOnLay2[iC1]<1200) {
478 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
479 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
480 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
481 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
482 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
483 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
484 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
485 else fNonRecons[fChipPredOnLay2[iC1]]++;
490 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
492 // The following excludes double associations
493 if (!fAssociationFlag[iC2]) {
495 // find the difference in angles
496 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
497 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
498 // take into account boundary condition
499 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
501 // find the difference in z (between linear projection from layer 1
502 // and the actual point: Dzeta= z1/r1*r2 -z2)
503 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
504 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
506 if (fHistOn && !lbkg) {
507 fhClustersDPhiAll->Fill(dPhi);
508 fhClustersDThetaAll->Fill(dTheta);
509 fhClustersDZetaAll->Fill(dZeta);
510 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
511 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
514 // make "elliptical" cut in Phi and Zeta!
515 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
516 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
520 //look for the minimum distance: the minimum is in iC2WithBestDist
521 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
522 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
526 iC2WithBestDist = iC2;
529 } // end of loop over clusters in layer 2
531 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
533 if (fHistOn && !lbkg) {
534 fhClustersDPhiAcc->Fill(dPhimin);
535 fhClustersDThetaAcc->Fill(dThetamin);
536 fhClustersDZetaAcc->Fill(dZetamin);
537 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
538 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
541 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
542 // flag the association
544 // store the tracklet
546 // use the theta from the clusters in the first layer
547 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
548 // use the phi from the clusters in the first layer
549 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
550 // Store the difference between phi1 and phi2
551 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
558 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
570 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
574 fTracklets[fNTracklets][3] = -2;
577 if (fHistOn && !lbkg) {
578 Float_t eta=fTracklets[fNTracklets][0];
579 eta= TMath::Tan(eta/2.);
580 eta=-TMath::Log(eta);
581 fhetaTracklets->Fill(eta);
582 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
585 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
586 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
589 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
592 nfClu2+=(Int_t)found; // this for debugging purpose
593 ntClu2++; // to check efficiency of the method FindChip
594 if(key<1200) { // the Chip has been found
595 if(fMC && !lbkg) { // this part only for MC
596 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
597 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
598 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
601 if(primary) fSuccessPP[key]++;
603 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
604 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
607 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
609 pe->UpDatePlaneEff(kTRUE,key); // success
610 fChipUpdatedInEvent[key]=kTRUE;
612 if(primary) fSuccessP[key]++;
613 if(secondary) fSuccessS[key]++;
617 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
618 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
620 if(primary) fSuccessP[key]++;
621 if(secondary) fSuccessS[key]++;
628 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
629 else if (fChipPredOnLay2[iC1]<1200) {
630 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
631 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
633 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
634 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
637 } // end of loop over clusters in layer 1
639 fNTracklets1=fNTracklets;
641 //###################################################################
643 // Second part : Interpolation to Layer 1
646 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
648 // here the control to check whether the efficiency of the chip traversed by this tracklet
649 // prediction has already been updated in this event using another tracklet prediction
650 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
652 // reset of variables for multiple candidates
653 Int_t iC1WithBestDist = 0; // reset
654 Float_t distmin = 100.; // just to put a huge number!
655 Float_t dPhimin = 0.; // Used for histograms only!
656 Float_t dThetamin = 0.; // Used for histograms only!
657 Float_t dZetamin = 0.; // Used for histograms only!
659 // in any case, if MC has been required, store statistics of primaries and secondaries
660 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
662 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
663 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
664 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
665 // do it always as a function of the chip number used to built the prediction
666 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
667 if (!found) {AliDebug(1,
668 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
670 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
671 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
672 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
673 { // this cluster is from a primary particle
674 fClusterPrimary[key]++;
676 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
677 } else { // this cluster is from a secondary particle
678 fClusterSecondary[key]++;
680 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
683 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
684 // (in case the prediction is reliable)
685 if( fChipPredOnLay1[iC2]<1200) {
686 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
687 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
688 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
689 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
690 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
691 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
692 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
693 else fNonRecons[fChipPredOnLay1[iC2]]++;
698 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
700 // The following excludes double associations
701 if (!fAssociationFlag1[iC1]) {
703 // find the difference in angles
704 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
705 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
706 // take into account boundary condition
707 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
710 // find the difference in z (between linear projection from layer 2
711 // and the actual point: Dzeta= z2/r2*r1 -z1)
712 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
713 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
716 if (fHistOn && !lbkg) {
717 fhClustersDPhiInterpAll->Fill(dPhi);
718 fhClustersDThetaInterpAll->Fill(dTheta);
719 fhClustersDZetaInterpAll->Fill(dZeta);
720 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
721 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
723 // make "elliptical" cut in Phi and Zeta!
724 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
725 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
729 //look for the minimum distance: the minimum is in iC1WithBestDist
730 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
731 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
735 iC1WithBestDist = iC1;
738 } // end of loop over clusters in layer 1
740 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
742 if (fHistOn && !lbkg) {
743 fhClustersDPhiInterpAcc->Fill(dPhimin);
744 fhClustersDThetaInterpAcc->Fill(dThetamin);
745 fhClustersDZetaInterpAcc->Fill(dZetamin);
746 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
747 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
750 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
751 // flag the association
753 // store the tracklet
755 // use the theta from the clusters in the first layer
756 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
757 // use the phi from the clusters in the first layer
758 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
759 // Store the difference between phi1 and phi2
760 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
767 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
779 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
783 fTracklets[fNTracklets][3] = -2;
786 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
787 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
790 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
793 nfClu1+=(Int_t)found; // this for debugging purpose
794 ntClu1++; // to check efficiency of the method FindChip
796 if(fMC && !lbkg) { // this part only for MC
797 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
798 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
799 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
800 if (label2 < 3) { // same label
802 if(primary) fSuccessPP[key]++;
804 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
805 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
808 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
810 pe->UpDatePlaneEff(kTRUE,key); // success
811 fChipUpdatedInEvent[key]=kTRUE;
813 if(primary) fSuccessP[key]++;
814 if(secondary) fSuccessS[key]++;
817 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
818 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
820 if(primary) fSuccessP[key]++;
821 if(secondary) fSuccessS[key]++;
828 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
829 else if (fChipPredOnLay1[iC2]<1200) {
830 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
831 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
833 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
834 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
837 } // end of loop over clusters in layer 2
839 AliDebug(1,Form("%d tracklets found", fNTracklets));
840 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
841 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
842 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
843 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
845 //____________________________________________________________________
846 Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer,const Float_t* vtx,
847 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
849 // Input: a) layer number in the range [0,1]
850 // b) vtx[3]: actual vertex
851 // c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
852 // d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
854 // Output: Unique key to locate a chip
855 // return: kTRUE if succesfull
857 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
858 Double_t r=GetRLayer(layer);
859 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
861 // set phiVtx in the range [0,2pi]
862 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
864 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
865 // of the intersection of the tracklet with the pixel layer.
866 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
868 if(TMath::Abs(thetaVtx)<1E-6) return kFALSE;
869 zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
871 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
872 Double_t vtxy[2]={vtx[0],vtx[1]};
873 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
874 // this method gives you two interceptions
875 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
876 // set phiAbs in the range [0,2pi]
877 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
878 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
879 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
880 // taking the closest one to phiVtx
881 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
882 } else phiAbs=phiVtx;
883 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
885 // now you need to locate the chip within the idet detector,
886 // starting from the local coordinates in such a detector
888 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
889 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
890 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
892 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
895 //______________________________________________________________________________
896 Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
898 // Return the average radius of a layer from Geometry
900 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
901 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
903 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
904 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
905 Double_t r=TMath::Sqrt(x*x + y*y);
907 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
908 r += TMath::Sqrt(x*x + y*y);
909 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
910 r += TMath::Sqrt(x*x + y*y);
911 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
912 r += TMath::Sqrt(x*x + y*y);
916 //______________________________________________________________________________
917 Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
918 Float_t &xloc, Float_t &zloc) {
919 // this method transform Global Cilindrical coordinates into local (i.e. module)
920 // cartesian coordinates
922 //Compute Cartesian Global Coordinate
923 Double_t xyzGlob[3],xyzLoc[3];
925 xyzGlob[0]=r*TMath::Cos(phi);
926 xyzGlob[1]=r*TMath::Sin(phi);
931 if(idet<0) return kFALSE;
933 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
934 Int_t lad = Int_t(idet/ndet) + 1;
935 Int_t det = idet - (lad-1)*ndet + 1;
937 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
939 xloc = (Float_t)xyzLoc[0];
940 zloc = (Float_t)xyzLoc[2];
944 //______________________________________________________________________________
945 Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
946 //--------------------------------------------------------------------
947 // This function finds the detector crossed by the track
948 // Input: layer in range [0,1]
949 // phi in ALICE absolute reference system
951 //--------------------------------------------------------------------
952 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
953 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
954 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
955 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
957 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
958 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
959 Double_t phiOffset=TMath::ATan2(y,x);
963 if (zOffset<0) // old geometry
964 dphi = -(phi-phiOffset);
966 dphi = phi-phiOffset;
968 if (dphi < 0) dphi += 2*TMath::Pi();
969 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
970 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
971 if (np>=nladders) np-=nladders;
972 if (np<0) np+=nladders;
974 Double_t dz=zOffset-z;
975 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
976 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
977 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
978 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
980 return np*ndetectors + nz;
982 //____________________________________________________________
983 Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
984 // this method find the intersection in xy between a tracklet (straight line) and
985 // a circonference (r=R), using polar coordinates.
987 Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
988 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
989 - R: radius of the circle
990 Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
992 Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
993 r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
994 In the same system, the equation of a semi-line is: phi=phiVtx;
995 Hence you get one interception only: P=(r,phiVtx)
996 Finally you want P in the ABSOLUTE ALICE reference system.
998 Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
999 Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
1000 Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
1001 Double_t cC=rO*rO-R*R;
1002 Double_t dDelta=bB*bB-4*cC;
1003 if(dDelta<0) return kFALSE;
1005 r1=(-bB-TMath::Sqrt(dDelta))/2;
1006 r2=(-bB+TMath::Sqrt(dDelta))/2;
1007 if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
1008 Double_t r=TMath::Max(r1,r2); // take the positive
1009 Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
1010 Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
1011 pvtx[0]=r*TMath::Cos(phiVtx);
1012 pvtx[1]=r*TMath::Sin(phiVtx);
1013 pP[0]=vtx[0]+pvtx[0];
1014 pP[1]=vtx[1]+pvtx[1];
1015 phi=TMath::ATan2(pP[1],pP[0]);
1018 //___________________________________________________________
1019 Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) const {
1021 // simple method to reduce all angles (in rad)
1025 while(angle >=2*TMath::Pi() || angle<0) {
1026 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
1027 if(angle < 0) angle+=2*TMath::Pi();
1031 //___________________________________________________________
1032 Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
1034 // This method check if a particle is primary; i.e.
1035 // it comes from the main vertex and it is a "stable" particle, according to
1036 // AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
1037 // a stable particle: it has no effect on this analysis).
1038 // This method can be called only for MC events, where Kinematics is available.
1039 // if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
1040 // returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
1041 // The latter (see below) try to verify if a primary particle is also "detectable".
1043 if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
1044 if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
1045 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
1046 // return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
1047 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
1048 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
1049 TParticle* part = stack->Particle(ipart);
1050 TParticle* part0 = stack->Particle(0); // first primary
1051 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
1052 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
1053 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
1055 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
1056 TParticlePDG* pdgPart = part->GetPDG();
1057 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
1059 Double_t distx = part->Vx() - part0->Vx();
1060 Double_t disty = part->Vy() - part0->Vy();
1061 Double_t distz = part->Vz() - part0->Vz();
1062 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
1064 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
1065 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
1066 return kFALSE; }// primary if within 500 microns from true Vertex
1068 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
1071 //_____________________________________________________________________________________________
1072 Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
1074 // This private method can be applied on MC particles (if stack is available),
1075 // provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
1077 // It define "detectable" a primary particle according to the following criteria:
1079 // - if no decay products can be found in the stack (note that this does not
1080 // means it is stable, since a particle is stored in stack if it has at least 1 hit in a
1081 // sensitive detector)
1082 // - if it has at least one decay daughter produced outside or just on the outer pixel layer
1083 // - if the last decay particle is an electron (or a muon) which is not produced in-between
1084 // the two pixel layers (this is likely to be a kink).
1085 if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
1086 if(!stack) {AliError("null pointer to MC stack"); return 0;}
1087 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
1089 TParticle* part = stack->Particle(ipart);
1090 //TParticle* part0 = stack->Particle(0); // first primary
1096 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
1097 // its real fate ! But you have to take it !
1098 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
1099 Int_t lastDau = part->GetLastDaughter();
1100 nDau = lastDau - firstDau + 1;
1101 Double_t distMax=0.;
1103 for(Int_t j=firstDau; j<=lastDau; j++) {
1104 dau = stack->Particle(j);
1105 Double_t distx = dau->Vx();
1106 Double_t disty = dau->Vy();
1107 //Double_t distz = dau->Vz();
1108 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
1109 if(distR<distMax) continue; // considere only the daughter produced at largest radius
1113 dau = stack->Particle(jmax);
1114 pdgDau=dau->GetPdgCode();
1115 if (pdgDau == 11 || pdgDau == 13 ) {
1116 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
1117 else nret =0; // delta-ray emission in material (keep it)
1119 else {// not ele or muon
1120 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
1124 //_________________________________________________________________
1125 void AliITSTrackleterSPDEff::InitPredictionMC() {
1127 // this method allocate memory for the MC related informations
1128 // all the counters are set to 0
1131 if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1132 fPredictionPrimary = new Int_t[1200];
1133 fPredictionSecondary = new Int_t[1200];
1134 fClusterPrimary = new Int_t[1200];
1135 fClusterSecondary = new Int_t[1200];
1136 fSuccessPP = new Int_t[1200];
1137 fSuccessTT = new Int_t[1200];
1138 fSuccessS = new Int_t[1200];
1139 fSuccessP = new Int_t[1200];
1140 fFailureS = new Int_t[1200];
1141 fFailureP = new Int_t[1200];
1142 fRecons = new Int_t[1200];
1143 fNonRecons = new Int_t[1200];
1144 for(Int_t i=0; i<1200; i++) {
1145 fPredictionPrimary[i]=0;
1146 fPredictionSecondary[i]=0;
1147 fPredictionSecondary[i]=0;
1148 fClusterSecondary[i]=0;
1160 //_________________________________________________________________
1161 void AliITSTrackleterSPDEff::DeletePredictionMC() {
1163 // this method deallocate memory for the MC related informations
1164 // all the counters are set to 0
1167 if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1168 if(fPredictionPrimary) {
1169 delete fPredictionPrimary; fPredictionPrimary=0;
1171 if(fPredictionSecondary) {
1172 delete fPredictionSecondary; fPredictionSecondary=0;
1174 if(fClusterPrimary) {
1175 delete fClusterPrimary; fClusterPrimary=0;
1177 if(fClusterSecondary) {
1178 delete fClusterSecondary; fClusterSecondary=0;
1181 delete fSuccessPP; fSuccessPP=0;
1184 delete fSuccessTT; fSuccessTT=0;
1187 delete fSuccessS; fSuccessS=0;
1190 delete fSuccessP; fSuccessP=0;
1193 delete fFailureS; fFailureS=0;
1196 delete fFailureP; fFailureP=0;
1199 delete fRecons; fRecons=0;
1202 delete fNonRecons; fNonRecons=0;
1206 //______________________________________________________________________
1207 Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
1209 // This method return the Data menmber fPredictionPrimary [1200].
1210 // You can call it only for MC events.
1211 // fPredictionPrimary[key] contains the number of tracklet predictions on the
1212 // given chip key built using a cluster on the other layer produced (at least)
1213 // from a primary particle.
1214 // Key refers to the chip crossed by the prediction
1217 if (!fMC) {CallWarningMC(); return 0;}
1218 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1219 return fPredictionPrimary[(Int_t)key];
1221 //______________________________________________________________________
1222 Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
1224 // This method return the Data menmber fPredictionSecondary [1200].
1225 // You can call it only for MC events.
1226 // fPredictionSecondary[key] contains the number of tracklet predictions on the
1227 // given chip key built using a cluster on the other layer produced (only)
1228 // from a secondary particle
1229 // Key refers to the chip crossed by the prediction
1232 if (!fMC) {CallWarningMC(); return 0;}
1233 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1234 return fPredictionSecondary[(Int_t)key];
1236 //______________________________________________________________________
1237 Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
1239 // This method return the Data menmber fClusterPrimary [1200].
1240 // You can call it only for MC events.
1241 // fClusterPrimary[key] contains the number of tracklet predictions
1242 // built using a cluster on that layer produced (only)
1243 // from a primary particle
1244 // Key refers to the chip used to build the prediction
1247 if (!fMC) {CallWarningMC(); return 0;}
1248 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1249 return fClusterPrimary[(Int_t)key];
1251 //______________________________________________________________________
1252 Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
1254 // This method return the Data menmber fClusterSecondary [1200].
1255 // You can call it only for MC events.
1256 // fClusterSecondary[key] contains the number of tracklet predictions
1257 // built using a cluster on that layer produced (only)
1258 // from a secondary particle
1259 // Key refers to the chip used to build the prediction
1261 if (!fMC) {CallWarningMC(); return 0;}
1262 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1263 return fClusterSecondary[(Int_t)key];
1265 //______________________________________________________________________
1266 Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1268 // This method return the Data menmber fSuccessPP [1200].
1269 // You can call it only for MC events.
1270 // fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1271 // with a cluster on the other layer) built by using the same primary particle
1272 // the unique chip key refers to the chip which get updated its efficiency
1274 if (!fMC) {CallWarningMC(); return 0;}
1275 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1276 return fSuccessPP[(Int_t)key];
1278 //______________________________________________________________________
1279 Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1281 // This method return the Data menmber fSuccessTT [1200].
1282 // You can call it only for MC events.
1283 // fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1284 // with a cluster on the other layer) built by using the same particle (whatever)
1285 // the unique chip key refers to the chip which get updated its efficiency
1287 if (!fMC) {CallWarningMC(); return 0;}
1288 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1289 return fSuccessTT[(Int_t)key];
1291 //______________________________________________________________________
1292 Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1294 // This method return the Data menmber fSuccessS [1200].
1295 // You can call it only for MC events.
1296 // fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1297 // with a cluster on the other layer) built by using a secondary particle
1298 // the unique chip key refers to the chip which get updated its efficiency
1300 if (!fMC) {CallWarningMC(); return 0;}
1301 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1302 return fSuccessS[(Int_t)key];
1304 //______________________________________________________________________
1305 Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1307 // This method return the Data menmber fSuccessP [1200].
1308 // You can call it only for MC events.
1309 // fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1310 // with a cluster on the other layer) built by using a primary particle
1311 // the unique chip key refers to the chip which get updated its efficiency
1313 if (!fMC) {CallWarningMC(); return 0;}
1314 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1315 return fSuccessP[(Int_t)key];
1317 //______________________________________________________________________
1318 Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1320 // This method return the Data menmber fFailureS [1200].
1321 // You can call it only for MC events.
1322 // fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1323 // with a cluster on the other layer) built by using a secondary particle
1324 // the unique chip key refers to the chip which get updated its efficiency
1326 if (!fMC) {CallWarningMC(); return 0;}
1327 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1328 return fFailureS[(Int_t)key];
1330 //______________________________________________________________________
1331 Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1333 // This method return the Data menmber fFailureP [1200].
1334 // You can call it only for MC events.
1335 // fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1336 // with a cluster on the other layer) built by using a primary particle
1337 // the unique chip key refers to the chip which get updated its efficiency
1339 if (!fMC) {CallWarningMC(); return 0;}
1340 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1341 return fFailureP[(Int_t)key];
1343 //_____________________________________________________________________
1344 Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1346 // This method return the Data menmber fRecons [1200].
1347 // You can call it only for MC events.
1348 // fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1349 // has an hit in the detector)
1350 // the unique chip key refers to the chip where fall the prediction
1352 if (!fMC) {CallWarningMC(); return 0;}
1353 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1354 return fRecons[(Int_t)key];
1356 //_____________________________________________________________________
1357 Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1359 // This method return the Data menmber fNonRecons [1200].
1360 // You can call it only for MC events.
1361 // fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1362 // has not any hit in the detector)
1363 // the unique chip key refers to the chip where fall the prediction
1365 if (!fMC) {CallWarningMC(); return 0;}
1366 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1367 return fNonRecons[(Int_t)key];
1369 //______________________________________________________________________
1370 void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1371 // Print out some class data values in Ascii Form to output stream
1373 // ostream *os Output stream where Ascii data is to be writen
1378 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindowL2 <<" "<< fZetaWindowL2
1379 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1380 << " " << fUpdateOncePerEventPlaneEff << " " << fMinContVtx
1381 << " " << fReflectClusterAroundZAxisForLayer0
1382 << " " << fReflectClusterAroundZAxisForLayer1;
1384 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1385 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1386 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1387 << " " << fUseOnlyStableParticle ;
1388 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1389 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1390 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1391 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
1392 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1393 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1394 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1395 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1396 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1397 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1398 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1399 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
1402 //______________________________________________________________________
1403 void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1404 // Read in some class data values in Ascii Form to output stream
1406 // istream *is Input stream where Ascii data is to be read in from
1413 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindowL2 >> fZetaWindowL2
1414 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1415 >> fUpdateOncePerEventPlaneEff >> fMinContVtx
1416 >> fReflectClusterAroundZAxisForLayer0
1417 >> fReflectClusterAroundZAxisForLayer1;
1418 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
1420 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1422 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1423 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1424 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1425 >> fUseOnlyStableParticle;
1426 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1427 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1428 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1429 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1430 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1431 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1432 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1433 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1434 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1435 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1436 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1437 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1441 //______________________________________________________________________
1442 ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1443 // Standard output streaming function
1445 // ostream &os output steam
1446 // AliITSTrackleterSPDEff &s class to be streamed.
1450 // ostream &os The stream pointer
1455 //______________________________________________________________________
1456 istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1457 // Standard inputput streaming function
1459 // istream &is input steam
1460 // AliITSTrackleterSPDEff &s class to be streamed.
1464 // ostream &os The stream pointer
1466 //printf("prova %d \n", (Int_t)s.GetMC());
1470 //______________________________________________________________________
1471 void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
1473 // This Method write into an either asci or root file
1474 // the used cuts and the statistics of the MC related quantities
1475 // The method SetMC() has to be called before
1476 // Input TString filename: name of file for output (it deletes already existing
1481 //if(!fMC) {CallWarningMC(); return;}
1482 if (!filename.Contains(".root")) {
1483 ofstream out(filename.Data(),ios::out | ios::binary);
1489 TFile* mcfile = TFile::Open(filename, "RECREATE");
1490 TH1F* cuts = new TH1F("cuts", "list of cuts", 11, 0, 11); // TH1I containing cuts
1491 cuts->SetBinContent(1,fPhiWindowL1);
1492 cuts->SetBinContent(2,fZetaWindowL1);
1493 cuts->SetBinContent(3,fPhiWindowL2);
1494 cuts->SetBinContent(4,fZetaWindowL2);
1495 cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
1496 cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
1497 cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
1498 cuts->SetBinContent(8,fMinContVtx);
1499 cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer0);
1500 cuts->SetBinContent(10,fReflectClusterAroundZAxisForLayer1);
1501 cuts->SetBinContent(11,fMC);
1504 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
1506 TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
1507 mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
1508 mc0->SetBinContent(2,fUseOnlySecondaryForPred);
1509 mc0->SetBinContent(3,fUseOnlySameParticle);
1510 mc0->SetBinContent(4,fUseOnlyDifferentParticle);
1511 mc0->SetBinContent(5,fUseOnlyStableParticle);
1515 mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
1516 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
1518 mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
1519 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
1521 mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
1522 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
1524 mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
1525 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
1527 mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
1528 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
1530 mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
1531 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
1533 mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
1534 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
1536 mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
1537 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
1539 mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
1540 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
1542 mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
1543 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
1545 mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
1546 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
1548 mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
1549 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
1557 //____________________________________________________________________
1558 void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
1560 // This Method read from an asci file (do not know why binary does not work)
1561 // the cuts to be used and the statistics of the MC related quantities
1562 // Input TString filename: name of input file for output
1563 // The method SetMC() has to be called before
1567 //if(!fMC) {CallWarningMC(); return;}
1568 if( gSystem->AccessPathName( filename.Data() ) ) {
1569 AliError( Form( "file (%s) not found", filename.Data() ) );
1573 if (!filename.Contains(".root")) {
1574 ifstream in(filename.Data(),ios::in | ios::binary);
1581 TFile *mcfile = TFile::Open(filename);
1582 TH1F *cuts = (TH1F*)mcfile->Get("cuts");
1583 fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
1584 fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
1585 fPhiWindowL2=(Float_t)cuts->GetBinContent(3);
1586 fZetaWindowL2=(Float_t)cuts->GetBinContent(4);
1587 fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
1588 fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
1589 fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
1590 fMinContVtx=(Int_t)cuts->GetBinContent(8);
1591 fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(9);
1592 fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(10);
1593 fMC=(Bool_t)cuts->GetBinContent(11);
1594 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1595 else { // only if file with MC predictions
1596 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1597 TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
1598 fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
1599 fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
1600 fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
1601 fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
1602 fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
1604 mc1 =(TH1I*)mcfile->Get("mc1");
1605 for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1606 mc1 =(TH1I*)mcfile->Get("mc2");
1607 for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1608 mc1 =(TH1I*)mcfile->Get("mc3");
1609 for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1610 mc1 =(TH1I*)mcfile->Get("mc4");
1611 for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1612 mc1 =(TH1I*)mcfile->Get("mc5");
1613 for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1614 mc1 =(TH1I*)mcfile->Get("mc6");
1615 for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
1616 mc1 =(TH1I*)mcfile->Get("mc7");
1617 for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1618 mc1 =(TH1I*)mcfile->Get("mc8");
1619 for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1620 mc1 =(TH1I*)mcfile->Get("mc9");
1621 for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1622 mc1 =(TH1I*)mcfile->Get("mc10");
1623 for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1624 mc1 =(TH1I*)mcfile->Get("mc11");
1625 for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1626 mc1 =(TH1I*)mcfile->Get("mc12");
1627 for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1633 //____________________________________________________________________
1634 Bool_t AliITSTrackleterSPDEff::SaveHists() {
1635 // This (private) method save the histograms on the output file
1636 // (only if fHistOn is TRUE).
1637 // Also the histograms from the base class are saved through the
1638 // AliITSMultReconstructor::SaveHists() call
1640 if (!GetHistOn()) return kFALSE;
1642 // AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1643 fhClustersDPhiAll->Write();
1644 fhClustersDThetaAll->Write();
1645 fhClustersDZetaAll->Write();
1646 fhDPhiVsDThetaAll->Write();
1647 fhDPhiVsDZetaAll->Write();
1649 fhClustersDPhiAcc->Write();
1650 fhClustersDThetaAcc->Write();
1651 fhClustersDZetaAcc->Write();
1652 fhDPhiVsDThetaAcc->Write();
1653 fhDPhiVsDZetaAcc->Write();
1655 fhetaTracklets->Write();
1656 fhphiTracklets->Write();
1657 fhetaClustersLay1->Write();
1658 fhphiClustersLay1->Write();
1660 fhClustersDPhiInterpAll->Write();
1661 fhClustersDThetaInterpAll->Write();
1662 fhClustersDZetaInterpAll->Write();
1663 fhDPhiVsDThetaInterpAll->Write();
1664 fhDPhiVsDZetaInterpAll->Write();
1666 fhClustersDPhiInterpAcc->Write();
1667 fhClustersDThetaInterpAcc->Write();
1668 fhClustersDZetaInterpAcc->Write();
1669 fhDPhiVsDThetaInterpAcc->Write();
1670 fhDPhiVsDZetaInterpAcc->Write();
1672 fhetaClustersLay2->Write();
1673 fhphiClustersLay2->Write();
1674 fhClustersInChip->Write();
1675 for (Int_t nhist=0;nhist<80;nhist++){
1676 fhClustersInModuleLay1[nhist]->Write();
1678 for (Int_t nhist=0;nhist<160;nhist++){
1679 fhClustersInModuleLay2[nhist]->Write();
1683 //__________________________________________________________
1684 Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1686 // Saves the histograms into a tree and saves the trees into a file
1687 // Also the histograms from the base class are saved
1689 if (!GetHistOn()) return kFALSE;
1690 if (!strcmp(filename.Data(),"")) {
1691 AliWarning("WriteHistosToFile: null output filename!");
1694 TFile *hFile=new TFile(filename.Data(),option,
1695 "The File containing the histos for SPD efficiency studies with tracklets");
1696 if(!SaveHists()) return kFALSE;
1701 //____________________________________________________________
1702 void AliITSTrackleterSPDEff::BookHistos() {
1704 // This method books addtitional histograms
1705 // w.r.t. those of the base class.
1706 // In particular, the differences of cluster coordinate between the two SPD
1707 // layers are computed in the interpolation phase
1709 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1711 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
1712 fhClustersDPhiAcc->SetDirectory(0);
1713 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
1714 fhClustersDThetaAcc->SetDirectory(0);
1715 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
1716 fhClustersDZetaAcc->SetDirectory(0);
1718 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
1719 fhDPhiVsDZetaAcc->SetDirectory(0);
1720 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
1721 fhDPhiVsDThetaAcc->SetDirectory(0);
1723 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
1724 fhClustersDPhiAll->SetDirectory(0);
1725 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
1726 fhClustersDThetaAll->SetDirectory(0);
1727 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
1728 fhClustersDZetaAll->SetDirectory(0);
1730 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
1731 fhDPhiVsDZetaAll->SetDirectory(0);
1732 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
1733 fhDPhiVsDThetaAll->SetDirectory(0);
1735 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
1736 fhetaTracklets->SetDirectory(0);
1737 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
1738 fhphiTracklets->SetDirectory(0);
1739 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
1740 fhetaClustersLay1->SetDirectory(0);
1741 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
1742 fhphiClustersLay1->SetDirectory(0);
1744 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1745 fhClustersDPhiInterpAcc->SetDirectory(0);
1746 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1747 fhClustersDThetaInterpAcc->SetDirectory(0);
1748 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1749 fhClustersDZetaInterpAcc->SetDirectory(0);
1751 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1752 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1753 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1754 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1756 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1757 fhClustersDPhiInterpAll->SetDirectory(0);
1758 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1759 fhClustersDThetaInterpAll->SetDirectory(0);
1760 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1761 fhClustersDZetaInterpAll->SetDirectory(0);
1763 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1764 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1765 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1766 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1768 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1769 fhetaClustersLay2->SetDirectory(0);
1770 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1771 fhphiClustersLay2->SetDirectory(0);
1772 fhClustersInChip = new TH1F("fhClustersInChip", "ClustersPerChip", 1200, -0.5, 1199.5);
1773 fhClustersInChip->SetDirectory(0);
1774 // each chip is divided 8(z) x 4(y), i.e. in 32 squares, each containing 4 columns and 64 rows.
1775 Float_t bz[160]; const Float_t kconv = 1.0E-04; // converts microns to cm.
1776 for(Int_t i=0;i<160;i++) bz[i] = 425.0; // most are 425 microns except below
1777 bz[ 31] = bz[ 32] = 625.0; // first chip boundry
1778 bz[ 63] = bz[ 64] = 625.0; // first chip boundry
1779 bz[ 95] = bz[ 96] = 625.0; // first chip boundry
1780 bz[127] = bz[128] = 625.0; // first chip boundry
1781 Double_t xbins[41]; // each bin in x (Z loc coordinate) includes 4 columns
1783 Float_t xmn,xmx,zmn,zmx;
1784 if(!fPlaneEffSPD->GetBlockBoundaries(0,xmn,xmx,zmn,zmx)) AliWarning("Could not book histo properly");
1785 xbins[0]=(Double_t)zmn;
1786 for(Int_t i=0;i<40;i++) {
1787 xbins[i+1]=xbins[i] + (bz[4*i]+bz[4*i+1]+bz[4*i+2]+bz[4*i+3])*kconv;
1789 TString histname="ClustersLay1_mod_",aux;
1790 fhClustersInModuleLay1 =new TH2F*[80];
1791 for (Int_t nhist=0;nhist<80;nhist++){
1795 fhClustersInModuleLay1[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1796 fhClustersInModuleLay1[nhist]->SetName(aux.Data());
1797 fhClustersInModuleLay1[nhist]->SetTitle(aux.Data());
1798 fhClustersInModuleLay1[nhist]->SetDirectory(0);
1800 histname="ClustersLay2_mod_";
1801 fhClustersInModuleLay2 =new TH2F*[160];
1802 for (Int_t nhist=0;nhist<160;nhist++){
1805 fhClustersInModuleLay2[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1806 fhClustersInModuleLay2[nhist]->SetName(aux.Data());
1807 fhClustersInModuleLay2[nhist]->SetTitle(aux.Data());
1808 fhClustersInModuleLay2[nhist]->SetDirectory(0);
1813 //____________________________________________________________
1814 void AliITSTrackleterSPDEff::DeleteHistos() {
1816 // Private method to delete Histograms from memory
1817 // it is called. e.g., by the destructor.
1819 // form AliITSMultReconstructor
1820 if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
1821 if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
1822 if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
1823 if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
1824 if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
1825 if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
1826 if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
1827 if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
1828 if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
1829 if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
1830 if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
1831 if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
1832 if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
1833 if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
1835 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1836 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1837 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1838 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1839 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1840 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1841 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1842 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1843 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1844 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1845 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1846 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1847 if(fhClustersInChip) {delete fhClustersInChip; fhClustersInChip=0;}
1848 if(fhClustersInModuleLay1) {
1849 for (Int_t i=0; i<80; i++ ) delete fhClustersInModuleLay1[i];
1850 delete [] fhClustersInModuleLay1; fhClustersInModuleLay1=0;
1852 if(fhClustersInModuleLay2) {
1853 for (Int_t i=0; i<160; i++ ) delete fhClustersInModuleLay2[i];
1854 delete [] fhClustersInModuleLay2; fhClustersInModuleLay2=0;
1857 //_______________________________________________________________
1858 Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1859 const Float_t* vtx, const AliStack *stack, TTree *ref) {
1860 // This (private) method can be used only for MC events, where both AliStack and the TrackReference
1862 // It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1864 // - Int_t layer (either 0 or 1): layer which you want to check if the tracklete can be
1866 // - Int_t iC : cluster index used to build the tracklet prediction
1867 // if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1868 // - Float_t* vtx: actual event vertex
1869 // - stack: pointer to Stack
1870 // - ref: pointer to TTRee of TrackReference
1871 Bool_t ret=kFALSE; // returned value
1872 Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1873 if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1874 if(!stack) {AliError("null pointer to MC stack"); return ret;}
1875 if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1876 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1877 if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1879 AliTrackReference *tref=0x0;
1880 Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1881 Int_t nentries = (Int_t)ref->GetEntries();
1882 TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1883 TBranch *br = ref->GetBranch("TrackReferences");
1884 br->SetAddress(&tcaRef);
1885 for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1886 br->GetEntry(itrack);
1887 Int_t nref=tcaRef->GetEntriesFast();
1888 if(nref>0) { //it is enough to look at the first one
1889 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1890 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1893 if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1894 br->GetEntry(imatch); // redundant, nevertheless ...
1895 Int_t nref=tcaRef->GetEntriesFast();
1896 for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1897 tref=(AliTrackReference*)tcaRef->At(iref);
1898 if(tref->R()>10) continue; // not SPD ref
1899 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1900 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1902 // compute the proper quantities for this tref, as was done for fClustersLay1/2
1903 Float_t x = tref->X() - vtx[0];
1904 Float_t y = tref->Y() - vtx[1];
1905 Float_t z = tref->Z() - vtx[2];
1907 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1909 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1910 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1911 trefLayExtr[2] = z; // Store z
1913 if(layer==1) { // try to see if it is reconstructable at the outer layer
1914 // find the difference in angles
1915 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1916 // take into account boundary condition
1917 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1919 // find the difference in z (between linear projection from layer 1
1920 // and the actual point: Dzeta= z1/r1*r2 -z2)
1921 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1922 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1924 // make "elliptical" cut in Phi and Zeta!
1925 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
1926 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
1927 if (d<1) {ret=kTRUE; break;}
1929 if(layer==0) { // try to see if it is reconstructable at the inner layer
1931 // find the difference in angles
1932 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1933 // take into account boundary condition
1934 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1936 // find the difference in z (between linear projection from layer 2
1937 // and the actual point: Dzeta= z2/r2*r1 -z1)
1938 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1939 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1941 // make "elliptical" cut in Phi and Zeta!
1942 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1943 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1944 if (d<1) {ret=kTRUE; break;};
1950 //_________________________________________________________________________
1951 void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
1953 // this method apply a rotation by 180 degree around the Z (beam) axis to all
1954 // the RecPoints in a given layer to be used to build tracklets.
1955 // **************** VERY IMPORTANT:: ***************
1956 // It must be called just after LoadClusterArrays, since afterwards the datamember
1957 // fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
1958 // instead of Cartesian
1960 if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1961 AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1963 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1964 fClustersLay1[iC1][0]*=-1;
1965 fClustersLay1[iC1][1]*=-1;
1969 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1970 fClustersLay2[iC2][0]*=-1;
1971 fClustersLay2[iC2][1]*=-1;
1976 //____________________________________________________________________________
1977 Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *esd){
1978 // This method is used to find the tracklets.
1979 // It is called from AliReconstruction
1980 // The vertex is supposed to be associated to the Tracker (i.e. to this) already
1981 // The cluster is supposed to be associated to the Tracker already
1982 // In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
1985 // apply cuts on the vertex quality
1986 const AliESDVertex *vertex = esd->GetVertex();
1987 if(vertex->GetNContributors()<fMinContVtx) return 0;
1989 AliRunLoader* runLoader = AliRunLoader::Instance();
1991 Error("Clusters2Tracks", "no run loader found");
1994 AliStack *pStack=0x0; TTree *tRefTree=0x0;
1996 runLoader->LoadKinematics("read");
1997 runLoader->LoadTrackRefs("read");
1998 pStack= runLoader->Stack();
1999 tRefTree= runLoader->TreeTR();
2001 Reconstruct(pStack,tRefTree);
2003 if (GetLightBkgStudyInParallel()) {
2004 AliStack *dummy1=0x0; TTree *dummy2=0x0;
2005 ReflectClusterAroundZAxisForLayer(1);
2006 Reconstruct(dummy1,dummy2,kTRUE);
2010 //____________________________________________________________________________
2011 Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
2013 // It is called from AliReconstruction
2019 if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
2020 if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
2021 if(GetLightBkgStudyInParallel()) {
2022 TString name="AliITSPlaneEffSPDtrackletBkg.root";
2023 TFile* pefile = TFile::Open(name, "RECREATE");
2024 rc*=fPlaneEffBkg->Write();
2029 //____________________________________________________________________
2031 AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
2033 // - gets the clusters from the cluster tree
2034 // - convert them into global coordinates
2035 // - store them in the internal arrays
2036 // - count the number of cluster-fired chips
2038 //AliDebug(1,"Loading clusters and cluster-fired chips ...");
2043 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
2044 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
2046 itsClusterBranch->SetAddress(&itsClusters);
2048 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
2049 Float_t cluGlo[3]={0.,0.,0.};
2051 // loop over the its subdetectors
2052 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
2054 if (!itsClusterTree->GetEvent(iIts))
2057 Int_t nClusters = itsClusters->GetEntriesFast();
2059 // number of clusters in each chip of the current module
2062 // loop over clusters
2063 while(nClusters--) {
2064 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
2066 layer = cluster->GetLayer();
2067 if (layer>1) continue;
2069 cluster->GetGlobalXYZ(cluGlo);
2070 Float_t x = cluGlo[0];
2071 Float_t y = cluGlo[1];
2072 Float_t z = cluGlo[2];
2075 fClustersLay1[fNClustersLay1][0] = x;
2076 fClustersLay1[fNClustersLay1][1] = y;
2077 fClustersLay1[fNClustersLay1][2] = z;
2079 for (Int_t i=0; i<3; i++)
2080 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
2083 Int_t det=cluster->GetDetectorIndex();
2084 if(det<0 || det>79) {AliError("Cluster with det. index out of boundaries"); return;}
2085 fhClustersInModuleLay1[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2089 fClustersLay2[fNClustersLay2][0] = x;
2090 fClustersLay2[fNClustersLay2][1] = y;
2091 fClustersLay2[fNClustersLay2][2] = z;
2093 for (Int_t i=0; i<3; i++)
2094 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
2097 Int_t det=cluster->GetDetectorIndex();
2098 if(det<0 || det>159) {AliError("Cluster with det. index out of boundaries"); return;}
2099 fhClustersInModuleLay2[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2103 }// end of cluster loop
2105 } // end of its "subdetector" loop
2107 itsClusters->Delete();
2111 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
2113 //_________________________________________________________________________
2115 AliITSTrackleterSPDEff::SetLightBkgStudyInParallel(Bool_t b) {
2117 // - set Bool_t fLightBackgroundStudyInParallel = b
2118 // a) if you set this kTRUE, then the estimation of the
2119 // SPD efficiency is done as usual for data, but in
2120 // parallel a light (i.e. without control histograms, etc.)
2121 // evaluation of combinatorial background is performed
2122 // with the usual ReflectClusterAroundZAxisForLayer method.
2123 // b) if you set this kFALSE, then you would not have a second
2124 // container for PlaneEfficiency statistics to be used for background
2125 // (fPlaneEffBkg=0). If you want to have a full evaluation of the
2126 // background (with all control histograms and additional data
2127 // members referring to the background) then you have to call the
2128 // method SetReflectClusterAroundZAxisForLayer(kTRUE) esplicitily
2129 fLightBkgStudyInParallel=b;
2130 if(fLightBkgStudyInParallel) {
2131 if(!fPlaneEffBkg) fPlaneEffBkg = new AliITSPlaneEffSPD();
2134 delete fPlaneEffBkg;
2138 //______________________________________________________________
2139 void AliITSTrackleterSPDEff::SetReflectClusterAroundZAxisForLayer(Int_t ilayer,Bool_t b){
2141 // method to study residual background:
2142 // Input b= KTRUE --> reflect the clusters
2143 // ilayer (either 0 or 1) --> which SPD layers should be reflected
2145 if(b) {AliInfo(Form("All clusters on layer %d will be rotated by 180 deg around z",ilayer));
2146 SetLightBkgStudyInParallel(kFALSE);}
2147 if(ilayer==0) fReflectClusterAroundZAxisForLayer0=b; // a rotation by 180degree around the Z axis
2148 else if(ilayer==1) fReflectClusterAroundZAxisForLayer1=b; // (x->-x; y->-y) to all RecPoints on a
2149 else AliInfo("Nothing done: input argument (ilayer) either 0 or 1"); // given layer is applied. In such a way