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>
36 #include <Riostream.h>
37 #include <TClonesArray.h>
39 #include "AliTracker.h"
40 #include "AliITSTrackleterSPDEff.h"
41 #include "AliITSgeomTGeo.h"
43 #include "AliITSPlaneEffSPD.h"
45 #include "AliTrackReference.h"
46 #include "AliRunLoader.h"
47 #include "AliITSReconstructor.h"
48 #include "AliITSRecPoint.h"
49 //____________________________________________________________________
50 ClassImp(AliITSTrackleterSPDEff)
53 //____________________________________________________________________
54 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
64 fOnlyOneTrackletPerC2(0),
71 fhClustersDThetaAcc(0),
72 fhClustersDZetaAcc(0),
74 fhClustersDThetaAll(0),
75 fhClustersDZetaAll(0),
91 fOnlyOneTrackletPerC1(0),
92 fUpdateOncePerEventPlaneEff(0),
93 fChipUpdatedInEvent(0),
95 fReflectClusterAroundZAxisForLayer0(kFALSE),
96 fReflectClusterAroundZAxisForLayer1(kFALSE),
98 fUseOnlyPrimaryForPred(1),
99 fUseOnlySecondaryForPred(0),
100 fUseOnlySameParticle(0),
101 fUseOnlyDifferentParticle(0),
102 fUseOnlyStableParticle(0),
103 fPredictionPrimary(0),
104 fPredictionSecondary(0),
106 fClusterSecondary(0),
115 fhClustersDPhiInterpAcc(0),
116 fhClustersDThetaInterpAcc(0),
117 fhClustersDZetaInterpAcc(0),
118 fhClustersDPhiInterpAll(0),
119 fhClustersDThetaInterpAll(0),
120 fhClustersDZetaInterpAll(0),
121 fhDPhiVsDThetaInterpAll(0),
122 fhDPhiVsDThetaInterpAcc(0),
123 fhDPhiVsDZetaInterpAll(0),
124 fhDPhiVsDZetaInterpAcc(0),
125 fhetaClustersLay2(0),
128 // default constructor
129 // from AliITSMultReconstructor
132 SetOnlyOneTrackletPerC2();
133 fClustersLay1 = new Float_t*[300000];
134 fClustersLay2 = new Float_t*[300000];
135 fTracklets = new Float_t*[300000];
136 fAssociationFlag = new Bool_t[300000];
140 SetOnlyOneTrackletPerC1();
142 fAssociationFlag1 = new Bool_t[300000];
143 fChipPredOnLay2 = new UInt_t[300000];
144 fChipPredOnLay1 = new UInt_t[300000];
145 fChipUpdatedInEvent = new Bool_t[1200];
147 for(Int_t i=0; i<300000; i++) {
148 // from AliITSMultReconstructor
149 fClustersLay1[i] = new Float_t[6];
150 fClustersLay2[i] = new Float_t[6];
151 fTracklets[i] = new Float_t[5];
152 fAssociationFlag[i] = kFALSE;
154 fAssociationFlag1[i] = kFALSE;
156 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
158 if (GetHistOn()) BookHistos();
160 fPlaneEffSPD = new AliITSPlaneEffSPD();
162 //______________________________________________________________________
163 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
165 // from AliITSMultReconstructor
166 fClustersLay1(mr.fClustersLay1),
167 fClustersLay2(mr.fClustersLay2),
168 fTracklets(mr.fTracklets),
169 fAssociationFlag(mr.fAssociationFlag),
170 fNClustersLay1(mr.fNClustersLay1),
171 fNClustersLay2(mr.fNClustersLay2),
172 fNTracklets(mr.fNTracklets),
173 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
174 fPhiWindow(mr.fPhiWindow),
175 fZetaWindow(mr.fZetaWindow),
176 fPhiOverlapCut(mr.fPhiOverlapCut),
177 fZetaOverlapCut(mr.fZetaOverlapCut),
179 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
180 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
181 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
182 fhClustersDPhiAll(mr.fhClustersDPhiAll),
183 fhClustersDThetaAll(mr.fhClustersDThetaAll),
184 fhClustersDZetaAll(mr.fhClustersDZetaAll),
185 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
186 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
187 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
188 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
189 fhetaTracklets(mr.fhetaTracklets),
190 fhphiTracklets(mr.fhphiTracklets),
191 fhetaClustersLay1(mr.fhetaClustersLay1),
192 fhphiClustersLay1(mr.fhphiClustersLay1),
194 fAssociationFlag1(mr.fAssociationFlag1),
195 fChipPredOnLay2(mr.fChipPredOnLay2),
196 fChipPredOnLay1(mr.fChipPredOnLay1),
197 fNTracklets1(mr.fNTracklets1),
198 fPhiWindowL1(mr.fPhiWindowL1),
199 fZetaWindowL1(mr.fZetaWindowL1),
200 fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
201 fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
202 fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
203 fPlaneEffSPD(mr.fPlaneEffSPD),
204 fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
205 fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
207 fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
208 fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
209 fUseOnlySameParticle(mr.fUseOnlySameParticle),
210 fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
211 fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
212 fPredictionPrimary(mr.fPredictionPrimary),
213 fPredictionSecondary(mr.fPredictionSecondary),
214 fClusterPrimary(mr.fClusterPrimary),
215 fClusterSecondary(mr.fClusterSecondary),
216 fSuccessPP(mr.fSuccessPP),
217 fSuccessTT(mr.fSuccessTT),
218 fSuccessS(mr.fSuccessS),
219 fSuccessP(mr.fSuccessP),
220 fFailureS(mr.fFailureS),
221 fFailureP(mr.fFailureP),
223 fNonRecons(mr.fNonRecons),
224 fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
225 fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
226 fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
227 fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
228 fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
229 fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
230 fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
231 fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
232 fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
233 fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
234 fhetaClustersLay2(mr.fhetaClustersLay2),
235 fhphiClustersLay2(mr.fhphiClustersLay2)
240 //______________________________________________________________________
241 AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
242 // Assignment operator
243 this->~AliITSTrackleterSPDEff();
244 new(this) AliITSTrackleterSPDEff(mr);
247 //______________________________________________________________________
248 AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
250 // from AliITSMultReconstructor
252 for(Int_t i=0; i<300000; i++) {
253 delete [] fClustersLay1[i];
254 delete [] fClustersLay2[i];
255 delete [] fTracklets[i];
257 delete [] fClustersLay1;
258 delete [] fClustersLay2;
259 delete [] fTracklets;
260 delete [] fAssociationFlag;
265 delete [] fAssociationFlag1;
267 delete [] fChipPredOnLay2;
268 delete [] fChipPredOnLay1;
270 delete [] fChipUpdatedInEvent;
272 delete [] fPredictionPrimary;
273 delete [] fPredictionSecondary;
274 delete [] fClusterPrimary;
275 delete [] fClusterSecondary;
276 delete [] fSuccessPP;
277 delete [] fSuccessTT;
283 delete [] fNonRecons;
288 //____________________________________________________________________
290 //AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t*, AliStack *pStack, TTree *tRef) {
291 AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef) {
293 // - you have to take care of the following, before of using Reconstruct
294 // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
295 // and convert the cluster coordinates to theta, phi (seen from the
296 // interaction vertex).
297 // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
298 // - Find the extrapolation/interpolation point.
299 // - Find the chip corresponding to that
300 // - Check if there is a cluster near that point
304 // retrieve the vertex position
306 vtx[0]=(Float_t)GetX();
307 vtx[1]=(Float_t)GetY();
308 vtx[2]=(Float_t)GetZ();
309 // to study residual background (i.e. contribution from TT' to measured efficiency)
310 if(fReflectClusterAroundZAxisForLayer0) ReflectClusterAroundZAxisForLayer(0);
311 if(fReflectClusterAroundZAxisForLayer1) ReflectClusterAroundZAxisForLayer(1);
313 if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
314 if(fMC && !tRef) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
316 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
317 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
318 Int_t nfClu1=0; Int_t ntClu1=0;
319 Int_t nfClu2=0; Int_t ntClu2=0;
321 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
322 // for this new event)
323 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
325 // find the tracklets
326 AliDebug(1,"Looking for tracklets... ");
327 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
329 //###########################################################
330 // Loop on layer 1 : finding theta, phi and z
332 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
333 Float_t x = fClustersLay1[iC1][0] - vtx[0];
334 Float_t y = fClustersLay1[iC1][1] - vtx[1];
335 Float_t z = fClustersLay1[iC1][2] - vtx[2];
337 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
339 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
340 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
341 fClustersLay1[iC1][2] = z; // Store z
343 // find the Radius and the chip corresponding to the extrapolation point
345 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
347 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
348 key=999999; // also some other actions should be taken if not Found
350 nfTraPred2+=(Int_t)found; // this for debugging purpose
351 ntTraPred2++; // to check efficiency of the method FindChip
352 fChipPredOnLay2[iC1] = key;
353 fAssociationFlag1[iC1] = kFALSE;
356 Float_t eta=fClustersLay1[iC1][0];
357 eta= TMath::Tan(eta/2.);
358 eta=-TMath::Log(eta);
359 fhetaClustersLay1->Fill(eta);
360 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
363 // Loop on layer 2 : finding theta, phi and r
364 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
365 Float_t x = fClustersLay2[iC2][0] - vtx[0];
366 Float_t y = fClustersLay2[iC2][1] - vtx[1];
367 Float_t z = fClustersLay2[iC2][2] - vtx[2];
369 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
371 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
372 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
373 fClustersLay2[iC2][2] = z; // Store z
375 // find the Radius and the chip corresponding to the extrapolation point
377 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
379 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
382 nfTraPred1+=(Int_t)found; // this for debugging purpose
383 ntTraPred1++; // to check efficiency of the method FindChip
384 fChipPredOnLay1[iC2] = key;
385 fAssociationFlag[iC2] = kFALSE;
388 Float_t eta=fClustersLay2[iC2][0];
389 eta= TMath::Tan(eta/2.);
390 eta=-TMath::Log(eta);
391 fhetaClustersLay2->Fill(eta);
392 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
396 //###########################################################
398 // First part : Extrapolation to Layer 2
401 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
403 // here the control to check whether the efficiency of the chip traversed by this tracklet
404 // prediction has already been updated in this event using another tracklet prediction
405 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
407 // reset of variables for multiple candidates
408 Int_t iC2WithBestDist = 0; // reset
409 Float_t distmin = 100.; // just to put a huge number!
410 Float_t dPhimin = 0.; // Used for histograms only!
411 Float_t dThetamin = 0.; // Used for histograms only!
412 Float_t dZetamin = 0.; // Used for histograms only!
414 // in any case, if MC has been required, store statistics of primaries and secondaries
415 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
417 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
418 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
419 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
420 // do it always as a function of the chip number used to built the prediction
421 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
422 if (!found) {AliWarning(
423 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
425 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
426 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
427 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
428 { // this cluster is from a primary particle
429 fClusterPrimary[key]++;
431 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
432 } else { // this cluster is from a secondary particle
433 fClusterSecondary[key]++;
435 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
438 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
439 // (in case the prediction is reliable)
440 if( fChipPredOnLay2[iC1]<1200) {
441 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
442 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
443 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
444 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
445 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
446 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
447 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
448 else fNonRecons[fChipPredOnLay2[iC1]]++;
453 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
455 // The following excludes double associations
456 if (!fAssociationFlag[iC2]) {
458 // find the difference in angles
459 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
460 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
461 // take into account boundary condition
462 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
464 // find the difference in z (between linear projection from layer 1
465 // and the actual point: Dzeta= z1/r1*r2 -z2)
466 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
467 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
470 fhClustersDPhiAll->Fill(dPhi);
471 fhClustersDThetaAll->Fill(dTheta);
472 fhClustersDZetaAll->Fill(dZeta);
473 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
474 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
477 // make "elliptical" cut in Phi and Zeta!
478 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
479 dZeta*dZeta/fZetaWindow/fZetaWindow);
483 //look for the minimum distance: the minimum is in iC2WithBestDist
484 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
485 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
489 iC2WithBestDist = iC2;
492 } // end of loop over clusters in layer 2
494 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
497 fhClustersDPhiAcc->Fill(dPhimin);
498 fhClustersDThetaAcc->Fill(dThetamin);
499 fhClustersDZetaAcc->Fill(dZetamin);
500 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
501 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
504 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
505 // flag the association
507 // store the tracklet
509 // use the theta from the clusters in the first layer
510 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
511 // use the phi from the clusters in the first layer
512 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
513 // Store the difference between phi1 and phi2
514 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
521 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
533 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
537 fTracklets[fNTracklets][3] = -2;
541 Float_t eta=fTracklets[fNTracklets][0];
542 eta= TMath::Tan(eta/2.);
543 eta=-TMath::Log(eta);
544 fhetaTracklets->Fill(eta);
545 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
548 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
549 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
552 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
555 nfClu2+=(Int_t)found; // this for debugging purpose
556 ntClu2++; // to check efficiency of the method FindChip
557 if(key<1200) { // the Chip has been found
558 if(fMC) { // this part only for MC
559 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
560 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
561 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
564 if(primary) fSuccessPP[key]++;
566 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
567 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
570 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
572 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
573 fChipUpdatedInEvent[key]=kTRUE;
575 if(primary) fSuccessP[key]++;
576 if(secondary) fSuccessS[key]++;
580 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
581 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
583 if(primary) fSuccessP[key]++;
584 if(secondary) fSuccessS[key]++;
591 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
592 else if (fChipPredOnLay2[iC1]<1200) {
593 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
594 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
596 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
597 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
600 } // end of loop over clusters in layer 1
602 fNTracklets1=fNTracklets;
604 //###################################################################
606 // Second part : Interpolation to Layer 1
609 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
611 // here the control to check whether the efficiency of the chip traversed by this tracklet
612 // prediction has already been updated in this event using another tracklet prediction
613 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
615 // reset of variables for multiple candidates
616 Int_t iC1WithBestDist = 0; // reset
617 Float_t distmin = 100.; // just to put a huge number!
618 Float_t dPhimin = 0.; // Used for histograms only!
619 Float_t dThetamin = 0.; // Used for histograms only!
620 Float_t dZetamin = 0.; // Used for histograms only!
622 // in any case, if MC has been required, store statistics of primaries and secondaries
623 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
625 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
626 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
627 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
628 // do it always as a function of the chip number used to built the prediction
629 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
630 if (!found) {AliWarning(
631 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
633 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
634 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
635 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
636 { // this cluster is from a primary particle
637 fClusterPrimary[key]++;
639 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
640 } else { // this cluster is from a secondary particle
641 fClusterSecondary[key]++;
643 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
646 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
647 // (in case the prediction is reliable)
648 if( fChipPredOnLay1[iC2]<1200) {
649 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
650 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
651 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
652 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
653 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
654 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
655 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
656 else fNonRecons[fChipPredOnLay1[iC2]]++;
661 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
663 // The following excludes double associations
664 if (!fAssociationFlag1[iC1]) {
666 // find the difference in angles
667 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
668 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
669 // take into account boundary condition
670 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
673 // find the difference in z (between linear projection from layer 2
674 // and the actual point: Dzeta= z2/r2*r1 -z1)
675 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
676 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
680 fhClustersDPhiInterpAll->Fill(dPhi);
681 fhClustersDThetaInterpAll->Fill(dTheta);
682 fhClustersDZetaInterpAll->Fill(dZeta);
683 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
684 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
686 // make "elliptical" cut in Phi and Zeta!
687 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
688 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
692 //look for the minimum distance: the minimum is in iC1WithBestDist
693 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
694 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
698 iC1WithBestDist = iC1;
701 } // end of loop over clusters in layer 1
703 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
706 fhClustersDPhiInterpAcc->Fill(dPhimin);
707 fhClustersDThetaInterpAcc->Fill(dThetamin);
708 fhClustersDZetaInterpAcc->Fill(dZetamin);
709 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
710 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
713 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
714 // flag the association
716 // store the tracklet
718 // use the theta from the clusters in the first layer
719 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
720 // use the phi from the clusters in the first layer
721 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
722 // Store the difference between phi1 and phi2
723 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
730 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
742 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
746 fTracklets[fNTracklets][3] = -2;
749 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
750 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
753 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
756 nfClu1+=(Int_t)found; // this for debugging purpose
757 ntClu1++; // to check efficiency of the method FindChip
759 if(fMC) { // this part only for MC
760 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
761 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
762 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
763 if (label2 < 3) { // same label
765 if(primary) fSuccessPP[key]++;
767 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
768 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
771 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
773 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
774 fChipUpdatedInEvent[key]=kTRUE;
776 if(primary) fSuccessP[key]++;
777 if(secondary) fSuccessS[key]++;
780 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
781 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
783 if(primary) fSuccessP[key]++;
784 if(secondary) fSuccessS[key]++;
791 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
792 else if (fChipPredOnLay1[iC2]<1200) {
793 fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
794 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
796 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
797 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
800 } // end of loop over clusters in layer 2
802 AliDebug(1,Form("%d tracklets found", fNTracklets));
803 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
804 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
805 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
806 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
808 //____________________________________________________________________
809 Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
810 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
812 // Input: a) layer number in the range [0,1]
813 // b) vtx[3]: actual vertex
814 // c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
815 // d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
817 // Output: Unique key to locate a chip
818 // return: kTRUE if succesfull
820 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
821 Double_t r=GetRLayer(layer);
822 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
824 // set phiVtx in the range [0,2pi]
825 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
827 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
828 // of the intersection of the tracklet with the pixel layer.
829 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
830 else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
831 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
832 Double_t vtxy[2]={vtx[0],vtx[1]};
833 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
834 // this method gives you two interceptions
835 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
836 // set phiAbs in the range [0,2pi]
837 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
838 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
839 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
840 // taking the closest one to phiVtx
841 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
842 } else phiAbs=phiVtx;
843 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
845 // now you need to locate the chip within the idet detector,
846 // starting from the local coordinates in such a detector
848 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
849 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
850 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
852 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
855 //______________________________________________________________________________
856 Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
858 // Return the average radius of a layer from Geometry
860 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
861 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
863 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
864 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
865 Double_t r=TMath::Sqrt(x*x + y*y);
867 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
868 r += TMath::Sqrt(x*x + y*y);
869 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
870 r += TMath::Sqrt(x*x + y*y);
871 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
872 r += TMath::Sqrt(x*x + y*y);
876 //______________________________________________________________________________
877 Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
878 Float_t &xloc, Float_t &zloc) {
879 // this method transform Global Cilindrical coordinates into local (i.e. module)
880 // cartesian coordinates
882 //Compute Cartesian Global Coordinate
883 Double_t xyzGlob[3],xyzLoc[3];
885 xyzGlob[0]=r*TMath::Cos(phi);
886 xyzGlob[1]=r*TMath::Sin(phi);
891 if(idet<0) return kFALSE;
893 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
894 Int_t lad = Int_t(idet/ndet) + 1;
895 Int_t det = idet - (lad-1)*ndet + 1;
897 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
899 xloc = (Float_t)xyzLoc[0];
900 zloc = (Float_t)xyzLoc[2];
904 //______________________________________________________________________________
905 Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
906 //--------------------------------------------------------------------
907 // This function finds the detector crossed by the track
908 // Input: layer in range [0,1]
909 // phi in ALICE absolute reference system
911 //--------------------------------------------------------------------
912 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
913 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
914 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
915 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
917 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
918 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
919 Double_t phiOffset=TMath::ATan2(y,x);
923 if (zOffset<0) // old geometry
924 dphi = -(phi-phiOffset);
926 dphi = phi-phiOffset;
928 if (dphi < 0) dphi += 2*TMath::Pi();
929 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
930 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
931 if (np>=nladders) np-=nladders;
932 if (np<0) np+=nladders;
934 Double_t dz=zOffset-z;
935 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
936 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
937 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
938 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
940 return np*ndetectors + nz;
942 //____________________________________________________________
943 Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
944 // this method find the intersection in xy between a tracklet (straight line) and
945 // a circonference (r=R), using polar coordinates.
947 Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
948 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
949 - R: radius of the circle
950 Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
952 Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
953 r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
954 In the same system, the equation of a semi-line is: phi=phiVtx;
955 Hence you get one interception only: P=(r,phiVtx)
956 Finally you want P in the ABSOLUTE ALICE reference system.
958 Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
959 Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
960 Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
961 Double_t cC=rO*rO-R*R;
962 Double_t dDelta=bB*bB-4*cC;
963 if(dDelta<0) return kFALSE;
965 r1=(-bB-TMath::Sqrt(dDelta))/2;
966 r2=(-bB+TMath::Sqrt(dDelta))/2;
967 if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
968 Double_t r=TMath::Max(r1,r2); // take the positive
969 Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
970 Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
971 pvtx[0]=r*TMath::Cos(phiVtx);
972 pvtx[1]=r*TMath::Sin(phiVtx);
973 pP[0]=vtx[0]+pvtx[0];
974 pP[1]=vtx[1]+pvtx[1];
975 phi=TMath::ATan2(pP[1],pP[0]);
978 //___________________________________________________________
979 Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
981 // simple method to reduce all angles (in rad)
985 while(angle >=2*TMath::Pi() || angle<0) {
986 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
987 if(angle < 0) angle+=2*TMath::Pi();
991 //___________________________________________________________
992 Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
994 // This method check if a particle is primary; i.e.
995 // it comes from the main vertex and it is a "stable" particle, according to
996 // AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
997 // a stable particle: it has no effect on this analysis).
998 // This method can be called only for MC events, where Kinematics is available.
999 // if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
1000 // returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
1001 // The latter (see below) try to verify if a primary particle is also "detectable".
1003 if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
1004 if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
1005 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
1006 // return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
1007 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
1008 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
1009 TParticle* part = stack->Particle(ipart);
1010 TParticle* part0 = stack->Particle(0); // first primary
1011 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
1012 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
1013 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
1015 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
1016 TParticlePDG* pdgPart = part->GetPDG();
1017 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
1019 Double_t distx = part->Vx() - part0->Vx();
1020 Double_t disty = part->Vy() - part0->Vy();
1021 Double_t distz = part->Vz() - part0->Vz();
1022 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
1024 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
1025 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
1026 return kFALSE; }// primary if within 500 microns from true Vertex
1028 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
1031 //_____________________________________________________________________________________________
1032 Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
1034 // This private method can be applied on MC particles (if stack is available),
1035 // provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
1037 // It define "detectable" a primary particle according to the following criteria:
1039 // - if no decay products can be found in the stack (note that this does not
1040 // means it is stable, since a particle is stored in stack if it has at least 1 hit in a
1041 // sensitive detector)
1042 // - if it has at least one decay daughter produced outside or just on the outer pixel layer
1043 // - if the last decay particle is an electron (or a muon) which is not produced in-between
1044 // the two pixel layers (this is likely to be a kink).
1045 if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
1046 if(!stack) {AliError("null pointer to MC stack"); return 0;}
1047 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
1049 TParticle* part = stack->Particle(ipart);
1050 //TParticle* part0 = stack->Particle(0); // first primary
1056 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
1057 // its real fate ! But you have to take it !
1058 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
1059 Int_t lastDau = part->GetLastDaughter();
1060 nDau = lastDau - firstDau + 1;
1061 Double_t distMax=0.;
1063 for(Int_t j=firstDau; j<=lastDau; j++) {
1064 dau = stack->Particle(j);
1065 Double_t distx = dau->Vx();
1066 Double_t disty = dau->Vy();
1067 //Double_t distz = dau->Vz();
1068 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
1069 if(distR<distMax) continue; // considere only the daughter produced at largest radius
1073 dau = stack->Particle(jmax);
1074 pdgDau=dau->GetPdgCode();
1075 if (pdgDau == 11 || pdgDau == 13 ) {
1076 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
1077 else nret =0; // delta-ray emission in material (keep it)
1079 else {// not ele or muon
1080 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
1084 //_________________________________________________________________
1085 void AliITSTrackleterSPDEff::InitPredictionMC() {
1087 // this method allocate memory for the MC related informations
1088 // all the counters are set to 0
1091 if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1092 fPredictionPrimary = new Int_t[1200];
1093 fPredictionSecondary = new Int_t[1200];
1094 fClusterPrimary = new Int_t[1200];
1095 fClusterSecondary = new Int_t[1200];
1096 fSuccessPP = new Int_t[1200];
1097 fSuccessTT = new Int_t[1200];
1098 fSuccessS = new Int_t[1200];
1099 fSuccessP = new Int_t[1200];
1100 fFailureS = new Int_t[1200];
1101 fFailureP = new Int_t[1200];
1102 fRecons = new Int_t[1200];
1103 fNonRecons = new Int_t[1200];
1104 for(Int_t i=0; i<1200; i++) {
1105 fPredictionPrimary[i]=0;
1106 fPredictionSecondary[i]=0;
1107 fPredictionSecondary[i]=0;
1108 fClusterSecondary[i]=0;
1120 //_________________________________________________________________
1121 void AliITSTrackleterSPDEff::DeletePredictionMC() {
1123 // this method deallocate memory for the MC related informations
1124 // all the counters are set to 0
1127 if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1128 if(fPredictionPrimary) {
1129 delete fPredictionPrimary; fPredictionPrimary=0;
1131 if(fPredictionSecondary) {
1132 delete fPredictionSecondary; fPredictionSecondary=0;
1134 if(fClusterPrimary) {
1135 delete fClusterPrimary; fClusterPrimary=0;
1137 if(fClusterSecondary) {
1138 delete fClusterSecondary; fClusterSecondary=0;
1141 delete fSuccessPP; fSuccessPP=0;
1144 delete fSuccessTT; fSuccessTT=0;
1147 delete fSuccessS; fSuccessS=0;
1150 delete fSuccessP; fSuccessP=0;
1153 delete fFailureS; fFailureS=0;
1156 delete fFailureP; fFailureP=0;
1159 delete fRecons; fRecons=0;
1162 delete fNonRecons; fNonRecons=0;
1166 //______________________________________________________________________
1167 Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
1169 // This method return the Data menmber fPredictionPrimary [1200].
1170 // You can call it only for MC events.
1171 // fPredictionPrimary[key] contains the number of tracklet predictions on the
1172 // given chip key built using a cluster on the other layer produced (at least)
1173 // from a primary particle.
1174 // Key refers to the chip crossed by the prediction
1177 if (!fMC) {CallWarningMC(); return 0;}
1178 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1179 return fPredictionPrimary[(Int_t)key];
1181 //______________________________________________________________________
1182 Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
1184 // This method return the Data menmber fPredictionSecondary [1200].
1185 // You can call it only for MC events.
1186 // fPredictionSecondary[key] contains the number of tracklet predictions on the
1187 // given chip key built using a cluster on the other layer produced (only)
1188 // from a secondary particle
1189 // Key refers to the chip crossed by the prediction
1192 if (!fMC) {CallWarningMC(); return 0;}
1193 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1194 return fPredictionSecondary[(Int_t)key];
1196 //______________________________________________________________________
1197 Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
1199 // This method return the Data menmber fClusterPrimary [1200].
1200 // You can call it only for MC events.
1201 // fClusterPrimary[key] contains the number of tracklet predictions
1202 // built using a cluster on that layer produced (only)
1203 // from a primary particle
1204 // Key refers to the chip used to build the prediction
1207 if (!fMC) {CallWarningMC(); return 0;}
1208 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1209 return fClusterPrimary[(Int_t)key];
1211 //______________________________________________________________________
1212 Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
1214 // This method return the Data menmber fClusterSecondary [1200].
1215 // You can call it only for MC events.
1216 // fClusterSecondary[key] contains the number of tracklet predictions
1217 // built using a cluster on that layer produced (only)
1218 // from a secondary particle
1219 // Key refers to the chip used to build the prediction
1221 if (!fMC) {CallWarningMC(); return 0;}
1222 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1223 return fClusterSecondary[(Int_t)key];
1225 //______________________________________________________________________
1226 Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1228 // This method return the Data menmber fSuccessPP [1200].
1229 // You can call it only for MC events.
1230 // fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1231 // with a cluster on the other layer) built by using the same primary particle
1232 // the unique chip key refers to the chip which get updated its efficiency
1234 if (!fMC) {CallWarningMC(); return 0;}
1235 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1236 return fSuccessPP[(Int_t)key];
1238 //______________________________________________________________________
1239 Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1241 // This method return the Data menmber fSuccessTT [1200].
1242 // You can call it only for MC events.
1243 // fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1244 // with a cluster on the other layer) built by using the same particle (whatever)
1245 // the unique chip key refers to the chip which get updated its efficiency
1247 if (!fMC) {CallWarningMC(); return 0;}
1248 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1249 return fSuccessTT[(Int_t)key];
1251 //______________________________________________________________________
1252 Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1254 // This method return the Data menmber fSuccessS [1200].
1255 // You can call it only for MC events.
1256 // fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1257 // with a cluster on the other layer) built by using a secondary particle
1258 // the unique chip key refers to the chip which get updated its efficiency
1260 if (!fMC) {CallWarningMC(); return 0;}
1261 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1262 return fSuccessS[(Int_t)key];
1264 //______________________________________________________________________
1265 Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1267 // This method return the Data menmber fSuccessP [1200].
1268 // You can call it only for MC events.
1269 // fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1270 // with a cluster on the other layer) built by using a primary particle
1271 // the unique chip key refers to the chip which get updated its efficiency
1273 if (!fMC) {CallWarningMC(); return 0;}
1274 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1275 return fSuccessP[(Int_t)key];
1277 //______________________________________________________________________
1278 Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1280 // This method return the Data menmber fFailureS [1200].
1281 // You can call it only for MC events.
1282 // fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1283 // with a cluster on the other layer) built by using a secondary particle
1284 // the unique chip key refers to the chip which get updated its efficiency
1286 if (!fMC) {CallWarningMC(); return 0;}
1287 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1288 return fFailureS[(Int_t)key];
1290 //______________________________________________________________________
1291 Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1293 // This method return the Data menmber fFailureP [1200].
1294 // You can call it only for MC events.
1295 // fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1296 // with a cluster on the other layer) built by using a primary particle
1297 // the unique chip key refers to the chip which get updated its efficiency
1299 if (!fMC) {CallWarningMC(); return 0;}
1300 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1301 return fFailureP[(Int_t)key];
1303 //_____________________________________________________________________
1304 Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1306 // This method return the Data menmber fRecons [1200].
1307 // You can call it only for MC events.
1308 // fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1309 // has an hit in the detector)
1310 // the unique chip key refers to the chip where fall the prediction
1312 if (!fMC) {CallWarningMC(); return 0;}
1313 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1314 return fRecons[(Int_t)key];
1316 //_____________________________________________________________________
1317 Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1319 // This method return the Data menmber fNonRecons [1200].
1320 // You can call it only for MC events.
1321 // fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1322 // has not any hit in the detector)
1323 // the unique chip key refers to the chip where fall the prediction
1325 if (!fMC) {CallWarningMC(); return 0;}
1326 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1327 return fNonRecons[(Int_t)key];
1329 //______________________________________________________________________
1330 void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1331 // Print out some class data values in Ascii Form to output stream
1333 // ostream *os Output stream where Ascii data is to be writen
1338 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow
1339 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
1340 << " " << fUpdateOncePerEventPlaneEff << " " << fReflectClusterAroundZAxisForLayer0
1341 << " " << fReflectClusterAroundZAxisForLayer1;
1343 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1344 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1345 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1346 << " " << fUseOnlyStableParticle ;
1347 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1348 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1349 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1350 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
1351 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1352 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1353 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1354 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1355 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1356 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1357 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1358 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
1361 //______________________________________________________________________
1362 void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1363 // Read in some class data values in Ascii Form to output stream
1365 // istream *is Input stream where Ascii data is to be read in from
1372 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow
1373 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
1374 >> fUpdateOncePerEventPlaneEff >> fReflectClusterAroundZAxisForLayer0
1375 >> fReflectClusterAroundZAxisForLayer1;
1376 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
1378 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1380 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1381 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1382 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1383 >> fUseOnlyStableParticle;
1384 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1385 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1386 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1387 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1388 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1389 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1390 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1391 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1392 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1393 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1394 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1395 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1399 //______________________________________________________________________
1400 ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1401 // Standard output streaming function
1403 // ostream &os output steam
1404 // AliITSTrackleterSPDEff &s class to be streamed.
1408 // ostream &os The stream pointer
1413 //______________________________________________________________________
1414 istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1415 // Standard inputput streaming function
1417 // istream &is input steam
1418 // AliITSTrackleterSPDEff &s class to be streamed.
1422 // ostream &os The stream pointer
1424 //printf("prova %d \n", (Int_t)s.GetMC());
1428 //______________________________________________________________________
1429 void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
1431 // This Method write into an either asci or root file
1432 // the used cuts and the statistics of the MC related quantities
1433 // The method SetMC() has to be called before
1434 // Input TString filename: name of file for output (it deletes already existing
1439 //if(!fMC) {CallWarningMC(); return;}
1440 if (!filename.Contains(".root")) {
1441 ofstream out(filename.Data(),ios::out | ios::binary);
1447 TFile* mcfile = TFile::Open(filename, "RECREATE");
1448 TH1F* cuts = new TH1F("cuts", "list of cuts", 10, 0, 10); // TH1I containing cuts
1449 cuts->SetBinContent(1,fPhiWindowL1);
1450 cuts->SetBinContent(2,fZetaWindowL1);
1451 cuts->SetBinContent(3,fPhiWindow);
1452 cuts->SetBinContent(4,fZetaWindow);
1453 cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
1454 cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
1455 cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
1456 cuts->SetBinContent(8,fReflectClusterAroundZAxisForLayer0);
1457 cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer1);
1458 cuts->SetBinContent(10,fMC);
1461 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
1463 TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
1464 mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
1465 mc0->SetBinContent(2,fUseOnlySecondaryForPred);
1466 mc0->SetBinContent(3,fUseOnlySameParticle);
1467 mc0->SetBinContent(4,fUseOnlyDifferentParticle);
1468 mc0->SetBinContent(5,fUseOnlyStableParticle);
1472 mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
1473 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
1475 mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
1476 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
1478 mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
1479 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
1481 mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
1482 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
1484 mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
1485 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
1487 mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
1488 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
1490 mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
1491 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
1493 mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
1494 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
1496 mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
1497 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
1499 mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
1500 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
1502 mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
1503 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
1505 mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
1506 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
1514 //____________________________________________________________________
1515 void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
1517 // This Method read from an asci file (do not know why binary does not work)
1518 // the cuts to be used and the statistics of the MC related quantities
1519 // Input TString filename: name of input file for output
1520 // The method SetMC() has to be called before
1524 //if(!fMC) {CallWarningMC(); return;}
1525 if( gSystem->AccessPathName( filename.Data() ) ) {
1526 AliError( Form( "file (%s) not found", filename.Data() ) );
1530 if (!filename.Contains(".root")) {
1531 ifstream in(filename.Data(),ios::in | ios::binary);
1538 TFile *mcfile = TFile::Open(filename);
1539 TH1F *cuts = (TH1F*)mcfile->Get("cuts");
1540 fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
1541 fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
1542 fPhiWindow=(Float_t)cuts->GetBinContent(3);
1543 fZetaWindow=(Float_t)cuts->GetBinContent(4);
1544 fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
1545 fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
1546 fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
1547 fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(8);
1548 fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(9);
1549 fMC=(Bool_t)cuts->GetBinContent(10);
1550 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1551 else { // only if file with MC predictions
1552 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1553 TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
1554 fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
1555 fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
1556 fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
1557 fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
1558 fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
1560 mc1 =(TH1I*)mcfile->Get("mc1");
1561 for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1562 mc1 =(TH1I*)mcfile->Get("mc2");
1563 for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1564 mc1 =(TH1I*)mcfile->Get("mc3");
1565 for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1566 mc1 =(TH1I*)mcfile->Get("mc4");
1567 for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1568 mc1 =(TH1I*)mcfile->Get("mc5");
1569 for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1570 mc1 =(TH1I*)mcfile->Get("mc6");
1571 for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
1572 mc1 =(TH1I*)mcfile->Get("mc7");
1573 for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1574 mc1 =(TH1I*)mcfile->Get("mc8");
1575 for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1576 mc1 =(TH1I*)mcfile->Get("mc9");
1577 for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1578 mc1 =(TH1I*)mcfile->Get("mc10");
1579 for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1580 mc1 =(TH1I*)mcfile->Get("mc11");
1581 for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1582 mc1 =(TH1I*)mcfile->Get("mc12");
1583 for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1589 //____________________________________________________________________
1590 Bool_t AliITSTrackleterSPDEff::SaveHists() {
1591 // This (private) method save the histograms on the output file
1592 // (only if fHistOn is TRUE).
1593 // Also the histograms from the base class are saved through the
1594 // AliITSMultReconstructor::SaveHists() call
1596 if (!GetHistOn()) return kFALSE;
1598 // AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1599 fhClustersDPhiAll->Write();
1600 fhClustersDThetaAll->Write();
1601 fhClustersDZetaAll->Write();
1602 fhDPhiVsDThetaAll->Write();
1603 fhDPhiVsDZetaAll->Write();
1605 fhClustersDPhiAcc->Write();
1606 fhClustersDThetaAcc->Write();
1607 fhClustersDZetaAcc->Write();
1608 fhDPhiVsDThetaAcc->Write();
1609 fhDPhiVsDZetaAcc->Write();
1611 fhetaTracklets->Write();
1612 fhphiTracklets->Write();
1613 fhetaClustersLay1->Write();
1614 fhphiClustersLay1->Write();
1616 fhClustersDPhiInterpAll->Write();
1617 fhClustersDThetaInterpAll->Write();
1618 fhClustersDZetaInterpAll->Write();
1619 fhDPhiVsDThetaInterpAll->Write();
1620 fhDPhiVsDZetaInterpAll->Write();
1622 fhClustersDPhiInterpAcc->Write();
1623 fhClustersDThetaInterpAcc->Write();
1624 fhClustersDZetaInterpAcc->Write();
1625 fhDPhiVsDThetaInterpAcc->Write();
1626 fhDPhiVsDZetaInterpAcc->Write();
1628 fhetaClustersLay2->Write();
1629 fhphiClustersLay2->Write();
1632 //__________________________________________________________
1633 Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1635 // Saves the histograms into a tree and saves the trees into a file
1636 // Also the histograms from the base class are saved
1638 if (!GetHistOn()) return kFALSE;
1639 if (!strcmp(filename.Data(),"")) {
1640 AliWarning("WriteHistosToFile: null output filename!");
1643 TFile *hFile=new TFile(filename.Data(),option,
1644 "The File containing the histos for SPD efficiency studies with tracklets");
1645 if(!SaveHists()) return kFALSE;
1650 //____________________________________________________________
1651 void AliITSTrackleterSPDEff::BookHistos() {
1653 // This method books addtitional histograms
1654 // w.r.t. those of the base class.
1655 // In particular, the differences of cluster coordinate between the two SPD
1656 // layers are computed in the interpolation phase
1658 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1660 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
1661 fhClustersDPhiAcc->SetDirectory(0);
1662 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
1663 fhClustersDThetaAcc->SetDirectory(0);
1664 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
1665 fhClustersDZetaAcc->SetDirectory(0);
1667 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
1668 fhDPhiVsDZetaAcc->SetDirectory(0);
1669 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
1670 fhDPhiVsDThetaAcc->SetDirectory(0);
1672 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
1673 fhClustersDPhiAll->SetDirectory(0);
1674 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
1675 fhClustersDThetaAll->SetDirectory(0);
1676 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
1677 fhClustersDZetaAll->SetDirectory(0);
1679 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
1680 fhDPhiVsDZetaAll->SetDirectory(0);
1681 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
1682 fhDPhiVsDThetaAll->SetDirectory(0);
1684 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
1685 fhetaTracklets->SetDirectory(0);
1686 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
1687 fhphiTracklets->SetDirectory(0);
1688 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
1689 fhetaClustersLay1->SetDirectory(0);
1690 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
1691 fhphiClustersLay1->SetDirectory(0);
1693 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1694 fhClustersDPhiInterpAcc->SetDirectory(0);
1695 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1696 fhClustersDThetaInterpAcc->SetDirectory(0);
1697 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1698 fhClustersDZetaInterpAcc->SetDirectory(0);
1700 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1701 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1702 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1703 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1705 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1706 fhClustersDPhiInterpAll->SetDirectory(0);
1707 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1708 fhClustersDThetaInterpAll->SetDirectory(0);
1709 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1710 fhClustersDZetaInterpAll->SetDirectory(0);
1712 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1713 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1714 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1715 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1717 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1718 fhetaClustersLay2->SetDirectory(0);
1719 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1720 fhphiClustersLay2->SetDirectory(0);
1723 //____________________________________________________________
1724 void AliITSTrackleterSPDEff::DeleteHistos() {
1726 // Private method to delete Histograms from memory
1727 // it is called. e.g., by the destructor.
1729 // form AliITSMultReconstructor
1730 if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
1731 if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
1732 if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
1733 if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
1734 if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
1735 if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
1736 if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
1737 if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
1738 if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
1739 if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
1740 if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
1741 if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
1742 if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
1743 if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
1745 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1746 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1747 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1748 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1749 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1750 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1751 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1752 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1753 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1754 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1755 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1756 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1758 //_______________________________________________________________
1759 Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1760 Float_t* vtx, AliStack *stack, TTree *ref) {
1761 // This (private) method can be used only for MC events, where both AliStack and the TrackReference
1763 // It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1765 // - Int_t layer (either 0 or 1): layer which you want to chech if the tracklete can be
1767 // - Int_t iC : cluster index used to build the tracklet prediction
1768 // if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1769 // - Float_t* vtx: actual event vertex
1770 // - stack: pointer to Stack
1771 // - ref: pointer to TTRee of TrackReference
1772 Bool_t ret=kFALSE; // returned value
1773 Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1774 if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1775 if(!stack) {AliError("null pointer to MC stack"); return ret;}
1776 if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1777 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1778 if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1780 AliTrackReference *tref=0x0;
1781 Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1782 Int_t nentries = (Int_t)ref->GetEntries();
1783 TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1784 TBranch *br = ref->GetBranch("TrackReferences");
1785 br->SetAddress(&tcaRef);
1786 for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1787 br->GetEntry(itrack);
1788 Int_t nref=tcaRef->GetEntriesFast();
1789 if(nref>0) { //it is enough to look at the first one
1790 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1791 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1794 if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1795 br->GetEntry(imatch); // redundant, nevertheless ...
1796 Int_t nref=tcaRef->GetEntriesFast();
1797 for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1798 tref=(AliTrackReference*)tcaRef->At(iref);
1799 if(tref->R()>10) continue; // not SPD ref
1800 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1801 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1803 // compute the proper quantities for this tref, as was done for fClustersLay1/2
1804 Float_t x = tref->X() - vtx[0];
1805 Float_t y = tref->Y() - vtx[1];
1806 Float_t z = tref->Z() - vtx[2];
1808 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1810 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1811 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1812 trefLayExtr[2] = z; // Store z
1814 if(layer==1) { // try to see if it is reconstructable at the outer layer
1815 // find the difference in angles
1816 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1817 // take into account boundary condition
1818 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1820 // find the difference in z (between linear projection from layer 1
1821 // and the actual point: Dzeta= z1/r1*r2 -z2)
1822 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1823 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1825 // make "elliptical" cut in Phi and Zeta!
1826 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
1827 dZeta*dZeta/fZetaWindow/fZetaWindow);
1828 if (d<1) {ret=kTRUE; break;}
1830 if(layer==0) { // try to see if it is reconstructable at the inner layer
1832 // find the difference in angles
1833 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1834 // take into account boundary condition
1835 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1837 // find the difference in z (between linear projection from layer 2
1838 // and the actual point: Dzeta= z2/r2*r1 -z1)
1839 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1840 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1842 // make "elliptical" cut in Phi and Zeta!
1843 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1844 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1845 if (d<1) {ret=kTRUE; break;};
1851 //_________________________________________________________________________
1852 void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
1854 // this method apply a rotation by 180 degree around the Z (beam) axis to all
1855 // the RecPoints in a given layer to be used to build tracklets.
1856 // **************** VERY IMPORTANT:: ***************
1857 // It must be called just after LoadClusterArrays, since afterwards the datamember
1858 // fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
1859 // instead of Cartesian
1861 if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1862 AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1864 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1865 fClustersLay1[iC1][0]*=-1;
1866 fClustersLay1[iC1][1]*=-1;
1870 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1871 fClustersLay2[iC2][0]*=-1;
1872 fClustersLay2[iC2][1]*=-1;
1877 //____________________________________________________________________________
1878 Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *){
1879 // This method is used to find the tracklets.
1880 // It is called from AliReconstruction
1881 // The vertex is supposed to be associated to the Tracker (i.e. to this) already
1882 // The cluster is supposed to be associated to the Tracker already
1883 // In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
1886 AliRunLoader* runLoader = AliRunLoader::Instance();
1888 Error("Clusters2Tracks", "no run loader found");
1891 AliStack *pStack=0x0; TTree *tRefTree=0x0;
1893 runLoader->LoadKinematics("read");
1894 runLoader->LoadTrackRefs("read");
1895 pStack= runLoader->Stack();
1896 tRefTree= runLoader->TreeTR();
1898 Reconstruct(pStack,tRefTree);
1901 //____________________________________________________________________________
1902 Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
1904 // It is called from AliReconstruction
1910 if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
1911 if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
1914 //____________________________________________________________________
1916 AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
1918 // - gets the clusters from the cluster tree
1919 // - convert them into global coordinates
1920 // - store them in the internal arrays
1921 // - count the number of cluster-fired chips
1923 //AliDebug(1,"Loading clusters and cluster-fired chips ...");
1928 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
1929 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
1931 itsClusterBranch->SetAddress(&itsClusters);
1933 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
1934 Float_t cluGlo[3]={0.,0.,0.};
1936 // loop over the its subdetectors
1937 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
1939 if (!itsClusterTree->GetEvent(iIts))
1942 Int_t nClusters = itsClusters->GetEntriesFast();
1944 // number of clusters in each chip of the current module
1947 // loop over clusters
1948 while(nClusters--) {
1949 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
1951 layer = cluster->GetLayer();
1952 if (layer>1) continue;
1954 cluster->GetGlobalXYZ(cluGlo);
1955 Float_t x = cluGlo[0];
1956 Float_t y = cluGlo[1];
1957 Float_t z = cluGlo[2];
1960 fClustersLay1[fNClustersLay1][0] = x;
1961 fClustersLay1[fNClustersLay1][1] = y;
1962 fClustersLay1[fNClustersLay1][2] = z;
1964 for (Int_t i=0; i<3; i++)
1965 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
1969 fClustersLay2[fNClustersLay2][0] = x;
1970 fClustersLay2[fNClustersLay2][1] = y;
1971 fClustersLay2[fNClustersLay2][2] = z;
1973 for (Int_t i=0; i<3; i++)
1974 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
1978 }// end of cluster loop
1980 } // end of its "subdetector" loop
1982 itsClusters->Delete();
1986 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));