1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 //____________________________________________________________________
17 // AliITSTrackleterSPDEff - find SPD chips efficiencies by using tracklets.
19 // This class has been developed from AliITSMultReconstructor (see
20 // it for more details). It is the class for the Trackleter used to estimate
21 // SPD plane efficiency.
22 // The trackleter prediction is built using the vertex and 1 cluster.
25 // Author : Giuseppe Eugenio Bruno, based on the skeleton of Reconstruct method provided by Tiziano Virgili
26 // email: giuseppe.bruno@ba.infn.it
28 //____________________________________________________________________
34 #include <TParticle.h>
35 #include <TParticlePDG.h>
37 #include <Riostream.h>
38 #include <TClonesArray.h>
40 #include "AliTracker.h"
41 #include "AliITSTrackleterSPDEff.h"
42 #include "AliITSgeomTGeo.h"
44 #include "AliITSPlaneEffSPD.h"
46 #include "AliTrackReference.h"
47 #include "AliRunLoader.h"
48 #include "AliITSReconstructor.h"
49 #include "AliITSRecPoint.h"
50 //____________________________________________________________________
51 ClassImp(AliITSTrackleterSPDEff)
54 //____________________________________________________________________
55 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
65 fOnlyOneTrackletPerC2(0),
72 fhClustersDThetaAcc(0),
73 fhClustersDZetaAcc(0),
75 fhClustersDThetaAll(0),
76 fhClustersDZetaAll(0),
92 fOnlyOneTrackletPerC1(0),
93 fUpdateOncePerEventPlaneEff(0),
94 fChipUpdatedInEvent(0),
96 fReflectClusterAroundZAxisForLayer0(kFALSE),
97 fReflectClusterAroundZAxisForLayer1(kFALSE),
99 fUseOnlyPrimaryForPred(1),
100 fUseOnlySecondaryForPred(0),
101 fUseOnlySameParticle(0),
102 fUseOnlyDifferentParticle(0),
103 fUseOnlyStableParticle(0),
104 fPredictionPrimary(0),
105 fPredictionSecondary(0),
107 fClusterSecondary(0),
116 fhClustersDPhiInterpAcc(0),
117 fhClustersDThetaInterpAcc(0),
118 fhClustersDZetaInterpAcc(0),
119 fhClustersDPhiInterpAll(0),
120 fhClustersDThetaInterpAll(0),
121 fhClustersDZetaInterpAll(0),
122 fhDPhiVsDThetaInterpAll(0),
123 fhDPhiVsDThetaInterpAcc(0),
124 fhDPhiVsDZetaInterpAll(0),
125 fhDPhiVsDZetaInterpAcc(0),
126 fhetaClustersLay2(0),
129 // default constructor
130 // from AliITSMultReconstructor
133 SetOnlyOneTrackletPerC2();
134 fClustersLay1 = new Float_t*[300000];
135 fClustersLay2 = new Float_t*[300000];
136 fTracklets = new Float_t*[300000];
137 fAssociationFlag = new Bool_t[300000];
141 SetOnlyOneTrackletPerC1();
143 fAssociationFlag1 = new Bool_t[300000];
144 fChipPredOnLay2 = new UInt_t[300000];
145 fChipPredOnLay1 = new UInt_t[300000];
146 fChipUpdatedInEvent = new Bool_t[1200];
148 for(Int_t i=0; i<300000; i++) {
149 // from AliITSMultReconstructor
150 fClustersLay1[i] = new Float_t[6];
151 fClustersLay2[i] = new Float_t[6];
152 fTracklets[i] = new Float_t[5];
153 fAssociationFlag[i] = kFALSE;
155 fAssociationFlag1[i] = kFALSE;
157 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
159 if (GetHistOn()) BookHistos();
161 fPlaneEffSPD = new AliITSPlaneEffSPD();
163 //______________________________________________________________________
164 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
166 // from AliITSMultReconstructor
167 fClustersLay1(mr.fClustersLay1),
168 fClustersLay2(mr.fClustersLay2),
169 fTracklets(mr.fTracklets),
170 fAssociationFlag(mr.fAssociationFlag),
171 fNClustersLay1(mr.fNClustersLay1),
172 fNClustersLay2(mr.fNClustersLay2),
173 fNTracklets(mr.fNTracklets),
174 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
175 fPhiWindow(mr.fPhiWindow),
176 fZetaWindow(mr.fZetaWindow),
177 fPhiOverlapCut(mr.fPhiOverlapCut),
178 fZetaOverlapCut(mr.fZetaOverlapCut),
180 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
181 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
182 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
183 fhClustersDPhiAll(mr.fhClustersDPhiAll),
184 fhClustersDThetaAll(mr.fhClustersDThetaAll),
185 fhClustersDZetaAll(mr.fhClustersDZetaAll),
186 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
187 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
188 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
189 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
190 fhetaTracklets(mr.fhetaTracklets),
191 fhphiTracklets(mr.fhphiTracklets),
192 fhetaClustersLay1(mr.fhetaClustersLay1),
193 fhphiClustersLay1(mr.fhphiClustersLay1),
195 fAssociationFlag1(mr.fAssociationFlag1),
196 fChipPredOnLay2(mr.fChipPredOnLay2),
197 fChipPredOnLay1(mr.fChipPredOnLay1),
198 fNTracklets1(mr.fNTracklets1),
199 fPhiWindowL1(mr.fPhiWindowL1),
200 fZetaWindowL1(mr.fZetaWindowL1),
201 fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
202 fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
203 fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
204 fPlaneEffSPD(mr.fPlaneEffSPD),
205 fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
206 fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
208 fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
209 fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
210 fUseOnlySameParticle(mr.fUseOnlySameParticle),
211 fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
212 fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
213 fPredictionPrimary(mr.fPredictionPrimary),
214 fPredictionSecondary(mr.fPredictionSecondary),
215 fClusterPrimary(mr.fClusterPrimary),
216 fClusterSecondary(mr.fClusterSecondary),
217 fSuccessPP(mr.fSuccessPP),
218 fSuccessTT(mr.fSuccessTT),
219 fSuccessS(mr.fSuccessS),
220 fSuccessP(mr.fSuccessP),
221 fFailureS(mr.fFailureS),
222 fFailureP(mr.fFailureP),
224 fNonRecons(mr.fNonRecons),
225 fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
226 fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
227 fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
228 fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
229 fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
230 fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
231 fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
232 fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
233 fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
234 fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
235 fhetaClustersLay2(mr.fhetaClustersLay2),
236 fhphiClustersLay2(mr.fhphiClustersLay2)
241 //______________________________________________________________________
242 AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
243 // Assignment operator
244 this->~AliITSTrackleterSPDEff();
245 new(this) AliITSTrackleterSPDEff(mr);
248 //______________________________________________________________________
249 AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
251 // from AliITSMultReconstructor
253 for(Int_t i=0; i<300000; i++) {
254 delete [] fClustersLay1[i];
255 delete [] fClustersLay2[i];
256 delete [] fTracklets[i];
258 delete [] fClustersLay1;
259 delete [] fClustersLay2;
260 delete [] fTracklets;
261 delete [] fAssociationFlag;
266 delete [] fAssociationFlag1;
268 delete [] fChipPredOnLay2;
269 delete [] fChipPredOnLay1;
271 delete [] fChipUpdatedInEvent;
273 delete [] fPredictionPrimary;
274 delete [] fPredictionSecondary;
275 delete [] fClusterPrimary;
276 delete [] fClusterSecondary;
277 delete [] fSuccessPP;
278 delete [] fSuccessTT;
284 delete [] fNonRecons;
289 //____________________________________________________________________
291 //AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t*, AliStack *pStack, TTree *tRef) {
292 AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef) {
294 // - you have to take care of the following, before of using Reconstruct
295 // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
296 // and convert the cluster coordinates to theta, phi (seen from the
297 // interaction vertex).
298 // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
299 // - Find the extrapolation/interpolation point.
300 // - Find the chip corresponding to that
301 // - Check if there is a cluster near that point
305 // retrieve the vertex position
307 vtx[0]=(Float_t)GetX();
308 vtx[1]=(Float_t)GetY();
309 vtx[2]=(Float_t)GetZ();
310 // to study residual background (i.e. contribution from TT' to measured efficiency)
311 if(fReflectClusterAroundZAxisForLayer0) ReflectClusterAroundZAxisForLayer(0);
312 if(fReflectClusterAroundZAxisForLayer1) ReflectClusterAroundZAxisForLayer(1);
314 if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
315 if(fMC && !tRef) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
317 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
318 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
319 Int_t nfClu1=0; Int_t ntClu1=0;
320 Int_t nfClu2=0; Int_t ntClu2=0;
322 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
323 // for this new event)
324 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
326 // find the tracklets
327 AliDebug(1,"Looking for tracklets... ");
328 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
330 //###########################################################
331 // Loop on layer 1 : finding theta, phi and z
333 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
334 Float_t x = fClustersLay1[iC1][0] - vtx[0];
335 Float_t y = fClustersLay1[iC1][1] - vtx[1];
336 Float_t z = fClustersLay1[iC1][2] - vtx[2];
338 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
340 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
341 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
342 fClustersLay1[iC1][2] = z; // Store z
344 // find the Radius and the chip corresponding to the extrapolation point
346 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
348 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
349 key=999999; // also some other actions should be taken if not Found
351 nfTraPred2+=(Int_t)found; // this for debugging purpose
352 ntTraPred2++; // to check efficiency of the method FindChip
353 fChipPredOnLay2[iC1] = key;
354 fAssociationFlag1[iC1] = kFALSE;
357 Float_t eta=fClustersLay1[iC1][0];
358 eta= TMath::Tan(eta/2.);
359 eta=-TMath::Log(eta);
360 fhetaClustersLay1->Fill(eta);
361 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
364 // Loop on layer 2 : finding theta, phi and r
365 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
366 Float_t x = fClustersLay2[iC2][0] - vtx[0];
367 Float_t y = fClustersLay2[iC2][1] - vtx[1];
368 Float_t z = fClustersLay2[iC2][2] - vtx[2];
370 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
372 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
373 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
374 fClustersLay2[iC2][2] = z; // Store z
376 // find the Radius and the chip corresponding to the extrapolation point
378 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
380 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
383 nfTraPred1+=(Int_t)found; // this for debugging purpose
384 ntTraPred1++; // to check efficiency of the method FindChip
385 fChipPredOnLay1[iC2] = key;
386 fAssociationFlag[iC2] = kFALSE;
389 Float_t eta=fClustersLay2[iC2][0];
390 eta= TMath::Tan(eta/2.);
391 eta=-TMath::Log(eta);
392 fhetaClustersLay2->Fill(eta);
393 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
397 //###########################################################
399 // First part : Extrapolation to Layer 2
402 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
404 // here the control to check whether the efficiency of the chip traversed by this tracklet
405 // prediction has already been updated in this event using another tracklet prediction
406 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
408 // reset of variables for multiple candidates
409 Int_t iC2WithBestDist = 0; // reset
410 Float_t distmin = 100.; // just to put a huge number!
411 Float_t dPhimin = 0.; // Used for histograms only!
412 Float_t dThetamin = 0.; // Used for histograms only!
413 Float_t dZetamin = 0.; // Used for histograms only!
415 // in any case, if MC has been required, store statistics of primaries and secondaries
416 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
418 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
419 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
420 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
421 // do it always as a function of the chip number used to built the prediction
422 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
423 if (!found) {AliWarning(
424 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
426 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
427 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
428 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
429 { // this cluster is from a primary particle
430 fClusterPrimary[key]++;
432 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
433 } else { // this cluster is from a secondary particle
434 fClusterSecondary[key]++;
436 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
439 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
440 // (in case the prediction is reliable)
441 if( fChipPredOnLay2[iC1]<1200) {
442 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
443 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
444 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
445 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
446 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
447 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
448 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
449 else fNonRecons[fChipPredOnLay2[iC1]]++;
454 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
456 // The following excludes double associations
457 if (!fAssociationFlag[iC2]) {
459 // find the difference in angles
460 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
461 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
462 // take into account boundary condition
463 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
465 // find the difference in z (between linear projection from layer 1
466 // and the actual point: Dzeta= z1/r1*r2 -z2)
467 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
468 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
471 fhClustersDPhiAll->Fill(dPhi);
472 fhClustersDThetaAll->Fill(dTheta);
473 fhClustersDZetaAll->Fill(dZeta);
474 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
475 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
478 // make "elliptical" cut in Phi and Zeta!
479 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
480 dZeta*dZeta/fZetaWindow/fZetaWindow);
484 //look for the minimum distance: the minimum is in iC2WithBestDist
485 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
486 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
490 iC2WithBestDist = iC2;
493 } // end of loop over clusters in layer 2
495 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
498 fhClustersDPhiAcc->Fill(dPhimin);
499 fhClustersDThetaAcc->Fill(dThetamin);
500 fhClustersDZetaAcc->Fill(dZetamin);
501 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
502 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
505 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
506 // flag the association
508 // store the tracklet
510 // use the theta from the clusters in the first layer
511 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
512 // use the phi from the clusters in the first layer
513 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
514 // Store the difference between phi1 and phi2
515 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
522 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
534 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
538 fTracklets[fNTracklets][3] = -2;
542 Float_t eta=fTracklets[fNTracklets][0];
543 eta= TMath::Tan(eta/2.);
544 eta=-TMath::Log(eta);
545 fhetaTracklets->Fill(eta);
546 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
549 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
550 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
553 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
556 nfClu2+=(Int_t)found; // this for debugging purpose
557 ntClu2++; // to check efficiency of the method FindChip
558 if(key<1200) { // the Chip has been found
559 if(fMC) { // this part only for MC
560 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
561 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
562 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
565 if(primary) fSuccessPP[key]++;
567 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
568 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
571 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
573 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
574 fChipUpdatedInEvent[key]=kTRUE;
576 if(primary) fSuccessP[key]++;
577 if(secondary) fSuccessS[key]++;
581 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
582 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
584 if(primary) fSuccessP[key]++;
585 if(secondary) fSuccessS[key]++;
592 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
593 else if (fChipPredOnLay2[iC1]<1200) {
594 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
595 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
597 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
598 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
601 } // end of loop over clusters in layer 1
603 fNTracklets1=fNTracklets;
605 //###################################################################
607 // Second part : Interpolation to Layer 1
610 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
612 // here the control to check whether the efficiency of the chip traversed by this tracklet
613 // prediction has already been updated in this event using another tracklet prediction
614 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
616 // reset of variables for multiple candidates
617 Int_t iC1WithBestDist = 0; // reset
618 Float_t distmin = 100.; // just to put a huge number!
619 Float_t dPhimin = 0.; // Used for histograms only!
620 Float_t dThetamin = 0.; // Used for histograms only!
621 Float_t dZetamin = 0.; // Used for histograms only!
623 // in any case, if MC has been required, store statistics of primaries and secondaries
624 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
626 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
627 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
628 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
629 // do it always as a function of the chip number used to built the prediction
630 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
631 if (!found) {AliWarning(
632 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
634 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
635 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
636 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
637 { // this cluster is from a primary particle
638 fClusterPrimary[key]++;
640 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
641 } else { // this cluster is from a secondary particle
642 fClusterSecondary[key]++;
644 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
647 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
648 // (in case the prediction is reliable)
649 if( fChipPredOnLay1[iC2]<1200) {
650 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
651 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
652 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
653 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
654 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
655 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
656 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
657 else fNonRecons[fChipPredOnLay1[iC2]]++;
662 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
664 // The following excludes double associations
665 if (!fAssociationFlag1[iC1]) {
667 // find the difference in angles
668 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
669 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
670 // take into account boundary condition
671 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
674 // find the difference in z (between linear projection from layer 2
675 // and the actual point: Dzeta= z2/r2*r1 -z1)
676 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
677 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
681 fhClustersDPhiInterpAll->Fill(dPhi);
682 fhClustersDThetaInterpAll->Fill(dTheta);
683 fhClustersDZetaInterpAll->Fill(dZeta);
684 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
685 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
687 // make "elliptical" cut in Phi and Zeta!
688 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
689 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
693 //look for the minimum distance: the minimum is in iC1WithBestDist
694 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
695 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
699 iC1WithBestDist = iC1;
702 } // end of loop over clusters in layer 1
704 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
707 fhClustersDPhiInterpAcc->Fill(dPhimin);
708 fhClustersDThetaInterpAcc->Fill(dThetamin);
709 fhClustersDZetaInterpAcc->Fill(dZetamin);
710 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
711 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
714 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
715 // flag the association
717 // store the tracklet
719 // use the theta from the clusters in the first layer
720 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
721 // use the phi from the clusters in the first layer
722 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
723 // Store the difference between phi1 and phi2
724 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
731 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
743 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
747 fTracklets[fNTracklets][3] = -2;
750 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
751 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
754 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
757 nfClu1+=(Int_t)found; // this for debugging purpose
758 ntClu1++; // to check efficiency of the method FindChip
760 if(fMC) { // this part only for MC
761 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
762 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
763 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
764 if (label2 < 3) { // same label
766 if(primary) fSuccessPP[key]++;
768 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
769 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
772 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
774 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
775 fChipUpdatedInEvent[key]=kTRUE;
777 if(primary) fSuccessP[key]++;
778 if(secondary) fSuccessS[key]++;
781 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
782 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
784 if(primary) fSuccessP[key]++;
785 if(secondary) fSuccessS[key]++;
792 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
793 else if (fChipPredOnLay1[iC2]<1200) {
794 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
795 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
797 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
798 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
801 } // end of loop over clusters in layer 2
803 AliDebug(1,Form("%d tracklets found", fNTracklets));
804 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
805 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
806 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
807 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
809 //____________________________________________________________________
810 Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
811 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
813 // Input: a) layer number in the range [0,1]
814 // b) vtx[3]: actual vertex
815 // c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
816 // d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
818 // Output: Unique key to locate a chip
819 // return: kTRUE if succesfull
821 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
822 Double_t r=GetRLayer(layer);
823 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
825 // set phiVtx in the range [0,2pi]
826 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
828 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
829 // of the intersection of the tracklet with the pixel layer.
830 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
831 else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
832 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
833 Double_t vtxy[2]={vtx[0],vtx[1]};
834 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
835 // this method gives you two interceptions
836 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
837 // set phiAbs in the range [0,2pi]
838 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
839 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
840 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
841 // taking the closest one to phiVtx
842 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
843 } else phiAbs=phiVtx;
844 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
846 // now you need to locate the chip within the idet detector,
847 // starting from the local coordinates in such a detector
849 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
850 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
851 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
853 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
856 //______________________________________________________________________________
857 Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
859 // Return the average radius of a layer from Geometry
861 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
862 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
864 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
865 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
866 Double_t r=TMath::Sqrt(x*x + y*y);
868 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
869 r += TMath::Sqrt(x*x + y*y);
870 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
871 r += TMath::Sqrt(x*x + y*y);
872 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
873 r += TMath::Sqrt(x*x + y*y);
877 //______________________________________________________________________________
878 Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
879 Float_t &xloc, Float_t &zloc) {
880 // this method transform Global Cilindrical coordinates into local (i.e. module)
881 // cartesian coordinates
883 //Compute Cartesian Global Coordinate
884 Double_t xyzGlob[3],xyzLoc[3];
886 xyzGlob[0]=r*TMath::Cos(phi);
887 xyzGlob[1]=r*TMath::Sin(phi);
892 if(idet<0) return kFALSE;
894 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
895 Int_t lad = Int_t(idet/ndet) + 1;
896 Int_t det = idet - (lad-1)*ndet + 1;
898 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
900 xloc = (Float_t)xyzLoc[0];
901 zloc = (Float_t)xyzLoc[2];
905 //______________________________________________________________________________
906 Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
907 //--------------------------------------------------------------------
908 // This function finds the detector crossed by the track
909 // Input: layer in range [0,1]
910 // phi in ALICE absolute reference system
912 //--------------------------------------------------------------------
913 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
914 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
915 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
916 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
918 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
919 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
920 Double_t phiOffset=TMath::ATan2(y,x);
924 if (zOffset<0) // old geometry
925 dphi = -(phi-phiOffset);
927 dphi = phi-phiOffset;
929 if (dphi < 0) dphi += 2*TMath::Pi();
930 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
931 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
932 if (np>=nladders) np-=nladders;
933 if (np<0) np+=nladders;
935 Double_t dz=zOffset-z;
936 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
937 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
938 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
939 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
941 return np*ndetectors + nz;
943 //____________________________________________________________
944 Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
945 // this method find the intersection in xy between a tracklet (straight line) and
946 // a circonference (r=R), using polar coordinates.
948 Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
949 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
950 - R: radius of the circle
951 Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
953 Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
954 r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
955 In the same system, the equation of a semi-line is: phi=phiVtx;
956 Hence you get one interception only: P=(r,phiVtx)
957 Finally you want P in the ABSOLUTE ALICE reference system.
959 Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
960 Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
961 Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
962 Double_t cC=rO*rO-R*R;
963 Double_t dDelta=bB*bB-4*cC;
964 if(dDelta<0) return kFALSE;
966 r1=(-bB-TMath::Sqrt(dDelta))/2;
967 r2=(-bB+TMath::Sqrt(dDelta))/2;
968 if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
969 Double_t r=TMath::Max(r1,r2); // take the positive
970 Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
971 Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
972 pvtx[0]=r*TMath::Cos(phiVtx);
973 pvtx[1]=r*TMath::Sin(phiVtx);
974 pP[0]=vtx[0]+pvtx[0];
975 pP[1]=vtx[1]+pvtx[1];
976 phi=TMath::ATan2(pP[1],pP[0]);
979 //___________________________________________________________
980 Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
982 // simple method to reduce all angles (in rad)
986 while(angle >=2*TMath::Pi() || angle<0) {
987 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
988 if(angle < 0) angle+=2*TMath::Pi();
992 //___________________________________________________________
993 Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
995 // This method check if a particle is primary; i.e.
996 // it comes from the main vertex and it is a "stable" particle, according to
997 // AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
998 // a stable particle: it has no effect on this analysis).
999 // This method can be called only for MC events, where Kinematics is available.
1000 // if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
1001 // returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
1002 // The latter (see below) try to verify if a primary particle is also "detectable".
1004 if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
1005 if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
1006 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
1007 // return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
1008 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
1009 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
1010 TParticle* part = stack->Particle(ipart);
1011 TParticle* part0 = stack->Particle(0); // first primary
1012 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
1013 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
1014 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
1016 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
1017 TParticlePDG* pdgPart = part->GetPDG();
1018 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
1020 Double_t distx = part->Vx() - part0->Vx();
1021 Double_t disty = part->Vy() - part0->Vy();
1022 Double_t distz = part->Vz() - part0->Vz();
1023 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
1025 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
1026 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
1027 return kFALSE; }// primary if within 500 microns from true Vertex
1029 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
1032 //_____________________________________________________________________________________________
1033 Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
1035 // This private method can be applied on MC particles (if stack is available),
1036 // provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
1038 // It define "detectable" a primary particle according to the following criteria:
1040 // - if no decay products can be found in the stack (note that this does not
1041 // means it is stable, since a particle is stored in stack if it has at least 1 hit in a
1042 // sensitive detector)
1043 // - if it has at least one decay daughter produced outside or just on the outer pixel layer
1044 // - if the last decay particle is an electron (or a muon) which is not produced in-between
1045 // the two pixel layers (this is likely to be a kink).
1046 if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
1047 if(!stack) {AliError("null pointer to MC stack"); return 0;}
1048 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
1050 TParticle* part = stack->Particle(ipart);
1051 //TParticle* part0 = stack->Particle(0); // first primary
1057 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
1058 // its real fate ! But you have to take it !
1059 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
1060 Int_t lastDau = part->GetLastDaughter();
1061 nDau = lastDau - firstDau + 1;
1062 Double_t distMax=0.;
1064 for(Int_t j=firstDau; j<=lastDau; j++) {
1065 dau = stack->Particle(j);
1066 Double_t distx = dau->Vx();
1067 Double_t disty = dau->Vy();
1068 //Double_t distz = dau->Vz();
1069 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
1070 if(distR<distMax) continue; // considere only the daughter produced at largest radius
1074 dau = stack->Particle(jmax);
1075 pdgDau=dau->GetPdgCode();
1076 if (pdgDau == 11 || pdgDau == 13 ) {
1077 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
1078 else nret =0; // delta-ray emission in material (keep it)
1080 else {// not ele or muon
1081 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
1085 //_________________________________________________________________
1086 void AliITSTrackleterSPDEff::InitPredictionMC() {
1088 // this method allocate memory for the MC related informations
1089 // all the counters are set to 0
1092 if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1093 fPredictionPrimary = new Int_t[1200];
1094 fPredictionSecondary = new Int_t[1200];
1095 fClusterPrimary = new Int_t[1200];
1096 fClusterSecondary = new Int_t[1200];
1097 fSuccessPP = new Int_t[1200];
1098 fSuccessTT = new Int_t[1200];
1099 fSuccessS = new Int_t[1200];
1100 fSuccessP = new Int_t[1200];
1101 fFailureS = new Int_t[1200];
1102 fFailureP = new Int_t[1200];
1103 fRecons = new Int_t[1200];
1104 fNonRecons = new Int_t[1200];
1105 for(Int_t i=0; i<1200; i++) {
1106 fPredictionPrimary[i]=0;
1107 fPredictionSecondary[i]=0;
1108 fPredictionSecondary[i]=0;
1109 fClusterSecondary[i]=0;
1121 //_________________________________________________________________
1122 void AliITSTrackleterSPDEff::DeletePredictionMC() {
1124 // this method deallocate memory for the MC related informations
1125 // all the counters are set to 0
1128 if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1129 if(fPredictionPrimary) {
1130 delete fPredictionPrimary; fPredictionPrimary=0;
1132 if(fPredictionSecondary) {
1133 delete fPredictionSecondary; fPredictionSecondary=0;
1135 if(fClusterPrimary) {
1136 delete fClusterPrimary; fClusterPrimary=0;
1138 if(fClusterSecondary) {
1139 delete fClusterSecondary; fClusterSecondary=0;
1142 delete fSuccessPP; fSuccessPP=0;
1145 delete fSuccessTT; fSuccessTT=0;
1148 delete fSuccessS; fSuccessS=0;
1151 delete fSuccessP; fSuccessP=0;
1154 delete fFailureS; fFailureS=0;
1157 delete fFailureP; fFailureP=0;
1160 delete fRecons; fRecons=0;
1163 delete fNonRecons; fNonRecons=0;
1167 //______________________________________________________________________
1168 Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
1170 // This method return the Data menmber fPredictionPrimary [1200].
1171 // You can call it only for MC events.
1172 // fPredictionPrimary[key] contains the number of tracklet predictions on the
1173 // given chip key built using a cluster on the other layer produced (at least)
1174 // from a primary particle.
1175 // Key refers to the chip crossed by the prediction
1178 if (!fMC) {CallWarningMC(); return 0;}
1179 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1180 return fPredictionPrimary[(Int_t)key];
1182 //______________________________________________________________________
1183 Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
1185 // This method return the Data menmber fPredictionSecondary [1200].
1186 // You can call it only for MC events.
1187 // fPredictionSecondary[key] contains the number of tracklet predictions on the
1188 // given chip key built using a cluster on the other layer produced (only)
1189 // from a secondary particle
1190 // Key refers to the chip crossed by the prediction
1193 if (!fMC) {CallWarningMC(); return 0;}
1194 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1195 return fPredictionSecondary[(Int_t)key];
1197 //______________________________________________________________________
1198 Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
1200 // This method return the Data menmber fClusterPrimary [1200].
1201 // You can call it only for MC events.
1202 // fClusterPrimary[key] contains the number of tracklet predictions
1203 // built using a cluster on that layer produced (only)
1204 // from a primary particle
1205 // Key refers to the chip used to build the prediction
1208 if (!fMC) {CallWarningMC(); return 0;}
1209 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1210 return fClusterPrimary[(Int_t)key];
1212 //______________________________________________________________________
1213 Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
1215 // This method return the Data menmber fClusterSecondary [1200].
1216 // You can call it only for MC events.
1217 // fClusterSecondary[key] contains the number of tracklet predictions
1218 // built using a cluster on that layer produced (only)
1219 // from a secondary particle
1220 // Key refers to the chip used to build the prediction
1222 if (!fMC) {CallWarningMC(); return 0;}
1223 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1224 return fClusterSecondary[(Int_t)key];
1226 //______________________________________________________________________
1227 Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1229 // This method return the Data menmber fSuccessPP [1200].
1230 // You can call it only for MC events.
1231 // fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1232 // with a cluster on the other layer) built by using the same primary particle
1233 // the unique chip key refers to the chip which get updated its efficiency
1235 if (!fMC) {CallWarningMC(); return 0;}
1236 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1237 return fSuccessPP[(Int_t)key];
1239 //______________________________________________________________________
1240 Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1242 // This method return the Data menmber fSuccessTT [1200].
1243 // You can call it only for MC events.
1244 // fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1245 // with a cluster on the other layer) built by using the same particle (whatever)
1246 // the unique chip key refers to the chip which get updated its efficiency
1248 if (!fMC) {CallWarningMC(); return 0;}
1249 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1250 return fSuccessTT[(Int_t)key];
1252 //______________________________________________________________________
1253 Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1255 // This method return the Data menmber fSuccessS [1200].
1256 // You can call it only for MC events.
1257 // fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1258 // with a cluster on the other layer) built by using a secondary particle
1259 // the unique chip key refers to the chip which get updated its efficiency
1261 if (!fMC) {CallWarningMC(); return 0;}
1262 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1263 return fSuccessS[(Int_t)key];
1265 //______________________________________________________________________
1266 Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1268 // This method return the Data menmber fSuccessP [1200].
1269 // You can call it only for MC events.
1270 // fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1271 // with a cluster on the other layer) built by using a primary particle
1272 // the unique chip key refers to the chip which get updated its efficiency
1274 if (!fMC) {CallWarningMC(); return 0;}
1275 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1276 return fSuccessP[(Int_t)key];
1278 //______________________________________________________________________
1279 Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1281 // This method return the Data menmber fFailureS [1200].
1282 // You can call it only for MC events.
1283 // fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1284 // with a cluster on the other layer) built by using a secondary particle
1285 // the unique chip key refers to the chip which get updated its efficiency
1287 if (!fMC) {CallWarningMC(); return 0;}
1288 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1289 return fFailureS[(Int_t)key];
1291 //______________________________________________________________________
1292 Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1294 // This method return the Data menmber fFailureP [1200].
1295 // You can call it only for MC events.
1296 // fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1297 // with a cluster on the other layer) built by using a primary particle
1298 // the unique chip key refers to the chip which get updated its efficiency
1300 if (!fMC) {CallWarningMC(); return 0;}
1301 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1302 return fFailureP[(Int_t)key];
1304 //_____________________________________________________________________
1305 Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1307 // This method return the Data menmber fRecons [1200].
1308 // You can call it only for MC events.
1309 // fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1310 // has an hit in the detector)
1311 // the unique chip key refers to the chip where fall the prediction
1313 if (!fMC) {CallWarningMC(); return 0;}
1314 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1315 return fRecons[(Int_t)key];
1317 //_____________________________________________________________________
1318 Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1320 // This method return the Data menmber fNonRecons [1200].
1321 // You can call it only for MC events.
1322 // fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1323 // has not any hit in the detector)
1324 // the unique chip key refers to the chip where fall the prediction
1326 if (!fMC) {CallWarningMC(); return 0;}
1327 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1328 return fNonRecons[(Int_t)key];
1330 //______________________________________________________________________
1331 void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1332 // Print out some class data values in Ascii Form to output stream
1334 // ostream *os Output stream where Ascii data is to be writen
1339 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow
1340 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1341 << " " << fUpdateOncePerEventPlaneEff << " " << fReflectClusterAroundZAxisForLayer0
1342 << " " << fReflectClusterAroundZAxisForLayer1;
1344 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1345 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1346 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1347 << " " << fUseOnlyStableParticle ;
1348 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1349 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1350 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1351 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
1352 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1353 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1354 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1355 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1356 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1357 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1358 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1359 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
1362 //______________________________________________________________________
1363 void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1364 // Read in some class data values in Ascii Form to output stream
1366 // istream *is Input stream where Ascii data is to be read in from
1373 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow
1374 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1375 >> fUpdateOncePerEventPlaneEff >> fReflectClusterAroundZAxisForLayer0
1376 >> fReflectClusterAroundZAxisForLayer1;
1377 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
1379 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1381 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1382 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1383 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1384 >> fUseOnlyStableParticle;
1385 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1386 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1387 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1388 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1389 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1390 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1391 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1392 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1393 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1394 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1395 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1396 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1400 //______________________________________________________________________
1401 ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1402 // Standard output streaming function
1404 // ostream &os output steam
1405 // AliITSTrackleterSPDEff &s class to be streamed.
1409 // ostream &os The stream pointer
1414 //______________________________________________________________________
1415 istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1416 // Standard inputput streaming function
1418 // istream &is input steam
1419 // AliITSTrackleterSPDEff &s class to be streamed.
1423 // ostream &os The stream pointer
1425 //printf("prova %d \n", (Int_t)s.GetMC());
1429 //______________________________________________________________________
1430 void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
1432 // This Method write into an either asci or root file
1433 // the used cuts and the statistics of the MC related quantities
1434 // The method SetMC() has to be called before
1435 // Input TString filename: name of file for output (it deletes already existing
1440 //if(!fMC) {CallWarningMC(); return;}
1441 if (!filename.Contains(".root")) {
1442 ofstream out(filename.Data(),ios::out | ios::binary);
1448 TFile* mcfile = TFile::Open(filename, "RECREATE");
1449 TH1F* cuts = new TH1F("cuts", "list of cuts", 10, 0, 10); // TH1I containing cuts
1450 cuts->SetBinContent(1,fPhiWindowL1);
1451 cuts->SetBinContent(2,fZetaWindowL1);
1452 cuts->SetBinContent(3,fPhiWindow);
1453 cuts->SetBinContent(4,fZetaWindow);
1454 cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
1455 cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
1456 cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
1457 cuts->SetBinContent(8,fReflectClusterAroundZAxisForLayer0);
1458 cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer1);
1459 cuts->SetBinContent(10,fMC);
1462 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
1464 TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
1465 mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
1466 mc0->SetBinContent(2,fUseOnlySecondaryForPred);
1467 mc0->SetBinContent(3,fUseOnlySameParticle);
1468 mc0->SetBinContent(4,fUseOnlyDifferentParticle);
1469 mc0->SetBinContent(5,fUseOnlyStableParticle);
1473 mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
1474 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
1476 mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
1477 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
1479 mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
1480 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
1482 mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
1483 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
1485 mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
1486 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
1488 mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
1489 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
1491 mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
1492 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
1494 mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
1495 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
1497 mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
1498 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
1500 mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
1501 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
1503 mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
1504 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
1506 mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
1507 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
1515 //____________________________________________________________________
1516 void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
1518 // This Method read from an asci file (do not know why binary does not work)
1519 // the cuts to be used and the statistics of the MC related quantities
1520 // Input TString filename: name of input file for output
1521 // The method SetMC() has to be called before
1525 //if(!fMC) {CallWarningMC(); return;}
1526 if( gSystem->AccessPathName( filename.Data() ) ) {
1527 AliError( Form( "file (%s) not found", filename.Data() ) );
1531 if (!filename.Contains(".root")) {
1532 ifstream in(filename.Data(),ios::in | ios::binary);
1539 TFile *mcfile = TFile::Open(filename);
1540 TH1F *cuts = (TH1F*)mcfile->Get("cuts");
1541 fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
1542 fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
1543 fPhiWindow=(Float_t)cuts->GetBinContent(3);
1544 fZetaWindow=(Float_t)cuts->GetBinContent(4);
1545 fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
1546 fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
1547 fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
1548 fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(8);
1549 fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(9);
1550 fMC=(Bool_t)cuts->GetBinContent(10);
1551 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1552 else { // only if file with MC predictions
1553 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1554 TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
1555 fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
1556 fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
1557 fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
1558 fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
1559 fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
1561 mc1 =(TH1I*)mcfile->Get("mc1");
1562 for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1563 mc1 =(TH1I*)mcfile->Get("mc2");
1564 for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1565 mc1 =(TH1I*)mcfile->Get("mc3");
1566 for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1567 mc1 =(TH1I*)mcfile->Get("mc4");
1568 for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1569 mc1 =(TH1I*)mcfile->Get("mc5");
1570 for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1571 mc1 =(TH1I*)mcfile->Get("mc6");
1572 for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
1573 mc1 =(TH1I*)mcfile->Get("mc7");
1574 for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1575 mc1 =(TH1I*)mcfile->Get("mc8");
1576 for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1577 mc1 =(TH1I*)mcfile->Get("mc9");
1578 for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1579 mc1 =(TH1I*)mcfile->Get("mc10");
1580 for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1581 mc1 =(TH1I*)mcfile->Get("mc11");
1582 for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1583 mc1 =(TH1I*)mcfile->Get("mc12");
1584 for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1590 //____________________________________________________________________
1591 Bool_t AliITSTrackleterSPDEff::SaveHists() {
1592 // This (private) method save the histograms on the output file
1593 // (only if fHistOn is TRUE).
1594 // Also the histograms from the base class are saved through the
1595 // AliITSMultReconstructor::SaveHists() call
1597 if (!GetHistOn()) return kFALSE;
1599 // AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1600 fhClustersDPhiAll->Write();
1601 fhClustersDThetaAll->Write();
1602 fhClustersDZetaAll->Write();
1603 fhDPhiVsDThetaAll->Write();
1604 fhDPhiVsDZetaAll->Write();
1606 fhClustersDPhiAcc->Write();
1607 fhClustersDThetaAcc->Write();
1608 fhClustersDZetaAcc->Write();
1609 fhDPhiVsDThetaAcc->Write();
1610 fhDPhiVsDZetaAcc->Write();
1612 fhetaTracklets->Write();
1613 fhphiTracklets->Write();
1614 fhetaClustersLay1->Write();
1615 fhphiClustersLay1->Write();
1617 fhClustersDPhiInterpAll->Write();
1618 fhClustersDThetaInterpAll->Write();
1619 fhClustersDZetaInterpAll->Write();
1620 fhDPhiVsDThetaInterpAll->Write();
1621 fhDPhiVsDZetaInterpAll->Write();
1623 fhClustersDPhiInterpAcc->Write();
1624 fhClustersDThetaInterpAcc->Write();
1625 fhClustersDZetaInterpAcc->Write();
1626 fhDPhiVsDThetaInterpAcc->Write();
1627 fhDPhiVsDZetaInterpAcc->Write();
1629 fhetaClustersLay2->Write();
1630 fhphiClustersLay2->Write();
1633 //__________________________________________________________
1634 Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1636 // Saves the histograms into a tree and saves the trees into a file
1637 // Also the histograms from the base class are saved
1639 if (!GetHistOn()) return kFALSE;
1640 if (!strcmp(filename.Data(),"")) {
1641 AliWarning("WriteHistosToFile: null output filename!");
1644 TFile *hFile=new TFile(filename.Data(),option,
1645 "The File containing the histos for SPD efficiency studies with tracklets");
1646 if(!SaveHists()) return kFALSE;
1651 //____________________________________________________________
1652 void AliITSTrackleterSPDEff::BookHistos() {
1654 // This method books addtitional histograms
1655 // w.r.t. those of the base class.
1656 // In particular, the differences of cluster coordinate between the two SPD
1657 // layers are computed in the interpolation phase
1659 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1661 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
1662 fhClustersDPhiAcc->SetDirectory(0);
1663 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
1664 fhClustersDThetaAcc->SetDirectory(0);
1665 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
1666 fhClustersDZetaAcc->SetDirectory(0);
1668 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
1669 fhDPhiVsDZetaAcc->SetDirectory(0);
1670 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
1671 fhDPhiVsDThetaAcc->SetDirectory(0);
1673 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
1674 fhClustersDPhiAll->SetDirectory(0);
1675 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
1676 fhClustersDThetaAll->SetDirectory(0);
1677 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
1678 fhClustersDZetaAll->SetDirectory(0);
1680 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
1681 fhDPhiVsDZetaAll->SetDirectory(0);
1682 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
1683 fhDPhiVsDThetaAll->SetDirectory(0);
1685 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
1686 fhetaTracklets->SetDirectory(0);
1687 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
1688 fhphiTracklets->SetDirectory(0);
1689 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
1690 fhetaClustersLay1->SetDirectory(0);
1691 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
1692 fhphiClustersLay1->SetDirectory(0);
1694 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1695 fhClustersDPhiInterpAcc->SetDirectory(0);
1696 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1697 fhClustersDThetaInterpAcc->SetDirectory(0);
1698 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1699 fhClustersDZetaInterpAcc->SetDirectory(0);
1701 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1702 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1703 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1704 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1706 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1707 fhClustersDPhiInterpAll->SetDirectory(0);
1708 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1709 fhClustersDThetaInterpAll->SetDirectory(0);
1710 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1711 fhClustersDZetaInterpAll->SetDirectory(0);
1713 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1714 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1715 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1716 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1718 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1719 fhetaClustersLay2->SetDirectory(0);
1720 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1721 fhphiClustersLay2->SetDirectory(0);
1724 //____________________________________________________________
1725 void AliITSTrackleterSPDEff::DeleteHistos() {
1727 // Private method to delete Histograms from memory
1728 // it is called. e.g., by the destructor.
1730 // form AliITSMultReconstructor
1731 if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
1732 if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
1733 if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
1734 if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
1735 if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
1736 if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
1737 if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
1738 if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
1739 if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
1740 if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
1741 if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
1742 if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
1743 if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
1744 if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
1746 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1747 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1748 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1749 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1750 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1751 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1752 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1753 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1754 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1755 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1756 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1757 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1759 //_______________________________________________________________
1760 Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1761 Float_t* vtx, AliStack *stack, TTree *ref) {
1762 // This (private) method can be used only for MC events, where both AliStack and the TrackReference
1764 // It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1766 // - Int_t layer (either 0 or 1): layer which you want to chech if the tracklete can be
1768 // - Int_t iC : cluster index used to build the tracklet prediction
1769 // if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1770 // - Float_t* vtx: actual event vertex
1771 // - stack: pointer to Stack
1772 // - ref: pointer to TTRee of TrackReference
1773 Bool_t ret=kFALSE; // returned value
1774 Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1775 if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1776 if(!stack) {AliError("null pointer to MC stack"); return ret;}
1777 if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1778 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1779 if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1781 AliTrackReference *tref=0x0;
1782 Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1783 Int_t nentries = (Int_t)ref->GetEntries();
1784 TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1785 TBranch *br = ref->GetBranch("TrackReferences");
1786 br->SetAddress(&tcaRef);
1787 for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1788 br->GetEntry(itrack);
1789 Int_t nref=tcaRef->GetEntriesFast();
1790 if(nref>0) { //it is enough to look at the first one
1791 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1792 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1795 if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1796 br->GetEntry(imatch); // redundant, nevertheless ...
1797 Int_t nref=tcaRef->GetEntriesFast();
1798 for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1799 tref=(AliTrackReference*)tcaRef->At(iref);
1800 if(tref->R()>10) continue; // not SPD ref
1801 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1802 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1804 // compute the proper quantities for this tref, as was done for fClustersLay1/2
1805 Float_t x = tref->X() - vtx[0];
1806 Float_t y = tref->Y() - vtx[1];
1807 Float_t z = tref->Z() - vtx[2];
1809 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1811 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1812 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1813 trefLayExtr[2] = z; // Store z
1815 if(layer==1) { // try to see if it is reconstructable at the outer layer
1816 // find the difference in angles
1817 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1818 // take into account boundary condition
1819 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1821 // find the difference in z (between linear projection from layer 1
1822 // and the actual point: Dzeta= z1/r1*r2 -z2)
1823 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1824 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1826 // make "elliptical" cut in Phi and Zeta!
1827 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
1828 dZeta*dZeta/fZetaWindow/fZetaWindow);
1829 if (d<1) {ret=kTRUE; break;}
1831 if(layer==0) { // try to see if it is reconstructable at the inner layer
1833 // find the difference in angles
1834 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1835 // take into account boundary condition
1836 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1838 // find the difference in z (between linear projection from layer 2
1839 // and the actual point: Dzeta= z2/r2*r1 -z1)
1840 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1841 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1843 // make "elliptical" cut in Phi and Zeta!
1844 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1845 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1846 if (d<1) {ret=kTRUE; break;};
1852 //_________________________________________________________________________
1853 void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
1855 // this method apply a rotation by 180 degree around the Z (beam) axis to all
1856 // the RecPoints in a given layer to be used to build tracklets.
1857 // **************** VERY IMPORTANT:: ***************
1858 // It must be called just after LoadClusterArrays, since afterwards the datamember
1859 // fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
1860 // instead of Cartesian
1862 if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1863 AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1865 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1866 fClustersLay1[iC1][0]*=-1;
1867 fClustersLay1[iC1][1]*=-1;
1871 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1872 fClustersLay2[iC2][0]*=-1;
1873 fClustersLay2[iC2][1]*=-1;
1878 //____________________________________________________________________________
1879 Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *){
1880 // This method is used to find the tracklets.
1881 // It is called from AliReconstruction
1882 // The vertex is supposed to be associated to the Tracker (i.e. to this) already
1883 // The cluster is supposed to be associated to the Tracker already
1884 // In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
1887 AliRunLoader* runLoader = AliRunLoader::Instance();
1889 Error("Clusters2Tracks", "no run loader found");
1892 AliStack *pStack=0x0; TTree *tRefTree=0x0;
1894 runLoader->LoadKinematics("read");
1895 runLoader->LoadTrackRefs("read");
1896 pStack= runLoader->Stack();
1897 tRefTree= runLoader->TreeTR();
1899 Reconstruct(pStack,tRefTree);
1902 //____________________________________________________________________________
1903 Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
1905 // It is called from AliReconstruction
1911 if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
1912 if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
1915 //____________________________________________________________________
1917 AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
1919 // - gets the clusters from the cluster tree
1920 // - convert them into global coordinates
1921 // - store them in the internal arrays
1922 // - count the number of cluster-fired chips
1924 //AliDebug(1,"Loading clusters and cluster-fired chips ...");
1929 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
1930 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
1932 itsClusterBranch->SetAddress(&itsClusters);
1934 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
1935 Float_t cluGlo[3]={0.,0.,0.};
1937 // loop over the its subdetectors
1938 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
1940 if (!itsClusterTree->GetEvent(iIts))
1943 Int_t nClusters = itsClusters->GetEntriesFast();
1945 // number of clusters in each chip of the current module
1948 // loop over clusters
1949 while(nClusters--) {
1950 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
1952 layer = cluster->GetLayer();
1953 if (layer>1) continue;
1955 cluster->GetGlobalXYZ(cluGlo);
1956 Float_t x = cluGlo[0];
1957 Float_t y = cluGlo[1];
1958 Float_t z = cluGlo[2];
1961 fClustersLay1[fNClustersLay1][0] = x;
1962 fClustersLay1[fNClustersLay1][1] = y;
1963 fClustersLay1[fNClustersLay1][2] = z;
1965 for (Int_t i=0; i<3; i++)
1966 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
1970 fClustersLay2[fNClustersLay2][0] = x;
1971 fClustersLay2[fNClustersLay2][1] = y;
1972 fClustersLay2[fNClustersLay2][2] = z;
1974 for (Int_t i=0; i<3; i++)
1975 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
1979 }// end of cluster loop
1981 } // end of its "subdetector" loop
1983 itsClusters->Delete();
1987 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));