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 SetOnlyOneTrackletPerC2();
142 fClustersLay1 = new Float_t*[300000];
143 fClustersLay2 = new Float_t*[300000];
144 fTracklets = new Float_t*[300000];
145 fAssociationFlag = new Bool_t[300000];
149 SetOnlyOneTrackletPerC1();
151 fAssociationFlag1 = new Bool_t[300000];
152 fChipPredOnLay2 = new UInt_t[300000];
153 fChipPredOnLay1 = new UInt_t[300000];
154 fChipUpdatedInEvent = new Bool_t[1200];
156 for(Int_t i=0; i<300000; i++) {
157 // from AliITSMultReconstructor
158 fClustersLay1[i] = new Float_t[6];
159 fClustersLay2[i] = new Float_t[6];
160 fTracklets[i] = new Float_t[5];
161 fAssociationFlag[i] = kFALSE;
163 fAssociationFlag1[i] = kFALSE;
165 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
167 if (GetHistOn()) BookHistos();
169 fPlaneEffSPD = new AliITSPlaneEffSPD();
170 SetLightBkgStudyInParallel();
172 //______________________________________________________________________
173 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
175 // from AliITSMultReconstructor
176 fClustersLay1(mr.fClustersLay1),
177 fClustersLay2(mr.fClustersLay2),
178 fTracklets(mr.fTracklets),
179 fAssociationFlag(mr.fAssociationFlag),
180 fNClustersLay1(mr.fNClustersLay1),
181 fNClustersLay2(mr.fNClustersLay2),
182 fNTracklets(mr.fNTracklets),
183 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
184 fPhiWindowL2(mr.fPhiWindowL2),
185 fZetaWindowL2(mr.fZetaWindowL2),
186 fPhiOverlapCut(mr.fPhiOverlapCut),
187 fZetaOverlapCut(mr.fZetaOverlapCut),
189 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
190 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
191 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
192 fhClustersDPhiAll(mr.fhClustersDPhiAll),
193 fhClustersDThetaAll(mr.fhClustersDThetaAll),
194 fhClustersDZetaAll(mr.fhClustersDZetaAll),
195 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
196 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
197 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
198 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
199 fhetaTracklets(mr.fhetaTracklets),
200 fhphiTracklets(mr.fhphiTracklets),
201 fhetaClustersLay1(mr.fhetaClustersLay1),
202 fhphiClustersLay1(mr.fhphiClustersLay1),
204 fAssociationFlag1(mr.fAssociationFlag1),
205 fChipPredOnLay2(mr.fChipPredOnLay2),
206 fChipPredOnLay1(mr.fChipPredOnLay1),
207 fNTracklets1(mr.fNTracklets1),
208 fPhiWindowL1(mr.fPhiWindowL1),
209 fZetaWindowL1(mr.fZetaWindowL1),
210 fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
211 fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
212 fMinContVtx(mr.fMinContVtx),
213 fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
214 fPlaneEffSPD(mr.fPlaneEffSPD),
215 fPlaneEffBkg(mr.fPlaneEffBkg),
216 fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
217 fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
218 fLightBkgStudyInParallel(mr.fLightBkgStudyInParallel),
220 fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
221 fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
222 fUseOnlySameParticle(mr.fUseOnlySameParticle),
223 fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
224 fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
225 fPredictionPrimary(mr.fPredictionPrimary),
226 fPredictionSecondary(mr.fPredictionSecondary),
227 fClusterPrimary(mr.fClusterPrimary),
228 fClusterSecondary(mr.fClusterSecondary),
229 fSuccessPP(mr.fSuccessPP),
230 fSuccessTT(mr.fSuccessTT),
231 fSuccessS(mr.fSuccessS),
232 fSuccessP(mr.fSuccessP),
233 fFailureS(mr.fFailureS),
234 fFailureP(mr.fFailureP),
236 fNonRecons(mr.fNonRecons),
237 fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
238 fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
239 fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
240 fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
241 fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
242 fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
243 fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
244 fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
245 fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
246 fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
247 fhetaClustersLay2(mr.fhetaClustersLay2),
248 fhphiClustersLay2(mr.fhphiClustersLay2),
249 fhClustersInChip(mr.fhClustersInChip),
250 fhClustersInModuleLay1(mr.fhClustersInModuleLay1),
251 fhClustersInModuleLay2(mr.fhClustersInModuleLay2)
256 //______________________________________________________________________
257 AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
258 // Assignment operator
259 this->~AliITSTrackleterSPDEff();
260 new(this) AliITSTrackleterSPDEff(mr);
263 //______________________________________________________________________
264 AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
266 // from AliITSMultReconstructor
268 for(Int_t i=0; i<300000; i++) {
269 delete [] fClustersLay1[i];
270 delete [] fClustersLay2[i];
271 delete [] fTracklets[i];
273 delete [] fClustersLay1;
274 delete [] fClustersLay2;
275 delete [] fTracklets;
276 delete [] fAssociationFlag;
281 delete [] fAssociationFlag1;
283 delete [] fChipPredOnLay2;
284 delete [] fChipPredOnLay1;
286 delete [] fChipUpdatedInEvent;
288 delete [] fPredictionPrimary;
289 delete [] fPredictionSecondary;
290 delete [] fClusterPrimary;
291 delete [] fClusterSecondary;
292 delete [] fSuccessPP;
293 delete [] fSuccessTT;
299 delete [] fNonRecons;
310 //____________________________________________________________________
312 AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef, Bool_t lbkg) {
314 // - you have to take care of the following, before of using Reconstruct
315 // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
316 // and convert the cluster coordinates to theta, phi (seen from the
317 // interaction vertex).
318 // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
319 // - Find the extrapolation/interpolation point.
320 // - Find the chip corresponding to that
321 // - Check if there is a cluster near that point
324 if(lbkg && !GetLightBkgStudyInParallel()) {
325 AliError("You asked for lightBackground in the Reconstruction without proper call to SetLightBkgStudyInParallel(1)");
328 AliITSPlaneEffSPD *pe;
335 // retrieve the vertex position
337 vtx[0]=(Float_t)GetX();
338 vtx[1]=(Float_t)GetY();
339 vtx[2]=(Float_t)GetZ();
340 // to study residual background (i.e. contribution from TT' to measured efficiency)
341 if(fReflectClusterAroundZAxisForLayer0 && !lbkg) ReflectClusterAroundZAxisForLayer(0);
342 if(fReflectClusterAroundZAxisForLayer1 && !lbkg) ReflectClusterAroundZAxisForLayer(1);
344 if(fMC && !pStack && !lbkg) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
345 if(fMC && !tRef && !lbkg) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
347 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
348 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
349 Int_t nfClu1=0; Int_t ntClu1=0;
350 Int_t nfClu2=0; Int_t ntClu2=0;
352 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
353 // for this new event)
354 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
356 // find the tracklets
357 AliDebug(1,"Looking for tracklets... ");
358 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
360 //###########################################################
361 // Loop on layer 1 : finding theta, phi and z
363 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
364 Float_t x = fClustersLay1[iC1][0] - vtx[0];
365 Float_t y = fClustersLay1[iC1][1] - vtx[1];
366 Float_t z = fClustersLay1[iC1][2] - vtx[2];
368 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
370 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
371 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
372 fClustersLay1[iC1][2] = z; // Store z
374 // find the Radius and the chip corresponding to the extrapolation point
376 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
378 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
379 key=999999; // also some other actions should be taken if not Found
381 nfTraPred2+=(Int_t)found; // this for debugging purpose
382 ntTraPred2++; // to check efficiency of the method FindChip
383 fChipPredOnLay2[iC1] = key;
384 fAssociationFlag1[iC1] = kFALSE;
386 if (fHistOn && !lbkg) {
387 Float_t eta=fClustersLay1[iC1][0];
388 eta= TMath::Tan(eta/2.);
389 eta=-TMath::Log(eta);
390 fhetaClustersLay1->Fill(eta);
391 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
392 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
395 // Loop on layer 2 : finding theta, phi and r
396 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
397 Float_t x = fClustersLay2[iC2][0] - vtx[0];
398 Float_t y = fClustersLay2[iC2][1] - vtx[1];
399 Float_t z = fClustersLay2[iC2][2] - vtx[2];
401 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
403 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
404 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
405 fClustersLay2[iC2][2] = z; // Store z
407 // find the Radius and the chip corresponding to the extrapolation point
409 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
411 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
414 nfTraPred1+=(Int_t)found; // this for debugging purpose
415 ntTraPred1++; // to check efficiency of the method FindChip
416 fChipPredOnLay1[iC2] = key;
417 fAssociationFlag[iC2] = kFALSE;
419 if (fHistOn && !lbkg) {
420 Float_t eta=fClustersLay2[iC2][0];
421 eta= TMath::Tan(eta/2.);
422 eta=-TMath::Log(eta);
423 fhetaClustersLay2->Fill(eta);
424 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
425 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
429 //###########################################################
431 // First part : Extrapolation to Layer 2
434 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
436 // here the control to check whether the efficiency of the chip traversed by this tracklet
437 // prediction has already been updated in this event using another tracklet prediction
438 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
440 // reset of variables for multiple candidates
441 Int_t iC2WithBestDist = 0; // reset
442 Float_t distmin = 100.; // just to put a huge number!
443 Float_t dPhimin = 0.; // Used for histograms only!
444 Float_t dThetamin = 0.; // Used for histograms only!
445 Float_t dZetamin = 0.; // Used for histograms only!
447 // in any case, if MC has been required, store statistics of primaries and secondaries
448 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
450 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
451 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
452 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
453 // do it always as a function of the chip number used to built the prediction
454 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
455 if (!found) {AliWarning(
456 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
458 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
459 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
460 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
461 { // this cluster is from a primary particle
462 fClusterPrimary[key]++;
464 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
465 } else { // this cluster is from a secondary particle
466 fClusterSecondary[key]++;
468 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
471 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
472 // (in case the prediction is reliable)
473 if( fChipPredOnLay2[iC1]<1200) {
474 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
475 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
476 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
477 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
478 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
479 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
480 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
481 else fNonRecons[fChipPredOnLay2[iC1]]++;
486 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
488 // The following excludes double associations
489 if (!fAssociationFlag[iC2]) {
491 // find the difference in angles
492 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
493 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
494 // take into account boundary condition
495 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
497 // find the difference in z (between linear projection from layer 1
498 // and the actual point: Dzeta= z1/r1*r2 -z2)
499 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
500 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
502 if (fHistOn && !lbkg) {
503 fhClustersDPhiAll->Fill(dPhi);
504 fhClustersDThetaAll->Fill(dTheta);
505 fhClustersDZetaAll->Fill(dZeta);
506 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
507 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
510 // make "elliptical" cut in Phi and Zeta!
511 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
512 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
516 //look for the minimum distance: the minimum is in iC2WithBestDist
517 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
518 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
522 iC2WithBestDist = iC2;
525 } // end of loop over clusters in layer 2
527 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
529 if (fHistOn && !lbkg) {
530 fhClustersDPhiAcc->Fill(dPhimin);
531 fhClustersDThetaAcc->Fill(dThetamin);
532 fhClustersDZetaAcc->Fill(dZetamin);
533 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
534 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
537 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
538 // flag the association
540 // store the tracklet
542 // use the theta from the clusters in the first layer
543 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
544 // use the phi from the clusters in the first layer
545 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
546 // Store the difference between phi1 and phi2
547 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
554 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
566 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
570 fTracklets[fNTracklets][3] = -2;
573 if (fHistOn && !lbkg) {
574 Float_t eta=fTracklets[fNTracklets][0];
575 eta= TMath::Tan(eta/2.);
576 eta=-TMath::Log(eta);
577 fhetaTracklets->Fill(eta);
578 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
581 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
582 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
585 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
588 nfClu2+=(Int_t)found; // this for debugging purpose
589 ntClu2++; // to check efficiency of the method FindChip
590 if(key<1200) { // the Chip has been found
591 if(fMC && !lbkg) { // this part only for MC
592 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
593 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
594 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
597 if(primary) fSuccessPP[key]++;
599 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
600 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
603 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
605 pe->UpDatePlaneEff(kTRUE,key); // success
606 fChipUpdatedInEvent[key]=kTRUE;
608 if(primary) fSuccessP[key]++;
609 if(secondary) fSuccessS[key]++;
613 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
614 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
616 if(primary) fSuccessP[key]++;
617 if(secondary) fSuccessS[key]++;
624 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
625 else if (fChipPredOnLay2[iC1]<1200) {
626 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
627 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
629 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
630 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
633 } // end of loop over clusters in layer 1
635 fNTracklets1=fNTracklets;
637 //###################################################################
639 // Second part : Interpolation to Layer 1
642 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
644 // here the control to check whether the efficiency of the chip traversed by this tracklet
645 // prediction has already been updated in this event using another tracklet prediction
646 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
648 // reset of variables for multiple candidates
649 Int_t iC1WithBestDist = 0; // reset
650 Float_t distmin = 100.; // just to put a huge number!
651 Float_t dPhimin = 0.; // Used for histograms only!
652 Float_t dThetamin = 0.; // Used for histograms only!
653 Float_t dZetamin = 0.; // Used for histograms only!
655 // in any case, if MC has been required, store statistics of primaries and secondaries
656 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
658 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
659 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
660 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
661 // do it always as a function of the chip number used to built the prediction
662 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
663 if (!found) {AliWarning(
664 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
666 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
667 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
668 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
669 { // this cluster is from a primary particle
670 fClusterPrimary[key]++;
672 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
673 } else { // this cluster is from a secondary particle
674 fClusterSecondary[key]++;
676 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
679 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
680 // (in case the prediction is reliable)
681 if( fChipPredOnLay1[iC2]<1200) {
682 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
683 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
684 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
685 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
686 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
687 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
688 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
689 else fNonRecons[fChipPredOnLay1[iC2]]++;
694 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
696 // The following excludes double associations
697 if (!fAssociationFlag1[iC1]) {
699 // find the difference in angles
700 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
701 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
702 // take into account boundary condition
703 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
706 // find the difference in z (between linear projection from layer 2
707 // and the actual point: Dzeta= z2/r2*r1 -z1)
708 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
709 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
712 if (fHistOn && !lbkg) {
713 fhClustersDPhiInterpAll->Fill(dPhi);
714 fhClustersDThetaInterpAll->Fill(dTheta);
715 fhClustersDZetaInterpAll->Fill(dZeta);
716 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
717 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
719 // make "elliptical" cut in Phi and Zeta!
720 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
721 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
725 //look for the minimum distance: the minimum is in iC1WithBestDist
726 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
727 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
731 iC1WithBestDist = iC1;
734 } // end of loop over clusters in layer 1
736 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
738 if (fHistOn && !lbkg) {
739 fhClustersDPhiInterpAcc->Fill(dPhimin);
740 fhClustersDThetaInterpAcc->Fill(dThetamin);
741 fhClustersDZetaInterpAcc->Fill(dZetamin);
742 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
743 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
746 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
747 // flag the association
749 // store the tracklet
751 // use the theta from the clusters in the first layer
752 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
753 // use the phi from the clusters in the first layer
754 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
755 // Store the difference between phi1 and phi2
756 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
763 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
775 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
779 fTracklets[fNTracklets][3] = -2;
782 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
783 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
786 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
789 nfClu1+=(Int_t)found; // this for debugging purpose
790 ntClu1++; // to check efficiency of the method FindChip
792 if(fMC && !lbkg) { // this part only for MC
793 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
794 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
795 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
796 if (label2 < 3) { // same label
798 if(primary) fSuccessPP[key]++;
800 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
801 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
804 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
806 pe->UpDatePlaneEff(kTRUE,key); // success
807 fChipUpdatedInEvent[key]=kTRUE;
809 if(primary) fSuccessP[key]++;
810 if(secondary) fSuccessS[key]++;
813 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
814 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
816 if(primary) fSuccessP[key]++;
817 if(secondary) fSuccessS[key]++;
824 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
825 else if (fChipPredOnLay1[iC2]<1200) {
826 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
827 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
829 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
830 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
833 } // end of loop over clusters in layer 2
835 AliDebug(1,Form("%d tracklets found", fNTracklets));
836 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
837 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
838 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
839 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
841 //____________________________________________________________________
842 Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
843 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
845 // Input: a) layer number in the range [0,1]
846 // b) vtx[3]: actual vertex
847 // c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
848 // d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
850 // Output: Unique key to locate a chip
851 // return: kTRUE if succesfull
853 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
854 Double_t r=GetRLayer(layer);
855 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
857 // set phiVtx in the range [0,2pi]
858 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
860 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
861 // of the intersection of the tracklet with the pixel layer.
862 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
864 if(TMath::Abs(thetaVtx)<1E-6) return kFALSE;
865 zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
867 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
868 Double_t vtxy[2]={vtx[0],vtx[1]};
869 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
870 // this method gives you two interceptions
871 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
872 // set phiAbs in the range [0,2pi]
873 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
874 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
875 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
876 // taking the closest one to phiVtx
877 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
878 } else phiAbs=phiVtx;
879 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
881 // now you need to locate the chip within the idet detector,
882 // starting from the local coordinates in such a detector
884 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
885 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
886 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
888 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
891 //______________________________________________________________________________
892 Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
894 // Return the average radius of a layer from Geometry
896 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
897 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
899 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
900 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
901 Double_t r=TMath::Sqrt(x*x + y*y);
903 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
904 r += TMath::Sqrt(x*x + y*y);
905 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
906 r += TMath::Sqrt(x*x + y*y);
907 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
908 r += TMath::Sqrt(x*x + y*y);
912 //______________________________________________________________________________
913 Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
914 Float_t &xloc, Float_t &zloc) {
915 // this method transform Global Cilindrical coordinates into local (i.e. module)
916 // cartesian coordinates
918 //Compute Cartesian Global Coordinate
919 Double_t xyzGlob[3],xyzLoc[3];
921 xyzGlob[0]=r*TMath::Cos(phi);
922 xyzGlob[1]=r*TMath::Sin(phi);
927 if(idet<0) return kFALSE;
929 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
930 Int_t lad = Int_t(idet/ndet) + 1;
931 Int_t det = idet - (lad-1)*ndet + 1;
933 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
935 xloc = (Float_t)xyzLoc[0];
936 zloc = (Float_t)xyzLoc[2];
940 //______________________________________________________________________________
941 Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
942 //--------------------------------------------------------------------
943 // This function finds the detector crossed by the track
944 // Input: layer in range [0,1]
945 // phi in ALICE absolute reference system
947 //--------------------------------------------------------------------
948 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
949 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
950 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
951 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
953 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
954 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
955 Double_t phiOffset=TMath::ATan2(y,x);
959 if (zOffset<0) // old geometry
960 dphi = -(phi-phiOffset);
962 dphi = phi-phiOffset;
964 if (dphi < 0) dphi += 2*TMath::Pi();
965 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
966 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
967 if (np>=nladders) np-=nladders;
968 if (np<0) np+=nladders;
970 Double_t dz=zOffset-z;
971 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
972 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
973 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
974 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
976 return np*ndetectors + nz;
978 //____________________________________________________________
979 Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
980 // this method find the intersection in xy between a tracklet (straight line) and
981 // a circonference (r=R), using polar coordinates.
983 Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
984 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
985 - R: radius of the circle
986 Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
988 Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
989 r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
990 In the same system, the equation of a semi-line is: phi=phiVtx;
991 Hence you get one interception only: P=(r,phiVtx)
992 Finally you want P in the ABSOLUTE ALICE reference system.
994 Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
995 Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
996 Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
997 Double_t cC=rO*rO-R*R;
998 Double_t dDelta=bB*bB-4*cC;
999 if(dDelta<0) return kFALSE;
1001 r1=(-bB-TMath::Sqrt(dDelta))/2;
1002 r2=(-bB+TMath::Sqrt(dDelta))/2;
1003 if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
1004 Double_t r=TMath::Max(r1,r2); // take the positive
1005 Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
1006 Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
1007 pvtx[0]=r*TMath::Cos(phiVtx);
1008 pvtx[1]=r*TMath::Sin(phiVtx);
1009 pP[0]=vtx[0]+pvtx[0];
1010 pP[1]=vtx[1]+pvtx[1];
1011 phi=TMath::ATan2(pP[1],pP[0]);
1014 //___________________________________________________________
1015 Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
1017 // simple method to reduce all angles (in rad)
1021 while(angle >=2*TMath::Pi() || angle<0) {
1022 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
1023 if(angle < 0) angle+=2*TMath::Pi();
1027 //___________________________________________________________
1028 Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
1030 // This method check if a particle is primary; i.e.
1031 // it comes from the main vertex and it is a "stable" particle, according to
1032 // AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
1033 // a stable particle: it has no effect on this analysis).
1034 // This method can be called only for MC events, where Kinematics is available.
1035 // if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
1036 // returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
1037 // The latter (see below) try to verify if a primary particle is also "detectable".
1039 if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
1040 if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
1041 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
1042 // return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
1043 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
1044 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
1045 TParticle* part = stack->Particle(ipart);
1046 TParticle* part0 = stack->Particle(0); // first primary
1047 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
1048 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
1049 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
1051 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
1052 TParticlePDG* pdgPart = part->GetPDG();
1053 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
1055 Double_t distx = part->Vx() - part0->Vx();
1056 Double_t disty = part->Vy() - part0->Vy();
1057 Double_t distz = part->Vz() - part0->Vz();
1058 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
1060 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
1061 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
1062 return kFALSE; }// primary if within 500 microns from true Vertex
1064 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
1067 //_____________________________________________________________________________________________
1068 Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
1070 // This private method can be applied on MC particles (if stack is available),
1071 // provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
1073 // It define "detectable" a primary particle according to the following criteria:
1075 // - if no decay products can be found in the stack (note that this does not
1076 // means it is stable, since a particle is stored in stack if it has at least 1 hit in a
1077 // sensitive detector)
1078 // - if it has at least one decay daughter produced outside or just on the outer pixel layer
1079 // - if the last decay particle is an electron (or a muon) which is not produced in-between
1080 // the two pixel layers (this is likely to be a kink).
1081 if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
1082 if(!stack) {AliError("null pointer to MC stack"); return 0;}
1083 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
1085 TParticle* part = stack->Particle(ipart);
1086 //TParticle* part0 = stack->Particle(0); // first primary
1092 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
1093 // its real fate ! But you have to take it !
1094 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
1095 Int_t lastDau = part->GetLastDaughter();
1096 nDau = lastDau - firstDau + 1;
1097 Double_t distMax=0.;
1099 for(Int_t j=firstDau; j<=lastDau; j++) {
1100 dau = stack->Particle(j);
1101 Double_t distx = dau->Vx();
1102 Double_t disty = dau->Vy();
1103 //Double_t distz = dau->Vz();
1104 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
1105 if(distR<distMax) continue; // considere only the daughter produced at largest radius
1109 dau = stack->Particle(jmax);
1110 pdgDau=dau->GetPdgCode();
1111 if (pdgDau == 11 || pdgDau == 13 ) {
1112 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
1113 else nret =0; // delta-ray emission in material (keep it)
1115 else {// not ele or muon
1116 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
1120 //_________________________________________________________________
1121 void AliITSTrackleterSPDEff::InitPredictionMC() {
1123 // this method allocate memory for the MC related informations
1124 // all the counters are set to 0
1127 if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1128 fPredictionPrimary = new Int_t[1200];
1129 fPredictionSecondary = new Int_t[1200];
1130 fClusterPrimary = new Int_t[1200];
1131 fClusterSecondary = new Int_t[1200];
1132 fSuccessPP = new Int_t[1200];
1133 fSuccessTT = new Int_t[1200];
1134 fSuccessS = new Int_t[1200];
1135 fSuccessP = new Int_t[1200];
1136 fFailureS = new Int_t[1200];
1137 fFailureP = new Int_t[1200];
1138 fRecons = new Int_t[1200];
1139 fNonRecons = new Int_t[1200];
1140 for(Int_t i=0; i<1200; i++) {
1141 fPredictionPrimary[i]=0;
1142 fPredictionSecondary[i]=0;
1143 fPredictionSecondary[i]=0;
1144 fClusterSecondary[i]=0;
1156 //_________________________________________________________________
1157 void AliITSTrackleterSPDEff::DeletePredictionMC() {
1159 // this method deallocate memory for the MC related informations
1160 // all the counters are set to 0
1163 if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1164 if(fPredictionPrimary) {
1165 delete fPredictionPrimary; fPredictionPrimary=0;
1167 if(fPredictionSecondary) {
1168 delete fPredictionSecondary; fPredictionSecondary=0;
1170 if(fClusterPrimary) {
1171 delete fClusterPrimary; fClusterPrimary=0;
1173 if(fClusterSecondary) {
1174 delete fClusterSecondary; fClusterSecondary=0;
1177 delete fSuccessPP; fSuccessPP=0;
1180 delete fSuccessTT; fSuccessTT=0;
1183 delete fSuccessS; fSuccessS=0;
1186 delete fSuccessP; fSuccessP=0;
1189 delete fFailureS; fFailureS=0;
1192 delete fFailureP; fFailureP=0;
1195 delete fRecons; fRecons=0;
1198 delete fNonRecons; fNonRecons=0;
1202 //______________________________________________________________________
1203 Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
1205 // This method return the Data menmber fPredictionPrimary [1200].
1206 // You can call it only for MC events.
1207 // fPredictionPrimary[key] contains the number of tracklet predictions on the
1208 // given chip key built using a cluster on the other layer produced (at least)
1209 // from a primary particle.
1210 // Key refers to the chip crossed by the prediction
1213 if (!fMC) {CallWarningMC(); return 0;}
1214 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1215 return fPredictionPrimary[(Int_t)key];
1217 //______________________________________________________________________
1218 Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
1220 // This method return the Data menmber fPredictionSecondary [1200].
1221 // You can call it only for MC events.
1222 // fPredictionSecondary[key] contains the number of tracklet predictions on the
1223 // given chip key built using a cluster on the other layer produced (only)
1224 // from a secondary particle
1225 // Key refers to the chip crossed by the prediction
1228 if (!fMC) {CallWarningMC(); return 0;}
1229 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1230 return fPredictionSecondary[(Int_t)key];
1232 //______________________________________________________________________
1233 Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
1235 // This method return the Data menmber fClusterPrimary [1200].
1236 // You can call it only for MC events.
1237 // fClusterPrimary[key] contains the number of tracklet predictions
1238 // built using a cluster on that layer produced (only)
1239 // from a primary particle
1240 // Key refers to the chip used to build the prediction
1243 if (!fMC) {CallWarningMC(); return 0;}
1244 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1245 return fClusterPrimary[(Int_t)key];
1247 //______________________________________________________________________
1248 Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
1250 // This method return the Data menmber fClusterSecondary [1200].
1251 // You can call it only for MC events.
1252 // fClusterSecondary[key] contains the number of tracklet predictions
1253 // built using a cluster on that layer produced (only)
1254 // from a secondary particle
1255 // Key refers to the chip used to build the prediction
1257 if (!fMC) {CallWarningMC(); return 0;}
1258 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1259 return fClusterSecondary[(Int_t)key];
1261 //______________________________________________________________________
1262 Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1264 // This method return the Data menmber fSuccessPP [1200].
1265 // You can call it only for MC events.
1266 // fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1267 // with a cluster on the other layer) built by using the same primary particle
1268 // the unique chip key refers to the chip which get updated its efficiency
1270 if (!fMC) {CallWarningMC(); return 0;}
1271 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1272 return fSuccessPP[(Int_t)key];
1274 //______________________________________________________________________
1275 Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1277 // This method return the Data menmber fSuccessTT [1200].
1278 // You can call it only for MC events.
1279 // fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1280 // with a cluster on the other layer) built by using the same particle (whatever)
1281 // the unique chip key refers to the chip which get updated its efficiency
1283 if (!fMC) {CallWarningMC(); return 0;}
1284 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1285 return fSuccessTT[(Int_t)key];
1287 //______________________________________________________________________
1288 Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1290 // This method return the Data menmber fSuccessS [1200].
1291 // You can call it only for MC events.
1292 // fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1293 // with a cluster on the other layer) built by using a secondary particle
1294 // the unique chip key refers to the chip which get updated its efficiency
1296 if (!fMC) {CallWarningMC(); return 0;}
1297 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1298 return fSuccessS[(Int_t)key];
1300 //______________________________________________________________________
1301 Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1303 // This method return the Data menmber fSuccessP [1200].
1304 // You can call it only for MC events.
1305 // fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1306 // with a cluster on the other layer) built by using a primary particle
1307 // the unique chip key refers to the chip which get updated its efficiency
1309 if (!fMC) {CallWarningMC(); return 0;}
1310 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1311 return fSuccessP[(Int_t)key];
1313 //______________________________________________________________________
1314 Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1316 // This method return the Data menmber fFailureS [1200].
1317 // You can call it only for MC events.
1318 // fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1319 // with a cluster on the other layer) built by using a secondary particle
1320 // the unique chip key refers to the chip which get updated its efficiency
1322 if (!fMC) {CallWarningMC(); return 0;}
1323 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1324 return fFailureS[(Int_t)key];
1326 //______________________________________________________________________
1327 Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1329 // This method return the Data menmber fFailureP [1200].
1330 // You can call it only for MC events.
1331 // fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1332 // with a cluster on the other layer) built by using a primary particle
1333 // the unique chip key refers to the chip which get updated its efficiency
1335 if (!fMC) {CallWarningMC(); return 0;}
1336 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1337 return fFailureP[(Int_t)key];
1339 //_____________________________________________________________________
1340 Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1342 // This method return the Data menmber fRecons [1200].
1343 // You can call it only for MC events.
1344 // fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1345 // has an hit in the detector)
1346 // the unique chip key refers to the chip where fall the prediction
1348 if (!fMC) {CallWarningMC(); return 0;}
1349 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1350 return fRecons[(Int_t)key];
1352 //_____________________________________________________________________
1353 Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1355 // This method return the Data menmber fNonRecons [1200].
1356 // You can call it only for MC events.
1357 // fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1358 // has not any hit in the detector)
1359 // the unique chip key refers to the chip where fall the prediction
1361 if (!fMC) {CallWarningMC(); return 0;}
1362 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1363 return fNonRecons[(Int_t)key];
1365 //______________________________________________________________________
1366 void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1367 // Print out some class data values in Ascii Form to output stream
1369 // ostream *os Output stream where Ascii data is to be writen
1374 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindowL2 <<" "<< fZetaWindowL2
1375 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1376 << " " << fUpdateOncePerEventPlaneEff << " " << fMinContVtx
1377 << " " << fReflectClusterAroundZAxisForLayer0
1378 << " " << fReflectClusterAroundZAxisForLayer1;
1380 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1381 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1382 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1383 << " " << fUseOnlyStableParticle ;
1384 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1385 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1386 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1387 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
1388 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1389 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1390 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1391 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1392 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1393 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1394 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1395 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
1398 //______________________________________________________________________
1399 void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1400 // Read in some class data values in Ascii Form to output stream
1402 // istream *is Input stream where Ascii data is to be read in from
1409 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindowL2 >> fZetaWindowL2
1410 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1411 >> fUpdateOncePerEventPlaneEff >> fMinContVtx
1412 >> fReflectClusterAroundZAxisForLayer0
1413 >> fReflectClusterAroundZAxisForLayer1;
1414 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
1416 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1418 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1419 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1420 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1421 >> fUseOnlyStableParticle;
1422 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1423 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1424 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1425 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1426 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1427 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1428 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1429 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1430 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1431 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1432 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1433 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1437 //______________________________________________________________________
1438 ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1439 // Standard output streaming function
1441 // ostream &os output steam
1442 // AliITSTrackleterSPDEff &s class to be streamed.
1446 // ostream &os The stream pointer
1451 //______________________________________________________________________
1452 istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1453 // Standard inputput streaming function
1455 // istream &is input steam
1456 // AliITSTrackleterSPDEff &s class to be streamed.
1460 // ostream &os The stream pointer
1462 //printf("prova %d \n", (Int_t)s.GetMC());
1466 //______________________________________________________________________
1467 void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
1469 // This Method write into an either asci or root file
1470 // the used cuts and the statistics of the MC related quantities
1471 // The method SetMC() has to be called before
1472 // Input TString filename: name of file for output (it deletes already existing
1477 //if(!fMC) {CallWarningMC(); return;}
1478 if (!filename.Contains(".root")) {
1479 ofstream out(filename.Data(),ios::out | ios::binary);
1485 TFile* mcfile = TFile::Open(filename, "RECREATE");
1486 TH1F* cuts = new TH1F("cuts", "list of cuts", 11, 0, 11); // TH1I containing cuts
1487 cuts->SetBinContent(1,fPhiWindowL1);
1488 cuts->SetBinContent(2,fZetaWindowL1);
1489 cuts->SetBinContent(3,fPhiWindowL2);
1490 cuts->SetBinContent(4,fZetaWindowL2);
1491 cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
1492 cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
1493 cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
1494 cuts->SetBinContent(8,fMinContVtx);
1495 cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer0);
1496 cuts->SetBinContent(10,fReflectClusterAroundZAxisForLayer1);
1497 cuts->SetBinContent(11,fMC);
1500 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
1502 TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
1503 mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
1504 mc0->SetBinContent(2,fUseOnlySecondaryForPred);
1505 mc0->SetBinContent(3,fUseOnlySameParticle);
1506 mc0->SetBinContent(4,fUseOnlyDifferentParticle);
1507 mc0->SetBinContent(5,fUseOnlyStableParticle);
1511 mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
1512 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
1514 mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
1515 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
1517 mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
1518 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
1520 mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
1521 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
1523 mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
1524 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
1526 mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
1527 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
1529 mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
1530 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
1532 mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
1533 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
1535 mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
1536 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
1538 mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
1539 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
1541 mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
1542 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
1544 mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
1545 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
1553 //____________________________________________________________________
1554 void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
1556 // This Method read from an asci file (do not know why binary does not work)
1557 // the cuts to be used and the statistics of the MC related quantities
1558 // Input TString filename: name of input file for output
1559 // The method SetMC() has to be called before
1563 //if(!fMC) {CallWarningMC(); return;}
1564 if( gSystem->AccessPathName( filename.Data() ) ) {
1565 AliError( Form( "file (%s) not found", filename.Data() ) );
1569 if (!filename.Contains(".root")) {
1570 ifstream in(filename.Data(),ios::in | ios::binary);
1577 TFile *mcfile = TFile::Open(filename);
1578 TH1F *cuts = (TH1F*)mcfile->Get("cuts");
1579 fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
1580 fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
1581 fPhiWindowL2=(Float_t)cuts->GetBinContent(3);
1582 fZetaWindowL2=(Float_t)cuts->GetBinContent(4);
1583 fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
1584 fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
1585 fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
1586 fMinContVtx=(Int_t)cuts->GetBinContent(8);
1587 fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(9);
1588 fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(10);
1589 fMC=(Bool_t)cuts->GetBinContent(11);
1590 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1591 else { // only if file with MC predictions
1592 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1593 TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
1594 fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
1595 fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
1596 fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
1597 fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
1598 fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
1600 mc1 =(TH1I*)mcfile->Get("mc1");
1601 for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1602 mc1 =(TH1I*)mcfile->Get("mc2");
1603 for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1604 mc1 =(TH1I*)mcfile->Get("mc3");
1605 for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1606 mc1 =(TH1I*)mcfile->Get("mc4");
1607 for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1608 mc1 =(TH1I*)mcfile->Get("mc5");
1609 for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1610 mc1 =(TH1I*)mcfile->Get("mc6");
1611 for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
1612 mc1 =(TH1I*)mcfile->Get("mc7");
1613 for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1614 mc1 =(TH1I*)mcfile->Get("mc8");
1615 for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1616 mc1 =(TH1I*)mcfile->Get("mc9");
1617 for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1618 mc1 =(TH1I*)mcfile->Get("mc10");
1619 for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1620 mc1 =(TH1I*)mcfile->Get("mc11");
1621 for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1622 mc1 =(TH1I*)mcfile->Get("mc12");
1623 for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1629 //____________________________________________________________________
1630 Bool_t AliITSTrackleterSPDEff::SaveHists() {
1631 // This (private) method save the histograms on the output file
1632 // (only if fHistOn is TRUE).
1633 // Also the histograms from the base class are saved through the
1634 // AliITSMultReconstructor::SaveHists() call
1636 if (!GetHistOn()) return kFALSE;
1638 // AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1639 fhClustersDPhiAll->Write();
1640 fhClustersDThetaAll->Write();
1641 fhClustersDZetaAll->Write();
1642 fhDPhiVsDThetaAll->Write();
1643 fhDPhiVsDZetaAll->Write();
1645 fhClustersDPhiAcc->Write();
1646 fhClustersDThetaAcc->Write();
1647 fhClustersDZetaAcc->Write();
1648 fhDPhiVsDThetaAcc->Write();
1649 fhDPhiVsDZetaAcc->Write();
1651 fhetaTracklets->Write();
1652 fhphiTracklets->Write();
1653 fhetaClustersLay1->Write();
1654 fhphiClustersLay1->Write();
1656 fhClustersDPhiInterpAll->Write();
1657 fhClustersDThetaInterpAll->Write();
1658 fhClustersDZetaInterpAll->Write();
1659 fhDPhiVsDThetaInterpAll->Write();
1660 fhDPhiVsDZetaInterpAll->Write();
1662 fhClustersDPhiInterpAcc->Write();
1663 fhClustersDThetaInterpAcc->Write();
1664 fhClustersDZetaInterpAcc->Write();
1665 fhDPhiVsDThetaInterpAcc->Write();
1666 fhDPhiVsDZetaInterpAcc->Write();
1668 fhetaClustersLay2->Write();
1669 fhphiClustersLay2->Write();
1670 fhClustersInChip->Write();
1671 for (Int_t nhist=0;nhist<80;nhist++){
1672 fhClustersInModuleLay1[nhist]->Write();
1674 for (Int_t nhist=0;nhist<160;nhist++){
1675 fhClustersInModuleLay2[nhist]->Write();
1679 //__________________________________________________________
1680 Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1682 // Saves the histograms into a tree and saves the trees into a file
1683 // Also the histograms from the base class are saved
1685 if (!GetHistOn()) return kFALSE;
1686 if (!strcmp(filename.Data(),"")) {
1687 AliWarning("WriteHistosToFile: null output filename!");
1690 TFile *hFile=new TFile(filename.Data(),option,
1691 "The File containing the histos for SPD efficiency studies with tracklets");
1692 if(!SaveHists()) return kFALSE;
1697 //____________________________________________________________
1698 void AliITSTrackleterSPDEff::BookHistos() {
1700 // This method books addtitional histograms
1701 // w.r.t. those of the base class.
1702 // In particular, the differences of cluster coordinate between the two SPD
1703 // layers are computed in the interpolation phase
1705 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1707 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
1708 fhClustersDPhiAcc->SetDirectory(0);
1709 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
1710 fhClustersDThetaAcc->SetDirectory(0);
1711 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
1712 fhClustersDZetaAcc->SetDirectory(0);
1714 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
1715 fhDPhiVsDZetaAcc->SetDirectory(0);
1716 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
1717 fhDPhiVsDThetaAcc->SetDirectory(0);
1719 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
1720 fhClustersDPhiAll->SetDirectory(0);
1721 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
1722 fhClustersDThetaAll->SetDirectory(0);
1723 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
1724 fhClustersDZetaAll->SetDirectory(0);
1726 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
1727 fhDPhiVsDZetaAll->SetDirectory(0);
1728 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
1729 fhDPhiVsDThetaAll->SetDirectory(0);
1731 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
1732 fhetaTracklets->SetDirectory(0);
1733 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
1734 fhphiTracklets->SetDirectory(0);
1735 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
1736 fhetaClustersLay1->SetDirectory(0);
1737 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
1738 fhphiClustersLay1->SetDirectory(0);
1740 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1741 fhClustersDPhiInterpAcc->SetDirectory(0);
1742 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1743 fhClustersDThetaInterpAcc->SetDirectory(0);
1744 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1745 fhClustersDZetaInterpAcc->SetDirectory(0);
1747 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1748 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1749 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1750 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1752 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1753 fhClustersDPhiInterpAll->SetDirectory(0);
1754 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1755 fhClustersDThetaInterpAll->SetDirectory(0);
1756 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1757 fhClustersDZetaInterpAll->SetDirectory(0);
1759 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1760 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1761 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1762 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1764 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1765 fhetaClustersLay2->SetDirectory(0);
1766 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1767 fhphiClustersLay2->SetDirectory(0);
1768 fhClustersInChip = new TH1F("fhClustersInChip", "ClustersPerChip", 1200, -0.5, 1199.5);
1769 fhClustersInChip->SetDirectory(0);
1770 // each chip is divided 8(z) x 4(y), i.e. in 32 squares, each containing 4 columns and 64 rows.
1771 Float_t bz[160]; const Float_t kconv = 1.0E-04; // converts microns to cm.
1772 for(Int_t i=0;i<160;i++) bz[i] = 425.0; // most are 425 microns except below
1773 bz[ 31] = bz[ 32] = 625.0; // first chip boundry
1774 bz[ 63] = bz[ 64] = 625.0; // first chip boundry
1775 bz[ 95] = bz[ 96] = 625.0; // first chip boundry
1776 bz[127] = bz[128] = 625.0; // first chip boundry
1777 Double_t xbins[41]; // each bin in x (Z loc coordinate) includes 4 columns
1779 Float_t xmn,xmx,zmn,zmx;
1780 if(!fPlaneEffSPD->GetBlockBoundaries(0,xmn,xmx,zmn,zmx)) AliWarning("Could not book histo properly");
1781 xbins[0]=(Double_t)zmn;
1782 for(Int_t i=0;i<40;i++) {
1783 xbins[i+1]=xbins[i] + (bz[4*i]+bz[4*i+1]+bz[4*i+2]+bz[4*i+3])*kconv;
1785 TString histname="ClustersLay1_mod_",aux;
1786 fhClustersInModuleLay1 =new TH2F*[80];
1787 for (Int_t nhist=0;nhist<80;nhist++){
1791 fhClustersInModuleLay1[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1792 fhClustersInModuleLay1[nhist]->SetName(aux.Data());
1793 fhClustersInModuleLay1[nhist]->SetTitle(aux.Data());
1794 fhClustersInModuleLay1[nhist]->SetDirectory(0);
1796 histname="ClustersLay2_mod_";
1797 fhClustersInModuleLay2 =new TH2F*[160];
1798 for (Int_t nhist=0;nhist<160;nhist++){
1801 fhClustersInModuleLay2[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1802 fhClustersInModuleLay2[nhist]->SetName(aux.Data());
1803 fhClustersInModuleLay2[nhist]->SetTitle(aux.Data());
1804 fhClustersInModuleLay2[nhist]->SetDirectory(0);
1809 //____________________________________________________________
1810 void AliITSTrackleterSPDEff::DeleteHistos() {
1812 // Private method to delete Histograms from memory
1813 // it is called. e.g., by the destructor.
1815 // form AliITSMultReconstructor
1816 if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
1817 if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
1818 if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
1819 if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
1820 if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
1821 if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
1822 if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
1823 if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
1824 if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
1825 if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
1826 if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
1827 if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
1828 if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
1829 if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
1831 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1832 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1833 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1834 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1835 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1836 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1837 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1838 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1839 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1840 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1841 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1842 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1843 if(fhClustersInChip) {delete fhClustersInChip; fhClustersInChip=0;}
1844 if(fhClustersInModuleLay1) {
1845 for (Int_t i=0; i<80; i++ ) delete fhClustersInModuleLay1[i];
1846 delete [] fhClustersInModuleLay1; fhClustersInModuleLay1=0;
1848 if(fhClustersInModuleLay2) {
1849 for (Int_t i=0; i<160; i++ ) delete fhClustersInModuleLay2[i];
1850 delete [] fhClustersInModuleLay2; fhClustersInModuleLay2=0;
1853 //_______________________________________________________________
1854 Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1855 Float_t* vtx, AliStack *stack, TTree *ref) {
1856 // This (private) method can be used only for MC events, where both AliStack and the TrackReference
1858 // It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1860 // - Int_t layer (either 0 or 1): layer which you want to check if the tracklete can be
1862 // - Int_t iC : cluster index used to build the tracklet prediction
1863 // if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1864 // - Float_t* vtx: actual event vertex
1865 // - stack: pointer to Stack
1866 // - ref: pointer to TTRee of TrackReference
1867 Bool_t ret=kFALSE; // returned value
1868 Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1869 if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1870 if(!stack) {AliError("null pointer to MC stack"); return ret;}
1871 if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1872 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1873 if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1875 AliTrackReference *tref=0x0;
1876 Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1877 Int_t nentries = (Int_t)ref->GetEntries();
1878 TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1879 TBranch *br = ref->GetBranch("TrackReferences");
1880 br->SetAddress(&tcaRef);
1881 for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1882 br->GetEntry(itrack);
1883 Int_t nref=tcaRef->GetEntriesFast();
1884 if(nref>0) { //it is enough to look at the first one
1885 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1886 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1889 if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1890 br->GetEntry(imatch); // redundant, nevertheless ...
1891 Int_t nref=tcaRef->GetEntriesFast();
1892 for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1893 tref=(AliTrackReference*)tcaRef->At(iref);
1894 if(tref->R()>10) continue; // not SPD ref
1895 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1896 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1898 // compute the proper quantities for this tref, as was done for fClustersLay1/2
1899 Float_t x = tref->X() - vtx[0];
1900 Float_t y = tref->Y() - vtx[1];
1901 Float_t z = tref->Z() - vtx[2];
1903 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1905 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1906 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1907 trefLayExtr[2] = z; // Store z
1909 if(layer==1) { // try to see if it is reconstructable at the outer layer
1910 // find the difference in angles
1911 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1912 // take into account boundary condition
1913 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1915 // find the difference in z (between linear projection from layer 1
1916 // and the actual point: Dzeta= z1/r1*r2 -z2)
1917 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1918 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1920 // make "elliptical" cut in Phi and Zeta!
1921 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
1922 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
1923 if (d<1) {ret=kTRUE; break;}
1925 if(layer==0) { // try to see if it is reconstructable at the inner layer
1927 // find the difference in angles
1928 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1929 // take into account boundary condition
1930 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1932 // find the difference in z (between linear projection from layer 2
1933 // and the actual point: Dzeta= z2/r2*r1 -z1)
1934 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1935 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1937 // make "elliptical" cut in Phi and Zeta!
1938 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1939 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1940 if (d<1) {ret=kTRUE; break;};
1946 //_________________________________________________________________________
1947 void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
1949 // this method apply a rotation by 180 degree around the Z (beam) axis to all
1950 // the RecPoints in a given layer to be used to build tracklets.
1951 // **************** VERY IMPORTANT:: ***************
1952 // It must be called just after LoadClusterArrays, since afterwards the datamember
1953 // fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
1954 // instead of Cartesian
1956 if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1957 AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1959 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1960 fClustersLay1[iC1][0]*=-1;
1961 fClustersLay1[iC1][1]*=-1;
1965 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1966 fClustersLay2[iC2][0]*=-1;
1967 fClustersLay2[iC2][1]*=-1;
1972 //____________________________________________________________________________
1973 Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *esd){
1974 // This method is used to find the tracklets.
1975 // It is called from AliReconstruction
1976 // The vertex is supposed to be associated to the Tracker (i.e. to this) already
1977 // The cluster is supposed to be associated to the Tracker already
1978 // In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
1981 // apply cuts on the vertex quality
1982 const AliESDVertex *vertex = esd->GetVertex();
1983 if(vertex->GetNContributors()<fMinContVtx) return 0;
1985 AliRunLoader* runLoader = AliRunLoader::Instance();
1987 Error("Clusters2Tracks", "no run loader found");
1990 AliStack *pStack=0x0; TTree *tRefTree=0x0;
1992 runLoader->LoadKinematics("read");
1993 runLoader->LoadTrackRefs("read");
1994 pStack= runLoader->Stack();
1995 tRefTree= runLoader->TreeTR();
1997 Reconstruct(pStack,tRefTree);
1999 if (GetLightBkgStudyInParallel()) {
2000 AliStack *dummy1=0x0; TTree *dummy2=0x0;
2001 ReflectClusterAroundZAxisForLayer(1);
2002 Reconstruct(dummy1,dummy2,kTRUE);
2006 //____________________________________________________________________________
2007 Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
2009 // It is called from AliReconstruction
2015 if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
2016 if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
2017 if(GetLightBkgStudyInParallel()) {
2018 TString name="AliITSPlaneEffSPDtrackletBkg.root";
2019 TFile* pefile = TFile::Open(name, "RECREATE");
2020 rc*=fPlaneEffBkg->Write();
2025 //____________________________________________________________________
2027 AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
2029 // - gets the clusters from the cluster tree
2030 // - convert them into global coordinates
2031 // - store them in the internal arrays
2032 // - count the number of cluster-fired chips
2034 //AliDebug(1,"Loading clusters and cluster-fired chips ...");
2039 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
2040 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
2042 itsClusterBranch->SetAddress(&itsClusters);
2044 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
2045 Float_t cluGlo[3]={0.,0.,0.};
2047 // loop over the its subdetectors
2048 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
2050 if (!itsClusterTree->GetEvent(iIts))
2053 Int_t nClusters = itsClusters->GetEntriesFast();
2055 // number of clusters in each chip of the current module
2058 // loop over clusters
2059 while(nClusters--) {
2060 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
2062 layer = cluster->GetLayer();
2063 if (layer>1) continue;
2065 cluster->GetGlobalXYZ(cluGlo);
2066 Float_t x = cluGlo[0];
2067 Float_t y = cluGlo[1];
2068 Float_t z = cluGlo[2];
2071 fClustersLay1[fNClustersLay1][0] = x;
2072 fClustersLay1[fNClustersLay1][1] = y;
2073 fClustersLay1[fNClustersLay1][2] = z;
2075 for (Int_t i=0; i<3; i++)
2076 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
2079 Int_t det=cluster->GetDetectorIndex();
2080 if(det<0 || det>79) {AliError("Cluster with det. index out of boundaries"); return;}
2081 fhClustersInModuleLay1[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2085 fClustersLay2[fNClustersLay2][0] = x;
2086 fClustersLay2[fNClustersLay2][1] = y;
2087 fClustersLay2[fNClustersLay2][2] = z;
2089 for (Int_t i=0; i<3; i++)
2090 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
2093 Int_t det=cluster->GetDetectorIndex();
2094 if(det<0 || det>159) {AliError("Cluster with det. index out of boundaries"); return;}
2095 fhClustersInModuleLay2[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2099 }// end of cluster loop
2101 } // end of its "subdetector" loop
2103 itsClusters->Delete();
2107 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
2109 //_________________________________________________________________________
2111 AliITSTrackleterSPDEff::SetLightBkgStudyInParallel(Bool_t b) {
2113 // - set Bool_t fLightBackgroundStudyInParallel = b
2114 // a) if you set this kTRUE, then the estimation of the
2115 // SPD efficiency is done as usual for data, but in
2116 // parallel a light (i.e. without control histograms, etc.)
2117 // evaluation of combinatorial background is performed
2118 // with the usual ReflectClusterAroundZAxisForLayer method.
2119 // b) if you set this kFALSE, then you would not have a second
2120 // container for PlaneEfficiency statistics to be used for background
2121 // (fPlaneEffBkg=0). If you want to have a full evaluation of the
2122 // background (with all control histograms and additional data
2123 // members referring to the background) then you have to call the
2124 // method SetReflectClusterAroundZAxisForLayer(kTRUE) esplicitily
2125 fLightBkgStudyInParallel=b;
2126 if(fLightBkgStudyInParallel) {
2127 if(!fPlaneEffBkg) fPlaneEffBkg = new AliITSPlaneEffSPD();
2130 delete fPlaneEffBkg;