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 //____________________________________________________________________
57 ClassImp(AliITSTrackleterSPDEff)
60 //____________________________________________________________________
61 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
71 fOnlyOneTrackletPerC2(0),
78 fhClustersDThetaAcc(0),
79 fhClustersDZetaAcc(0),
81 fhClustersDThetaAll(0),
82 fhClustersDZetaAll(0),
98 fOnlyOneTrackletPerC1(0),
99 fUpdateOncePerEventPlaneEff(0),
101 fChipUpdatedInEvent(0),
104 fReflectClusterAroundZAxisForLayer0(kFALSE),
105 fReflectClusterAroundZAxisForLayer1(kFALSE),
106 fLightBkgStudyInParallel(kFALSE),
108 fUseOnlyPrimaryForPred(0),
109 fUseOnlySecondaryForPred(0),
110 fUseOnlySameParticle(0),
111 fUseOnlyDifferentParticle(0),
112 fUseOnlyStableParticle(1),
113 fPredictionPrimary(0),
114 fPredictionSecondary(0),
116 fClusterSecondary(0),
125 fhClustersDPhiInterpAcc(0),
126 fhClustersDThetaInterpAcc(0),
127 fhClustersDZetaInterpAcc(0),
128 fhClustersDPhiInterpAll(0),
129 fhClustersDThetaInterpAll(0),
130 fhClustersDZetaInterpAll(0),
131 fhDPhiVsDThetaInterpAll(0),
132 fhDPhiVsDThetaInterpAcc(0),
133 fhDPhiVsDZetaInterpAll(0),
134 fhDPhiVsDZetaInterpAcc(0),
135 fhetaClustersLay2(0),
136 fhphiClustersLay2(0),
138 fhClustersInModuleLay1(0),
139 fhClustersInModuleLay2(0)
141 // default constructor
142 // from AliITSMultReconstructor
145 //______________________________________________________________________
146 void AliITSTrackleterSPDEff::Init() {
149 SetOnlyOneTrackletPerC2();
150 fClustersLay1 = new Float_t*[300000];
151 fClustersLay2 = new Float_t*[300000];
152 fTracklets = new Float_t*[300000];
153 fAssociationFlag = new Bool_t[300000];
157 SetOnlyOneTrackletPerC1();
159 fAssociationFlag1 = new Bool_t[300000];
160 fChipPredOnLay2 = new UInt_t[300000];
161 fChipPredOnLay1 = new UInt_t[300000];
162 fChipUpdatedInEvent = new Bool_t[1200];
164 for(Int_t i=0; i<300000; i++) {
165 // from AliITSMultReconstructor
166 fClustersLay1[i] = new Float_t[6];
167 fClustersLay2[i] = new Float_t[6];
168 fTracklets[i] = new Float_t[5];
169 fAssociationFlag[i] = kFALSE;
171 fAssociationFlag1[i] = kFALSE;
173 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
175 if (GetHistOn()) BookHistos();
177 fPlaneEffSPD = new AliITSPlaneEffSPD();
178 SetLightBkgStudyInParallel();
180 //______________________________________________________________________
181 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
183 // from AliITSMultReconstructor
184 fClustersLay1(mr.fClustersLay1),
185 fClustersLay2(mr.fClustersLay2),
186 fTracklets(mr.fTracklets),
187 fAssociationFlag(mr.fAssociationFlag),
188 fNClustersLay1(mr.fNClustersLay1),
189 fNClustersLay2(mr.fNClustersLay2),
190 fNTracklets(mr.fNTracklets),
191 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
192 fPhiWindowL2(mr.fPhiWindowL2),
193 fZetaWindowL2(mr.fZetaWindowL2),
194 fPhiOverlapCut(mr.fPhiOverlapCut),
195 fZetaOverlapCut(mr.fZetaOverlapCut),
197 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
198 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
199 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
200 fhClustersDPhiAll(mr.fhClustersDPhiAll),
201 fhClustersDThetaAll(mr.fhClustersDThetaAll),
202 fhClustersDZetaAll(mr.fhClustersDZetaAll),
203 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
204 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
205 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
206 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
207 fhetaTracklets(mr.fhetaTracklets),
208 fhphiTracklets(mr.fhphiTracklets),
209 fhetaClustersLay1(mr.fhetaClustersLay1),
210 fhphiClustersLay1(mr.fhphiClustersLay1),
212 fAssociationFlag1(mr.fAssociationFlag1),
213 fChipPredOnLay2(mr.fChipPredOnLay2),
214 fChipPredOnLay1(mr.fChipPredOnLay1),
215 fNTracklets1(mr.fNTracklets1),
216 fPhiWindowL1(mr.fPhiWindowL1),
217 fZetaWindowL1(mr.fZetaWindowL1),
218 fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
219 fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
220 fMinContVtx(mr.fMinContVtx),
221 fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
222 fPlaneEffSPD(mr.fPlaneEffSPD),
223 fPlaneEffBkg(mr.fPlaneEffBkg),
224 fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
225 fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
226 fLightBkgStudyInParallel(mr.fLightBkgStudyInParallel),
228 fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
229 fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
230 fUseOnlySameParticle(mr.fUseOnlySameParticle),
231 fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
232 fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
233 fPredictionPrimary(mr.fPredictionPrimary),
234 fPredictionSecondary(mr.fPredictionSecondary),
235 fClusterPrimary(mr.fClusterPrimary),
236 fClusterSecondary(mr.fClusterSecondary),
237 fSuccessPP(mr.fSuccessPP),
238 fSuccessTT(mr.fSuccessTT),
239 fSuccessS(mr.fSuccessS),
240 fSuccessP(mr.fSuccessP),
241 fFailureS(mr.fFailureS),
242 fFailureP(mr.fFailureP),
244 fNonRecons(mr.fNonRecons),
245 fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
246 fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
247 fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
248 fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
249 fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
250 fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
251 fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
252 fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
253 fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
254 fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
255 fhetaClustersLay2(mr.fhetaClustersLay2),
256 fhphiClustersLay2(mr.fhphiClustersLay2),
257 fhClustersInChip(mr.fhClustersInChip),
258 fhClustersInModuleLay1(mr.fhClustersInModuleLay1),
259 fhClustersInModuleLay2(mr.fhClustersInModuleLay2)
264 //______________________________________________________________________
265 AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
266 // Assignment operator
267 this->~AliITSTrackleterSPDEff();
268 new(this) AliITSTrackleterSPDEff(mr);
271 //______________________________________________________________________
272 AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
274 // from AliITSMultReconstructor
276 for(Int_t i=0; i<300000; i++) {
277 delete [] fClustersLay1[i];
278 delete [] fClustersLay2[i];
279 delete [] fTracklets[i];
281 delete [] fClustersLay1;
282 delete [] fClustersLay2;
283 delete [] fTracklets;
284 delete [] fAssociationFlag;
289 delete [] fAssociationFlag1;
291 delete [] fChipPredOnLay2;
292 delete [] fChipPredOnLay1;
294 delete [] fChipUpdatedInEvent;
296 delete [] fPredictionPrimary;
297 delete [] fPredictionSecondary;
298 delete [] fClusterPrimary;
299 delete [] fClusterSecondary;
300 delete [] fSuccessPP;
301 delete [] fSuccessTT;
307 delete [] fNonRecons;
318 //____________________________________________________________________
320 AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef, Bool_t lbkg) {
322 // - you have to take care of the following, before of using Reconstruct
323 // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
324 // and convert the cluster coordinates to theta, phi (seen from the
325 // interaction vertex).
326 // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
327 // - Find the extrapolation/interpolation point.
328 // - Find the chip corresponding to that
329 // - Check if there is a cluster near that point
332 if(lbkg && !GetLightBkgStudyInParallel()) {
333 AliError("You asked for lightBackground in the Reconstruction without proper call to SetLightBkgStudyInParallel(1)");
336 AliITSPlaneEffSPD *pe;
343 // retrieve the vertex position
345 vtx[0]=(Float_t)GetX();
346 vtx[1]=(Float_t)GetY();
347 vtx[2]=(Float_t)GetZ();
348 // to study residual background (i.e. contribution from TT' to measured efficiency)
349 if(fReflectClusterAroundZAxisForLayer0 && !lbkg) ReflectClusterAroundZAxisForLayer(0);
350 if(fReflectClusterAroundZAxisForLayer1 && !lbkg) ReflectClusterAroundZAxisForLayer(1);
352 if(fMC && !pStack && !lbkg) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
353 if(fMC && !tRef && !lbkg) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
355 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
356 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
357 Int_t nfClu1=0; Int_t ntClu1=0;
358 Int_t nfClu2=0; Int_t ntClu2=0;
360 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
361 // for this new event)
362 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
364 // find the tracklets
365 AliDebug(1,"Looking for tracklets... ");
366 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
368 //###########################################################
369 // Loop on layer 1 : finding theta, phi and z
371 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
372 Float_t x = fClustersLay1[iC1][0] - vtx[0];
373 Float_t y = fClustersLay1[iC1][1] - vtx[1];
374 Float_t z = fClustersLay1[iC1][2] - vtx[2];
376 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
378 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
379 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
380 fClustersLay1[iC1][2] = z; // Store z
382 // find the Radius and the chip corresponding to the extrapolation point
384 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
386 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
387 key=999999; // also some other actions should be taken if not Found
389 nfTraPred2+=(Int_t)found; // this for debugging purpose
390 ntTraPred2++; // to check efficiency of the method FindChip
391 fChipPredOnLay2[iC1] = key;
392 fAssociationFlag1[iC1] = kFALSE;
394 if (fHistOn && !lbkg) {
395 Float_t eta=fClustersLay1[iC1][0];
396 eta= TMath::Tan(eta/2.);
397 eta=-TMath::Log(eta);
398 fhetaClustersLay1->Fill(eta);
399 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
400 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
403 // Loop on layer 2 : finding theta, phi and r
404 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
405 Float_t x = fClustersLay2[iC2][0] - vtx[0];
406 Float_t y = fClustersLay2[iC2][1] - vtx[1];
407 Float_t z = fClustersLay2[iC2][2] - vtx[2];
409 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
411 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
412 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
413 fClustersLay2[iC2][2] = z; // Store z
415 // find the Radius and the chip corresponding to the extrapolation point
417 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
419 AliDebug(1,Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
422 nfTraPred1+=(Int_t)found; // this for debugging purpose
423 ntTraPred1++; // to check efficiency of the method FindChip
424 fChipPredOnLay1[iC2] = key;
425 fAssociationFlag[iC2] = kFALSE;
427 if (fHistOn && !lbkg) {
428 Float_t eta=fClustersLay2[iC2][0];
429 eta= TMath::Tan(eta/2.);
430 eta=-TMath::Log(eta);
431 fhetaClustersLay2->Fill(eta);
432 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
433 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
437 //###########################################################
439 // First part : Extrapolation to Layer 2
442 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
444 // here the control to check whether the efficiency of the chip traversed by this tracklet
445 // prediction has already been updated in this event using another tracklet prediction
446 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
448 // reset of variables for multiple candidates
449 Int_t iC2WithBestDist = 0; // reset
450 Float_t distmin = 100.; // just to put a huge number!
451 Float_t dPhimin = 0.; // Used for histograms only!
452 Float_t dThetamin = 0.; // Used for histograms only!
453 Float_t dZetamin = 0.; // Used for histograms only!
455 // in any case, if MC has been required, store statistics of primaries and secondaries
456 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
458 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
459 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
460 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
461 // do it always as a function of the chip number used to built the prediction
462 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
463 if (!found) {AliDebug(1,
464 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
466 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
467 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
468 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
469 { // this cluster is from a primary particle
470 fClusterPrimary[key]++;
472 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
473 } else { // this cluster is from a secondary particle
474 fClusterSecondary[key]++;
476 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
479 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
480 // (in case the prediction is reliable)
481 if( fChipPredOnLay2[iC1]<1200) {
482 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
483 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
484 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
485 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
486 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
487 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
488 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
489 else fNonRecons[fChipPredOnLay2[iC1]]++;
494 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
496 // The following excludes double associations
497 if (!fAssociationFlag[iC2]) {
499 // find the difference in angles
500 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
501 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
502 // take into account boundary condition
503 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
505 // find the difference in z (between linear projection from layer 1
506 // and the actual point: Dzeta= z1/r1*r2 -z2)
507 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
508 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
510 if (fHistOn && !lbkg) {
511 fhClustersDPhiAll->Fill(dPhi);
512 fhClustersDThetaAll->Fill(dTheta);
513 fhClustersDZetaAll->Fill(dZeta);
514 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
515 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
518 // make "elliptical" cut in Phi and Zeta!
519 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
520 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
524 //look for the minimum distance: the minimum is in iC2WithBestDist
525 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
526 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
530 iC2WithBestDist = iC2;
533 } // end of loop over clusters in layer 2
535 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
537 if (fHistOn && !lbkg) {
538 fhClustersDPhiAcc->Fill(dPhimin);
539 fhClustersDThetaAcc->Fill(dThetamin);
540 fhClustersDZetaAcc->Fill(dZetamin);
541 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
542 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
545 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
546 // flag the association
548 // store the tracklet
550 // use the theta from the clusters in the first layer
551 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
552 // use the phi from the clusters in the first layer
553 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
554 // Store the difference between phi1 and phi2
555 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
562 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
574 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
578 fTracklets[fNTracklets][3] = -2;
581 if (fHistOn && !lbkg) {
582 Float_t eta=fTracklets[fNTracklets][0];
583 eta= TMath::Tan(eta/2.);
584 eta=-TMath::Log(eta);
585 fhetaTracklets->Fill(eta);
586 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
589 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
590 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
593 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
596 nfClu2+=(Int_t)found; // this for debugging purpose
597 ntClu2++; // to check efficiency of the method FindChip
598 if(key<1200) { // the Chip has been found
599 if(fMC && !lbkg) { // this part only for MC
600 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
601 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
602 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
605 if(primary) fSuccessPP[key]++;
607 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
608 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
611 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
613 pe->UpDatePlaneEff(kTRUE,key); // success
614 fChipUpdatedInEvent[key]=kTRUE;
616 if(primary) fSuccessP[key]++;
617 if(secondary) fSuccessS[key]++;
621 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
622 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
624 if(primary) fSuccessP[key]++;
625 if(secondary) fSuccessS[key]++;
632 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
633 else if (fChipPredOnLay2[iC1]<1200) {
634 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
635 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
637 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
638 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
641 } // end of loop over clusters in layer 1
643 fNTracklets1=fNTracklets;
645 //###################################################################
647 // Second part : Interpolation to Layer 1
650 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
652 // here the control to check whether the efficiency of the chip traversed by this tracklet
653 // prediction has already been updated in this event using another tracklet prediction
654 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
656 // reset of variables for multiple candidates
657 Int_t iC1WithBestDist = 0; // reset
658 Float_t distmin = 100.; // just to put a huge number!
659 Float_t dPhimin = 0.; // Used for histograms only!
660 Float_t dThetamin = 0.; // Used for histograms only!
661 Float_t dZetamin = 0.; // Used for histograms only!
663 // in any case, if MC has been required, store statistics of primaries and secondaries
664 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
666 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
667 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
668 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
669 // do it always as a function of the chip number used to built the prediction
670 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
671 if (!found) {AliDebug(1,
672 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
674 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
675 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
676 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
677 { // this cluster is from a primary particle
678 fClusterPrimary[key]++;
680 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
681 } else { // this cluster is from a secondary particle
682 fClusterSecondary[key]++;
684 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
687 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
688 // (in case the prediction is reliable)
689 if( fChipPredOnLay1[iC2]<1200) {
690 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
691 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
692 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
693 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
694 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
695 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
696 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
697 else fNonRecons[fChipPredOnLay1[iC2]]++;
702 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
704 // The following excludes double associations
705 if (!fAssociationFlag1[iC1]) {
707 // find the difference in angles
708 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
709 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
710 // take into account boundary condition
711 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
714 // find the difference in z (between linear projection from layer 2
715 // and the actual point: Dzeta= z2/r2*r1 -z1)
716 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
717 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
720 if (fHistOn && !lbkg) {
721 fhClustersDPhiInterpAll->Fill(dPhi);
722 fhClustersDThetaInterpAll->Fill(dTheta);
723 fhClustersDZetaInterpAll->Fill(dZeta);
724 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
725 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
727 // make "elliptical" cut in Phi and Zeta!
728 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
729 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
733 //look for the minimum distance: the minimum is in iC1WithBestDist
734 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
735 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
739 iC1WithBestDist = iC1;
742 } // end of loop over clusters in layer 1
744 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
746 if (fHistOn && !lbkg) {
747 fhClustersDPhiInterpAcc->Fill(dPhimin);
748 fhClustersDThetaInterpAcc->Fill(dThetamin);
749 fhClustersDZetaInterpAcc->Fill(dZetamin);
750 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
751 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
754 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
755 // flag the association
757 // store the tracklet
759 // use the theta from the clusters in the first layer
760 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
761 // use the phi from the clusters in the first layer
762 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
763 // Store the difference between phi1 and phi2
764 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
771 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
783 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
787 fTracklets[fNTracklets][3] = -2;
790 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
791 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
794 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
797 nfClu1+=(Int_t)found; // this for debugging purpose
798 ntClu1++; // to check efficiency of the method FindChip
800 if(fMC && !lbkg) { // this part only for MC
801 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
802 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
803 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
804 if (label2 < 3) { // same label
806 if(primary) fSuccessPP[key]++;
808 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
809 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
812 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
814 pe->UpDatePlaneEff(kTRUE,key); // success
815 fChipUpdatedInEvent[key]=kTRUE;
817 if(primary) fSuccessP[key]++;
818 if(secondary) fSuccessS[key]++;
821 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
822 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
824 if(primary) fSuccessP[key]++;
825 if(secondary) fSuccessS[key]++;
832 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
833 else if (fChipPredOnLay1[iC2]<1200) {
834 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
835 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
837 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
838 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
841 } // end of loop over clusters in layer 2
843 AliDebug(1,Form("%d tracklets found", fNTracklets));
844 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
845 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
846 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
847 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
849 //____________________________________________________________________
850 Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer,const Float_t* vtx,
851 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
853 // Input: a) layer number in the range [0,1]
854 // b) vtx[3]: actual vertex
855 // c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
856 // d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
858 // Output: Unique key to locate a chip
859 // return: kTRUE if succesfull
861 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
862 Double_t r=GetRLayer(layer);
863 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
865 // set phiVtx in the range [0,2pi]
866 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
868 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
869 // of the intersection of the tracklet with the pixel layer.
870 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
872 if(TMath::Abs(thetaVtx)<1E-6) return kFALSE;
873 zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
875 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
876 Double_t vtxy[2]={vtx[0],vtx[1]};
877 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
878 // this method gives you two interceptions
879 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
880 // set phiAbs in the range [0,2pi]
881 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
882 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
883 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
884 // taking the closest one to phiVtx
885 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
886 } else phiAbs=phiVtx;
887 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
889 // now you need to locate the chip within the idet detector,
890 // starting from the local coordinates in such a detector
892 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
893 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
894 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
896 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
899 //______________________________________________________________________________
900 Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
902 // Return the average radius of a layer from Geometry
904 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
905 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
907 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
908 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
909 Double_t r=TMath::Sqrt(x*x + y*y);
911 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
912 r += TMath::Sqrt(x*x + y*y);
913 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
914 r += TMath::Sqrt(x*x + y*y);
915 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
916 r += TMath::Sqrt(x*x + y*y);
920 //______________________________________________________________________________
921 Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
922 Float_t &xloc, Float_t &zloc) {
923 // this method transform Global Cilindrical coordinates into local (i.e. module)
924 // cartesian coordinates
926 //Compute Cartesian Global Coordinate
927 Double_t xyzGlob[3],xyzLoc[3];
929 xyzGlob[0]=r*TMath::Cos(phi);
930 xyzGlob[1]=r*TMath::Sin(phi);
935 if(idet<0) return kFALSE;
937 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
938 Int_t lad = Int_t(idet/ndet) + 1;
939 Int_t det = idet - (lad-1)*ndet + 1;
941 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
943 xloc = (Float_t)xyzLoc[0];
944 zloc = (Float_t)xyzLoc[2];
948 //______________________________________________________________________________
949 Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
950 //--------------------------------------------------------------------
951 // This function finds the detector crossed by the track
952 // Input: layer in range [0,1]
953 // phi in ALICE absolute reference system
955 //--------------------------------------------------------------------
956 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
957 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
958 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
959 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
961 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
962 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
963 Double_t phiOffset=TMath::ATan2(y,x);
967 if (zOffset<0) // old geometry
968 dphi = -(phi-phiOffset);
970 dphi = phi-phiOffset;
972 if (dphi < 0) dphi += 2*TMath::Pi();
973 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
974 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
975 if (np>=nladders) np-=nladders;
976 if (np<0) np+=nladders;
978 Double_t dz=zOffset-z;
979 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
980 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
981 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
982 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
984 return np*ndetectors + nz;
986 //____________________________________________________________
987 Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
988 // this method find the intersection in xy between a tracklet (straight line) and
989 // a circonference (r=R), using polar coordinates.
991 Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
992 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
993 - R: radius of the circle
994 Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
996 Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
997 r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
998 In the same system, the equation of a semi-line is: phi=phiVtx;
999 Hence you get one interception only: P=(r,phiVtx)
1000 Finally you want P in the ABSOLUTE ALICE reference system.
1002 Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
1003 Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
1004 Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
1005 Double_t cC=rO*rO-R*R;
1006 Double_t dDelta=bB*bB-4*cC;
1007 if(dDelta<0) return kFALSE;
1009 r1=(-bB-TMath::Sqrt(dDelta))/2;
1010 r2=(-bB+TMath::Sqrt(dDelta))/2;
1011 if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
1012 Double_t r=TMath::Max(r1,r2); // take the positive
1013 Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
1014 Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
1015 pvtx[0]=r*TMath::Cos(phiVtx);
1016 pvtx[1]=r*TMath::Sin(phiVtx);
1017 pP[0]=vtx[0]+pvtx[0];
1018 pP[1]=vtx[1]+pvtx[1];
1019 phi=TMath::ATan2(pP[1],pP[0]);
1022 //___________________________________________________________
1023 Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) const {
1025 // simple method to reduce all angles (in rad)
1029 while(angle >=2*TMath::Pi() || angle<0) {
1030 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
1031 if(angle < 0) angle+=2*TMath::Pi();
1035 //___________________________________________________________
1036 Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
1038 // This method check if a particle is primary; i.e.
1039 // it comes from the main vertex and it is a "stable" particle, according to
1040 // AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
1041 // a stable particle: it has no effect on this analysis).
1042 // This method can be called only for MC events, where Kinematics is available.
1043 // if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
1044 // returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
1045 // The latter (see below) try to verify if a primary particle is also "detectable".
1047 if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
1048 if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
1049 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
1050 // return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
1051 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
1052 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
1053 TParticle* part = stack->Particle(ipart);
1054 TParticle* part0 = stack->Particle(0); // first primary
1055 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
1056 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
1057 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
1059 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
1060 TParticlePDG* pdgPart = part->GetPDG();
1061 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
1063 Double_t distx = part->Vx() - part0->Vx();
1064 Double_t disty = part->Vy() - part0->Vy();
1065 Double_t distz = part->Vz() - part0->Vz();
1066 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
1068 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
1069 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
1070 return kFALSE; }// primary if within 500 microns from true Vertex
1072 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
1075 //_____________________________________________________________________________________________
1076 Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
1078 // This private method can be applied on MC particles (if stack is available),
1079 // provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
1081 // It define "detectable" a primary particle according to the following criteria:
1083 // - if no decay products can be found in the stack (note that this does not
1084 // means it is stable, since a particle is stored in stack if it has at least 1 hit in a
1085 // sensitive detector)
1086 // - if it has at least one decay daughter produced outside or just on the outer pixel layer
1087 // - if the last decay particle is an electron (or a muon) which is not produced in-between
1088 // the two pixel layers (this is likely to be a kink).
1089 if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
1090 if(!stack) {AliError("null pointer to MC stack"); return 0;}
1091 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
1093 TParticle* part = stack->Particle(ipart);
1094 //TParticle* part0 = stack->Particle(0); // first primary
1100 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
1101 // its real fate ! But you have to take it !
1102 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
1103 Int_t lastDau = part->GetLastDaughter();
1104 nDau = lastDau - firstDau + 1;
1105 Double_t distMax=0.;
1107 for(Int_t j=firstDau; j<=lastDau; j++) {
1108 dau = stack->Particle(j);
1109 Double_t distx = dau->Vx();
1110 Double_t disty = dau->Vy();
1111 //Double_t distz = dau->Vz();
1112 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
1113 if(distR<distMax) continue; // considere only the daughter produced at largest radius
1117 dau = stack->Particle(jmax);
1118 pdgDau=dau->GetPdgCode();
1119 if (pdgDau == 11 || pdgDau == 13 ) {
1120 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
1121 else nret =0; // delta-ray emission in material (keep it)
1123 else {// not ele or muon
1124 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
1128 //_________________________________________________________________
1129 void AliITSTrackleterSPDEff::InitPredictionMC() {
1131 // this method allocate memory for the MC related informations
1132 // all the counters are set to 0
1135 if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1136 fPredictionPrimary = new Int_t[1200];
1137 fPredictionSecondary = new Int_t[1200];
1138 fClusterPrimary = new Int_t[1200];
1139 fClusterSecondary = new Int_t[1200];
1140 fSuccessPP = new Int_t[1200];
1141 fSuccessTT = new Int_t[1200];
1142 fSuccessS = new Int_t[1200];
1143 fSuccessP = new Int_t[1200];
1144 fFailureS = new Int_t[1200];
1145 fFailureP = new Int_t[1200];
1146 fRecons = new Int_t[1200];
1147 fNonRecons = new Int_t[1200];
1148 for(Int_t i=0; i<1200; i++) {
1149 fPredictionPrimary[i]=0;
1150 fPredictionSecondary[i]=0;
1151 fPredictionSecondary[i]=0;
1152 fClusterSecondary[i]=0;
1164 //_________________________________________________________________
1165 void AliITSTrackleterSPDEff::DeletePredictionMC() {
1167 // this method deallocate memory for the MC related informations
1168 // all the counters are set to 0
1171 if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1172 if(fPredictionPrimary) {
1173 delete fPredictionPrimary; fPredictionPrimary=0;
1175 if(fPredictionSecondary) {
1176 delete fPredictionSecondary; fPredictionSecondary=0;
1178 if(fClusterPrimary) {
1179 delete fClusterPrimary; fClusterPrimary=0;
1181 if(fClusterSecondary) {
1182 delete fClusterSecondary; fClusterSecondary=0;
1185 delete fSuccessPP; fSuccessPP=0;
1188 delete fSuccessTT; fSuccessTT=0;
1191 delete fSuccessS; fSuccessS=0;
1194 delete fSuccessP; fSuccessP=0;
1197 delete fFailureS; fFailureS=0;
1200 delete fFailureP; fFailureP=0;
1203 delete fRecons; fRecons=0;
1206 delete fNonRecons; fNonRecons=0;
1210 //______________________________________________________________________
1211 Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
1213 // This method return the Data menmber fPredictionPrimary [1200].
1214 // You can call it only for MC events.
1215 // fPredictionPrimary[key] contains the number of tracklet predictions on the
1216 // given chip key built using a cluster on the other layer produced (at least)
1217 // from a primary particle.
1218 // Key refers to the chip crossed by the prediction
1221 if (!fMC) {CallWarningMC(); return 0;}
1222 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1223 return fPredictionPrimary[(Int_t)key];
1225 //______________________________________________________________________
1226 Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
1228 // This method return the Data menmber fPredictionSecondary [1200].
1229 // You can call it only for MC events.
1230 // fPredictionSecondary[key] contains the number of tracklet predictions on the
1231 // given chip key built using a cluster on the other layer produced (only)
1232 // from a secondary particle
1233 // Key refers to the chip crossed by the prediction
1236 if (!fMC) {CallWarningMC(); return 0;}
1237 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1238 return fPredictionSecondary[(Int_t)key];
1240 //______________________________________________________________________
1241 Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
1243 // This method return the Data menmber fClusterPrimary [1200].
1244 // You can call it only for MC events.
1245 // fClusterPrimary[key] contains the number of tracklet predictions
1246 // built using a cluster on that layer produced (only)
1247 // from a primary particle
1248 // Key refers to the chip used to build the prediction
1251 if (!fMC) {CallWarningMC(); return 0;}
1252 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1253 return fClusterPrimary[(Int_t)key];
1255 //______________________________________________________________________
1256 Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
1258 // This method return the Data menmber fClusterSecondary [1200].
1259 // You can call it only for MC events.
1260 // fClusterSecondary[key] contains the number of tracklet predictions
1261 // built using a cluster on that layer produced (only)
1262 // from a secondary particle
1263 // Key refers to the chip used to build the prediction
1265 if (!fMC) {CallWarningMC(); return 0;}
1266 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1267 return fClusterSecondary[(Int_t)key];
1269 //______________________________________________________________________
1270 Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1272 // This method return the Data menmber fSuccessPP [1200].
1273 // You can call it only for MC events.
1274 // fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1275 // with a cluster on the other layer) built by using the same primary particle
1276 // the unique chip key refers to the chip which get updated its efficiency
1278 if (!fMC) {CallWarningMC(); return 0;}
1279 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1280 return fSuccessPP[(Int_t)key];
1282 //______________________________________________________________________
1283 Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1285 // This method return the Data menmber fSuccessTT [1200].
1286 // You can call it only for MC events.
1287 // fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1288 // with a cluster on the other layer) built by using the same particle (whatever)
1289 // the unique chip key refers to the chip which get updated its efficiency
1291 if (!fMC) {CallWarningMC(); return 0;}
1292 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1293 return fSuccessTT[(Int_t)key];
1295 //______________________________________________________________________
1296 Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1298 // This method return the Data menmber fSuccessS [1200].
1299 // You can call it only for MC events.
1300 // fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1301 // with a cluster on the other layer) built by using a secondary particle
1302 // the unique chip key refers to the chip which get updated its efficiency
1304 if (!fMC) {CallWarningMC(); return 0;}
1305 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1306 return fSuccessS[(Int_t)key];
1308 //______________________________________________________________________
1309 Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1311 // This method return the Data menmber fSuccessP [1200].
1312 // You can call it only for MC events.
1313 // fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1314 // with a cluster on the other layer) built by using a primary particle
1315 // the unique chip key refers to the chip which get updated its efficiency
1317 if (!fMC) {CallWarningMC(); return 0;}
1318 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1319 return fSuccessP[(Int_t)key];
1321 //______________________________________________________________________
1322 Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1324 // This method return the Data menmber fFailureS [1200].
1325 // You can call it only for MC events.
1326 // fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1327 // with a cluster on the other layer) built by using a secondary particle
1328 // the unique chip key refers to the chip which get updated its efficiency
1330 if (!fMC) {CallWarningMC(); return 0;}
1331 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1332 return fFailureS[(Int_t)key];
1334 //______________________________________________________________________
1335 Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1337 // This method return the Data menmber fFailureP [1200].
1338 // You can call it only for MC events.
1339 // fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1340 // with a cluster on the other layer) built by using a primary particle
1341 // the unique chip key refers to the chip which get updated its efficiency
1343 if (!fMC) {CallWarningMC(); return 0;}
1344 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1345 return fFailureP[(Int_t)key];
1347 //_____________________________________________________________________
1348 Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1350 // This method return the Data menmber fRecons [1200].
1351 // You can call it only for MC events.
1352 // fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1353 // has an hit in the detector)
1354 // the unique chip key refers to the chip where fall the prediction
1356 if (!fMC) {CallWarningMC(); return 0;}
1357 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1358 return fRecons[(Int_t)key];
1360 //_____________________________________________________________________
1361 Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1363 // This method return the Data menmber fNonRecons [1200].
1364 // You can call it only for MC events.
1365 // fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1366 // has not any hit in the detector)
1367 // the unique chip key refers to the chip where fall the prediction
1369 if (!fMC) {CallWarningMC(); return 0;}
1370 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1371 return fNonRecons[(Int_t)key];
1373 //______________________________________________________________________
1374 void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1375 // Print out some class data values in Ascii Form to output stream
1377 // ostream *os Output stream where Ascii data is to be writen
1382 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindowL2 <<" "<< fZetaWindowL2
1383 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1384 << " " << fUpdateOncePerEventPlaneEff << " " << fMinContVtx
1385 << " " << fReflectClusterAroundZAxisForLayer0
1386 << " " << fReflectClusterAroundZAxisForLayer1;
1388 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1389 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1390 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1391 << " " << fUseOnlyStableParticle ;
1392 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1393 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1394 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1395 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
1396 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1397 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1398 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1399 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1400 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1401 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1402 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1403 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
1406 //______________________________________________________________________
1407 void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1408 // Read in some class data values in Ascii Form to output stream
1410 // istream *is Input stream where Ascii data is to be read in from
1417 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindowL2 >> fZetaWindowL2
1418 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1419 >> fUpdateOncePerEventPlaneEff >> fMinContVtx
1420 >> fReflectClusterAroundZAxisForLayer0
1421 >> fReflectClusterAroundZAxisForLayer1;
1422 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
1424 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1426 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1427 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1428 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1429 >> fUseOnlyStableParticle;
1430 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1431 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1432 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1433 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1434 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1435 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1436 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1437 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1438 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1439 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1440 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1441 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1445 //______________________________________________________________________
1446 ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1447 // Standard output streaming function
1449 // ostream &os output steam
1450 // AliITSTrackleterSPDEff &s class to be streamed.
1454 // ostream &os The stream pointer
1459 //______________________________________________________________________
1460 istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1461 // Standard inputput streaming function
1463 // istream &is input steam
1464 // AliITSTrackleterSPDEff &s class to be streamed.
1468 // ostream &os The stream pointer
1470 //printf("prova %d \n", (Int_t)s.GetMC());
1474 //______________________________________________________________________
1475 void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
1477 // This Method write into an either asci or root file
1478 // the used cuts and the statistics of the MC related quantities
1479 // The method SetMC() has to be called before
1480 // Input TString filename: name of file for output (it deletes already existing
1485 //if(!fMC) {CallWarningMC(); return;}
1486 if (!filename.Contains(".root")) {
1487 ofstream out(filename.Data(),ios::out | ios::binary);
1493 TFile* mcfile = TFile::Open(filename, "RECREATE");
1494 TH1F* cuts = new TH1F("cuts", "list of cuts", 11, 0, 11); // TH1I containing cuts
1495 cuts->SetBinContent(1,fPhiWindowL1);
1496 cuts->SetBinContent(2,fZetaWindowL1);
1497 cuts->SetBinContent(3,fPhiWindowL2);
1498 cuts->SetBinContent(4,fZetaWindowL2);
1499 cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
1500 cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
1501 cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
1502 cuts->SetBinContent(8,fMinContVtx);
1503 cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer0);
1504 cuts->SetBinContent(10,fReflectClusterAroundZAxisForLayer1);
1505 cuts->SetBinContent(11,fMC);
1508 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
1510 TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
1511 mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
1512 mc0->SetBinContent(2,fUseOnlySecondaryForPred);
1513 mc0->SetBinContent(3,fUseOnlySameParticle);
1514 mc0->SetBinContent(4,fUseOnlyDifferentParticle);
1515 mc0->SetBinContent(5,fUseOnlyStableParticle);
1519 mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
1520 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
1522 mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
1523 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
1525 mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
1526 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
1528 mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
1529 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
1531 mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
1532 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
1534 mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
1535 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
1537 mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
1538 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
1540 mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
1541 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
1543 mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
1544 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
1546 mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
1547 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
1549 mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
1550 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
1552 mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
1553 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
1561 //____________________________________________________________________
1562 void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
1564 // This Method read from an asci file (do not know why binary does not work)
1565 // the cuts to be used and the statistics of the MC related quantities
1566 // Input TString filename: name of input file for output
1567 // The method SetMC() has to be called before
1571 //if(!fMC) {CallWarningMC(); return;}
1572 if( gSystem->AccessPathName( filename.Data() ) ) {
1573 AliError( Form( "file (%s) not found", filename.Data() ) );
1577 if (!filename.Contains(".root")) {
1578 ifstream in(filename.Data(),ios::in | ios::binary);
1585 TFile *mcfile = TFile::Open(filename);
1586 TH1F *cuts = (TH1F*)mcfile->Get("cuts");
1587 fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
1588 fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
1589 fPhiWindowL2=(Float_t)cuts->GetBinContent(3);
1590 fZetaWindowL2=(Float_t)cuts->GetBinContent(4);
1591 fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
1592 fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
1593 fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
1594 fMinContVtx=(Int_t)cuts->GetBinContent(8);
1595 fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(9);
1596 fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(10);
1597 fMC=(Bool_t)cuts->GetBinContent(11);
1598 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1599 else { // only if file with MC predictions
1600 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1601 TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
1602 fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
1603 fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
1604 fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
1605 fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
1606 fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
1608 mc1 =(TH1I*)mcfile->Get("mc1");
1609 for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1610 mc1 =(TH1I*)mcfile->Get("mc2");
1611 for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1612 mc1 =(TH1I*)mcfile->Get("mc3");
1613 for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1614 mc1 =(TH1I*)mcfile->Get("mc4");
1615 for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1616 mc1 =(TH1I*)mcfile->Get("mc5");
1617 for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1618 mc1 =(TH1I*)mcfile->Get("mc6");
1619 for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
1620 mc1 =(TH1I*)mcfile->Get("mc7");
1621 for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1622 mc1 =(TH1I*)mcfile->Get("mc8");
1623 for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1624 mc1 =(TH1I*)mcfile->Get("mc9");
1625 for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1626 mc1 =(TH1I*)mcfile->Get("mc10");
1627 for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1628 mc1 =(TH1I*)mcfile->Get("mc11");
1629 for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1630 mc1 =(TH1I*)mcfile->Get("mc12");
1631 for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1637 //____________________________________________________________________
1638 Bool_t AliITSTrackleterSPDEff::SaveHists() {
1639 // This (private) method save the histograms on the output file
1640 // (only if fHistOn is TRUE).
1641 // Also the histograms from the base class are saved through the
1642 // AliITSMultReconstructor::SaveHists() call
1644 if (!GetHistOn()) return kFALSE;
1646 // AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1647 fhClustersDPhiAll->Write();
1648 fhClustersDThetaAll->Write();
1649 fhClustersDZetaAll->Write();
1650 fhDPhiVsDThetaAll->Write();
1651 fhDPhiVsDZetaAll->Write();
1653 fhClustersDPhiAcc->Write();
1654 fhClustersDThetaAcc->Write();
1655 fhClustersDZetaAcc->Write();
1656 fhDPhiVsDThetaAcc->Write();
1657 fhDPhiVsDZetaAcc->Write();
1659 fhetaTracklets->Write();
1660 fhphiTracklets->Write();
1661 fhetaClustersLay1->Write();
1662 fhphiClustersLay1->Write();
1664 fhClustersDPhiInterpAll->Write();
1665 fhClustersDThetaInterpAll->Write();
1666 fhClustersDZetaInterpAll->Write();
1667 fhDPhiVsDThetaInterpAll->Write();
1668 fhDPhiVsDZetaInterpAll->Write();
1670 fhClustersDPhiInterpAcc->Write();
1671 fhClustersDThetaInterpAcc->Write();
1672 fhClustersDZetaInterpAcc->Write();
1673 fhDPhiVsDThetaInterpAcc->Write();
1674 fhDPhiVsDZetaInterpAcc->Write();
1676 fhetaClustersLay2->Write();
1677 fhphiClustersLay2->Write();
1678 fhClustersInChip->Write();
1679 for (Int_t nhist=0;nhist<80;nhist++){
1680 fhClustersInModuleLay1[nhist]->Write();
1682 for (Int_t nhist=0;nhist<160;nhist++){
1683 fhClustersInModuleLay2[nhist]->Write();
1687 //__________________________________________________________
1688 Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1690 // Saves the histograms into a tree and saves the trees into a file
1691 // Also the histograms from the base class are saved
1693 if (!GetHistOn()) return kFALSE;
1694 if (!strcmp(filename.Data(),"")) {
1695 AliWarning("WriteHistosToFile: null output filename!");
1698 TFile *hFile=new TFile(filename.Data(),option,
1699 "The File containing the histos for SPD efficiency studies with tracklets");
1700 if(!SaveHists()) return kFALSE;
1705 //____________________________________________________________
1706 void AliITSTrackleterSPDEff::BookHistos() {
1708 // This method books addtitional histograms
1709 // w.r.t. those of the base class.
1710 // In particular, the differences of cluster coordinate between the two SPD
1711 // layers are computed in the interpolation phase
1713 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1715 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
1716 fhClustersDPhiAcc->SetDirectory(0);
1717 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
1718 fhClustersDThetaAcc->SetDirectory(0);
1719 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
1720 fhClustersDZetaAcc->SetDirectory(0);
1722 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
1723 fhDPhiVsDZetaAcc->SetDirectory(0);
1724 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
1725 fhDPhiVsDThetaAcc->SetDirectory(0);
1727 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
1728 fhClustersDPhiAll->SetDirectory(0);
1729 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
1730 fhClustersDThetaAll->SetDirectory(0);
1731 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
1732 fhClustersDZetaAll->SetDirectory(0);
1734 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
1735 fhDPhiVsDZetaAll->SetDirectory(0);
1736 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
1737 fhDPhiVsDThetaAll->SetDirectory(0);
1739 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
1740 fhetaTracklets->SetDirectory(0);
1741 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
1742 fhphiTracklets->SetDirectory(0);
1743 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
1744 fhetaClustersLay1->SetDirectory(0);
1745 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
1746 fhphiClustersLay1->SetDirectory(0);
1748 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1749 fhClustersDPhiInterpAcc->SetDirectory(0);
1750 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1751 fhClustersDThetaInterpAcc->SetDirectory(0);
1752 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1753 fhClustersDZetaInterpAcc->SetDirectory(0);
1755 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1756 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1757 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1758 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1760 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1761 fhClustersDPhiInterpAll->SetDirectory(0);
1762 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1763 fhClustersDThetaInterpAll->SetDirectory(0);
1764 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1765 fhClustersDZetaInterpAll->SetDirectory(0);
1767 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1768 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1769 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1770 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1772 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1773 fhetaClustersLay2->SetDirectory(0);
1774 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1775 fhphiClustersLay2->SetDirectory(0);
1776 fhClustersInChip = new TH1F("fhClustersInChip", "ClustersPerChip", 1200, -0.5, 1199.5);
1777 fhClustersInChip->SetDirectory(0);
1778 // each chip is divided 8(z) x 4(y), i.e. in 32 squares, each containing 4 columns and 64 rows.
1779 Float_t bz[160]; const Float_t kconv = 1.0E-04; // converts microns to cm.
1780 for(Int_t i=0;i<160;i++) bz[i] = 425.0; // most are 425 microns except below
1781 bz[ 31] = bz[ 32] = 625.0; // first chip boundry
1782 bz[ 63] = bz[ 64] = 625.0; // first chip boundry
1783 bz[ 95] = bz[ 96] = 625.0; // first chip boundry
1784 bz[127] = bz[128] = 625.0; // first chip boundry
1785 Double_t xbins[41]; // each bin in x (Z loc coordinate) includes 4 columns
1787 Float_t xmn,xmx,zmn,zmx;
1788 if(!fPlaneEffSPD->GetBlockBoundaries(0,xmn,xmx,zmn,zmx)) AliWarning("Could not book histo properly");
1789 xbins[0]=(Double_t)zmn;
1790 for(Int_t i=0;i<40;i++) {
1791 xbins[i+1]=xbins[i] + (bz[4*i]+bz[4*i+1]+bz[4*i+2]+bz[4*i+3])*kconv;
1793 TString histname="ClustersLay1_mod_",aux;
1794 fhClustersInModuleLay1 =new TH2F*[80];
1795 for (Int_t nhist=0;nhist<80;nhist++){
1799 fhClustersInModuleLay1[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1800 fhClustersInModuleLay1[nhist]->SetName(aux.Data());
1801 fhClustersInModuleLay1[nhist]->SetTitle(aux.Data());
1802 fhClustersInModuleLay1[nhist]->SetDirectory(0);
1804 histname="ClustersLay2_mod_";
1805 fhClustersInModuleLay2 =new TH2F*[160];
1806 for (Int_t nhist=0;nhist<160;nhist++){
1809 fhClustersInModuleLay2[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1810 fhClustersInModuleLay2[nhist]->SetName(aux.Data());
1811 fhClustersInModuleLay2[nhist]->SetTitle(aux.Data());
1812 fhClustersInModuleLay2[nhist]->SetDirectory(0);
1817 //____________________________________________________________
1818 void AliITSTrackleterSPDEff::DeleteHistos() {
1820 // Private method to delete Histograms from memory
1821 // it is called. e.g., by the destructor.
1823 // form AliITSMultReconstructor
1824 if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
1825 if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
1826 if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
1827 if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
1828 if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
1829 if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
1830 if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
1831 if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
1832 if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
1833 if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
1834 if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
1835 if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
1836 if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
1837 if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
1839 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1840 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1841 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1842 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1843 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1844 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1845 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1846 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1847 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1848 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1849 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1850 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1851 if(fhClustersInChip) {delete fhClustersInChip; fhClustersInChip=0;}
1852 if(fhClustersInModuleLay1) {
1853 for (Int_t i=0; i<80; i++ ) delete fhClustersInModuleLay1[i];
1854 delete [] fhClustersInModuleLay1; fhClustersInModuleLay1=0;
1856 if(fhClustersInModuleLay2) {
1857 for (Int_t i=0; i<160; i++ ) delete fhClustersInModuleLay2[i];
1858 delete [] fhClustersInModuleLay2; fhClustersInModuleLay2=0;
1861 //_______________________________________________________________
1862 Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1863 const Float_t* vtx, const AliStack *stack, TTree *ref) {
1864 // This (private) method can be used only for MC events, where both AliStack and the TrackReference
1866 // It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1868 // - Int_t layer (either 0 or 1): layer which you want to check if the tracklete can be
1870 // - Int_t iC : cluster index used to build the tracklet prediction
1871 // if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1872 // - Float_t* vtx: actual event vertex
1873 // - stack: pointer to Stack
1874 // - ref: pointer to TTRee of TrackReference
1875 Bool_t ret=kFALSE; // returned value
1876 Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1877 if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1878 if(!stack) {AliError("null pointer to MC stack"); return ret;}
1879 if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1880 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1881 if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1883 AliTrackReference *tref=0x0;
1884 Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1885 Int_t nentries = (Int_t)ref->GetEntries();
1886 TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1887 TBranch *br = ref->GetBranch("TrackReferences");
1888 br->SetAddress(&tcaRef);
1889 for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1890 br->GetEntry(itrack);
1891 Int_t nref=tcaRef->GetEntriesFast();
1892 if(nref>0) { //it is enough to look at the first one
1893 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1894 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1897 if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1898 br->GetEntry(imatch); // redundant, nevertheless ...
1899 Int_t nref=tcaRef->GetEntriesFast();
1900 for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1901 tref=(AliTrackReference*)tcaRef->At(iref);
1902 if(tref->R()>10) continue; // not SPD ref
1903 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1904 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1906 // compute the proper quantities for this tref, as was done for fClustersLay1/2
1907 Float_t x = tref->X() - vtx[0];
1908 Float_t y = tref->Y() - vtx[1];
1909 Float_t z = tref->Z() - vtx[2];
1911 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1913 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1914 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1915 trefLayExtr[2] = z; // Store z
1917 if(layer==1) { // try to see if it is reconstructable at the outer layer
1918 // find the difference in angles
1919 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1920 // take into account boundary condition
1921 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1923 // find the difference in z (between linear projection from layer 1
1924 // and the actual point: Dzeta= z1/r1*r2 -z2)
1925 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1926 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1928 // make "elliptical" cut in Phi and Zeta!
1929 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
1930 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
1931 if (d<1) {ret=kTRUE; break;}
1933 if(layer==0) { // try to see if it is reconstructable at the inner layer
1935 // find the difference in angles
1936 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1937 // take into account boundary condition
1938 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1940 // find the difference in z (between linear projection from layer 2
1941 // and the actual point: Dzeta= z2/r2*r1 -z1)
1942 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1943 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1945 // make "elliptical" cut in Phi and Zeta!
1946 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1947 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1948 if (d<1) {ret=kTRUE; break;};
1954 //_________________________________________________________________________
1955 void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
1957 // this method apply a rotation by 180 degree around the Z (beam) axis to all
1958 // the RecPoints in a given layer to be used to build tracklets.
1959 // **************** VERY IMPORTANT:: ***************
1960 // It must be called just after LoadClusterArrays, since afterwards the datamember
1961 // fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
1962 // instead of Cartesian
1964 if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1965 AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1967 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1968 fClustersLay1[iC1][0]*=-1;
1969 fClustersLay1[iC1][1]*=-1;
1973 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1974 fClustersLay2[iC2][0]*=-1;
1975 fClustersLay2[iC2][1]*=-1;
1980 //____________________________________________________________________________
1981 Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *esd){
1982 // This method is used to find the tracklets.
1983 // It is called from AliReconstruction
1984 // The vertex is supposed to be associated to the Tracker (i.e. to this) already
1985 // The cluster is supposed to be associated to the Tracker already
1986 // In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
1989 // apply cuts on the vertex quality
1990 const AliESDVertex *vertex = esd->GetVertex();
1991 if(vertex->GetNContributors()<fMinContVtx) return 0;
1993 AliRunLoader* runLoader = AliRunLoader::Instance();
1995 Error("Clusters2Tracks", "no run loader found");
1998 AliStack *pStack=0x0; TTree *tRefTree=0x0;
2000 runLoader->LoadKinematics("read");
2001 runLoader->LoadTrackRefs("read");
2002 pStack= runLoader->Stack();
2003 tRefTree= runLoader->TreeTR();
2005 Reconstruct(pStack,tRefTree);
2007 if (GetLightBkgStudyInParallel()) {
2008 AliStack *dummy1=0x0; TTree *dummy2=0x0;
2009 ReflectClusterAroundZAxisForLayer(1);
2010 Reconstruct(dummy1,dummy2,kTRUE);
2014 //____________________________________________________________________________
2015 Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
2017 // It is called from AliReconstruction
2023 if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
2024 if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
2025 if(GetLightBkgStudyInParallel()) {
2026 TString name="AliITSPlaneEffSPDtrackletBkg.root";
2027 TFile* pefile = TFile::Open(name, "RECREATE");
2028 rc*=fPlaneEffBkg->Write();
2033 //____________________________________________________________________
2035 AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
2037 // - gets the clusters from the cluster tree
2038 // - convert them into global coordinates
2039 // - store them in the internal arrays
2040 // - count the number of cluster-fired chips
2042 //AliDebug(1,"Loading clusters and cluster-fired chips ...");
2047 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
2048 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
2050 itsClusterBranch->SetAddress(&itsClusters);
2052 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
2053 Float_t cluGlo[3]={0.,0.,0.};
2055 // loop over the its subdetectors
2056 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
2058 if (!itsClusterTree->GetEvent(iIts))
2061 Int_t nClusters = itsClusters->GetEntriesFast();
2063 // number of clusters in each chip of the current module
2066 // loop over clusters
2067 while(nClusters--) {
2068 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
2070 layer = cluster->GetLayer();
2071 if (layer>1) continue;
2073 cluster->GetGlobalXYZ(cluGlo);
2074 Float_t x = cluGlo[0];
2075 Float_t y = cluGlo[1];
2076 Float_t z = cluGlo[2];
2079 fClustersLay1[fNClustersLay1][0] = x;
2080 fClustersLay1[fNClustersLay1][1] = y;
2081 fClustersLay1[fNClustersLay1][2] = z;
2083 for (Int_t i=0; i<3; i++)
2084 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
2087 Int_t det=cluster->GetDetectorIndex();
2088 if(det<0 || det>79) {AliError("Cluster with det. index out of boundaries"); return;}
2089 fhClustersInModuleLay1[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2093 fClustersLay2[fNClustersLay2][0] = x;
2094 fClustersLay2[fNClustersLay2][1] = y;
2095 fClustersLay2[fNClustersLay2][2] = z;
2097 for (Int_t i=0; i<3; i++)
2098 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
2101 Int_t det=cluster->GetDetectorIndex();
2102 if(det<0 || det>159) {AliError("Cluster with det. index out of boundaries"); return;}
2103 fhClustersInModuleLay2[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2107 }// end of cluster loop
2109 } // end of its "subdetector" loop
2111 itsClusters->Delete();
2115 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
2117 //_________________________________________________________________________
2119 AliITSTrackleterSPDEff::SetLightBkgStudyInParallel(Bool_t b) {
2121 // - set Bool_t fLightBackgroundStudyInParallel = b
2122 // a) if you set this kTRUE, then the estimation of the
2123 // SPD efficiency is done as usual for data, but in
2124 // parallel a light (i.e. without control histograms, etc.)
2125 // evaluation of combinatorial background is performed
2126 // with the usual ReflectClusterAroundZAxisForLayer method.
2127 // b) if you set this kFALSE, then you would not have a second
2128 // container for PlaneEfficiency statistics to be used for background
2129 // (fPlaneEffBkg=0). If you want to have a full evaluation of the
2130 // background (with all control histograms and additional data
2131 // members referring to the background) then you have to call the
2132 // method SetReflectClusterAroundZAxisForLayer(kTRUE) esplicitily
2133 fLightBkgStudyInParallel=b;
2134 if(fLightBkgStudyInParallel) {
2135 if(!fPlaneEffBkg) fPlaneEffBkg = new AliITSPlaneEffSPD();
2138 delete fPlaneEffBkg;
2142 //______________________________________________________________
2143 void AliITSTrackleterSPDEff::SetReflectClusterAroundZAxisForLayer(Int_t ilayer,Bool_t b){
2145 // method to study residual background:
2146 // Input b= KTRUE --> reflect the clusters
2147 // ilayer (either 0 or 1) --> which SPD layers should be reflected
2149 if(b) {AliInfo(Form("All clusters on layer %d will be rotated by 180 deg around z",ilayer));
2150 SetLightBkgStudyInParallel(kFALSE);}
2151 if(ilayer==0) fReflectClusterAroundZAxisForLayer0=b; // a rotation by 180degree around the Z axis
2152 else if(ilayer==1) fReflectClusterAroundZAxisForLayer1=b; // (x->-x; y->-y) to all RecPoints on a
2153 else AliInfo("Nothing done: input argument (ilayer) either 0 or 1"); // given layer is applied. In such a way