]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSTrackleterSPDEff.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
ESD can init the Bfield via AliESDRun::InitMagneticField()
[u/mrichter/AliRoot.git] / ITS / AliITSTrackleterSPDEff.cxx
CommitLineData
275a301c 1/**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
275a301c 15//____________________________________________________________________
84161aec 16//
275a301c 17// AliITSTrackleterSPDEff - find SPD chips efficiencies by using tracklets.
84161aec 18//
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.
275a301c 22// The trackleter prediction is built using the vertex and 1 cluster.
23//
84161aec 24//
275a301c 25// Author : Giuseppe Eugenio Bruno, based on the skeleton of Reconstruct method provided by Tiziano Virgili
26// email: giuseppe.bruno@ba.infn.it
84161aec 27//
275a301c 28//____________________________________________________________________
29
84161aec 30/* $Id$ */
31
275a301c 32#include <TFile.h>
a3b31967 33#include <TTree.h>
275a301c 34#include <TParticle.h>
78178a72 35#include <TParticlePDG.h>
275a301c 36#include <TSystem.h>
37#include <Riostream.h>
a3b31967 38#include <TClonesArray.h>
275a301c 39
58e8dc31 40#include "AliTracker.h"
275a301c 41#include "AliITSTrackleterSPDEff.h"
42#include "AliITSgeomTGeo.h"
43#include "AliLog.h"
44#include "AliITSPlaneEffSPD.h"
45#include "AliStack.h"
a3b31967 46#include "AliTrackReference.h"
58e8dc31 47#include "AliRunLoader.h"
48#include "AliITSReconstructor.h"
49#include "AliITSRecPoint.h"
275a301c 50//____________________________________________________________________
51ClassImp(AliITSTrackleterSPDEff)
52
53
54//____________________________________________________________________
55AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
58e8dc31 56AliTracker(),
57//
58fClustersLay1(0),
59fClustersLay2(0),
60fTracklets(0),
61fAssociationFlag(0),
62fNClustersLay1(0),
63fNClustersLay2(0),
64fNTracklets(0),
65fOnlyOneTrackletPerC2(0),
7284b2b2 66fPhiWindowL2(0),
67fZetaWindowL2(0),
58e8dc31 68fPhiOverlapCut(0),
69fZetaOverlapCut(0),
70fHistOn(0),
71fhClustersDPhiAcc(0),
72fhClustersDThetaAcc(0),
73fhClustersDZetaAcc(0),
74fhClustersDPhiAll(0),
75fhClustersDThetaAll(0),
76fhClustersDZetaAll(0),
77fhDPhiVsDThetaAll(0),
78fhDPhiVsDThetaAcc(0),
79fhDPhiVsDZetaAll(0),
80fhDPhiVsDZetaAcc(0),
81fhetaTracklets(0),
82fhphiTracklets(0),
83fhetaClustersLay1(0),
84fhphiClustersLay1(0),
85//
275a301c 86fAssociationFlag1(0),
87fChipPredOnLay2(0),
88fChipPredOnLay1(0),
89fNTracklets1(0),
90fPhiWindowL1(0),
91fZetaWindowL1(0),
92fOnlyOneTrackletPerC1(0),
a3b31967 93fUpdateOncePerEventPlaneEff(0),
94fChipUpdatedInEvent(0),
275a301c 95fPlaneEffSPD(0),
03ee9629 96fPlaneEffBkg(0),
0fce916f 97fReflectClusterAroundZAxisForLayer0(kFALSE),
98fReflectClusterAroundZAxisForLayer1(kFALSE),
03ee9629 99fLightBkgStudyInParallel(kFALSE),
275a301c 100fMC(0),
95fcf795 101fUseOnlyPrimaryForPred(0),
275a301c 102fUseOnlySecondaryForPred(0),
103fUseOnlySameParticle(0),
104fUseOnlyDifferentParticle(0),
95fcf795 105fUseOnlyStableParticle(1),
275a301c 106fPredictionPrimary(0),
107fPredictionSecondary(0),
108fClusterPrimary(0),
109fClusterSecondary(0),
a3b31967 110fSuccessPP(0),
111fSuccessTT(0),
112fSuccessS(0),
113fSuccessP(0),
114fFailureS(0),
115fFailureP(0),
116fRecons(0),
117fNonRecons(0),
275a301c 118fhClustersDPhiInterpAcc(0),
119fhClustersDThetaInterpAcc(0),
120fhClustersDZetaInterpAcc(0),
121fhClustersDPhiInterpAll(0),
122fhClustersDThetaInterpAll(0),
123fhClustersDZetaInterpAll(0),
124fhDPhiVsDThetaInterpAll(0),
125fhDPhiVsDThetaInterpAcc(0),
126fhDPhiVsDZetaInterpAll(0),
127fhDPhiVsDZetaInterpAcc(0),
128fhetaClustersLay2(0),
17d531c2 129fhphiClustersLay2(0),
130fhClustersInChip(0),
131fhClustersInModuleLay1(0),
132fhClustersInModuleLay2(0)
275a301c 133{
84161aec 134 // default constructor
58e8dc31 135// from AliITSMultReconstructor
7284b2b2 136 SetPhiWindowL2();
137 SetZetaWindowL2();
58e8dc31 138 SetOnlyOneTrackletPerC2();
139 fClustersLay1 = new Float_t*[300000];
140 fClustersLay2 = new Float_t*[300000];
141 fTracklets = new Float_t*[300000];
142 fAssociationFlag = new Bool_t[300000];
143//
275a301c 144 SetPhiWindowL1();
145 SetZetaWindowL1();
146 SetOnlyOneTrackletPerC1();
147
148 fAssociationFlag1 = new Bool_t[300000];
149 fChipPredOnLay2 = new UInt_t[300000];
150 fChipPredOnLay1 = new UInt_t[300000];
a3b31967 151 fChipUpdatedInEvent = new Bool_t[1200];
275a301c 152
153 for(Int_t i=0; i<300000; i++) {
58e8dc31 154 // from AliITSMultReconstructor
155 fClustersLay1[i] = new Float_t[6];
156 fClustersLay2[i] = new Float_t[6];
157 fTracklets[i] = new Float_t[5];
158 fAssociationFlag[i] = kFALSE;
159 //
275a301c 160 fAssociationFlag1[i] = kFALSE;
161 }
a3b31967 162 for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
275a301c 163
164 if (GetHistOn()) BookHistos();
165
166 fPlaneEffSPD = new AliITSPlaneEffSPD();
03ee9629 167 SetLightBkgStudyInParallel();
275a301c 168}
169//______________________________________________________________________
58e8dc31 170AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
171AliTracker(mr),
172// from AliITSMultReconstructor
173fClustersLay1(mr.fClustersLay1),
174fClustersLay2(mr.fClustersLay2),
175fTracklets(mr.fTracklets),
176fAssociationFlag(mr.fAssociationFlag),
177fNClustersLay1(mr.fNClustersLay1),
178fNClustersLay2(mr.fNClustersLay2),
179fNTracklets(mr.fNTracklets),
180fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
7284b2b2 181fPhiWindowL2(mr.fPhiWindowL2),
182fZetaWindowL2(mr.fZetaWindowL2),
58e8dc31 183fPhiOverlapCut(mr.fPhiOverlapCut),
184fZetaOverlapCut(mr.fZetaOverlapCut),
185fHistOn(mr.fHistOn),
186fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
187fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
188fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
189fhClustersDPhiAll(mr.fhClustersDPhiAll),
190fhClustersDThetaAll(mr.fhClustersDThetaAll),
191fhClustersDZetaAll(mr.fhClustersDZetaAll),
192fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
193fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
194fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
195fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
196fhetaTracklets(mr.fhetaTracklets),
197fhphiTracklets(mr.fhphiTracklets),
198fhetaClustersLay1(mr.fhetaClustersLay1),
199fhphiClustersLay1(mr.fhphiClustersLay1),
200//
275a301c 201fAssociationFlag1(mr.fAssociationFlag1),
202fChipPredOnLay2(mr.fChipPredOnLay2),
203fChipPredOnLay1(mr.fChipPredOnLay1),
204fNTracklets1(mr.fNTracklets1),
205fPhiWindowL1(mr.fPhiWindowL1),
206fZetaWindowL1(mr.fZetaWindowL1),
207fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
a3b31967 208fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
209fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
275a301c 210fPlaneEffSPD(mr.fPlaneEffSPD),
03ee9629 211fPlaneEffBkg(mr.fPlaneEffBkg),
0fce916f 212fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
213fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
03ee9629 214fLightBkgStudyInParallel(mr.fLightBkgStudyInParallel),
275a301c 215fMC(mr.fMC),
216fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
217fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
218fUseOnlySameParticle(mr.fUseOnlySameParticle),
219fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
220fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
221fPredictionPrimary(mr.fPredictionPrimary),
222fPredictionSecondary(mr.fPredictionSecondary),
223fClusterPrimary(mr.fClusterPrimary),
224fClusterSecondary(mr.fClusterSecondary),
a3b31967 225fSuccessPP(mr.fSuccessPP),
226fSuccessTT(mr.fSuccessTT),
227fSuccessS(mr.fSuccessS),
228fSuccessP(mr.fSuccessP),
229fFailureS(mr.fFailureS),
230fFailureP(mr.fFailureP),
231fRecons(mr.fRecons),
232fNonRecons(mr.fNonRecons),
275a301c 233fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
234fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
235fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
236fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
237fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
238fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
239fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
240fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
241fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
242fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
243fhetaClustersLay2(mr.fhetaClustersLay2),
17d531c2 244fhphiClustersLay2(mr.fhphiClustersLay2),
245fhClustersInChip(mr.fhClustersInChip),
246fhClustersInModuleLay1(mr.fhClustersInModuleLay1),
247fhClustersInModuleLay2(mr.fhClustersInModuleLay2)
275a301c 248{
249 // Copy constructor
250}
251
252//______________________________________________________________________
253AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
254 // Assignment operator
255 this->~AliITSTrackleterSPDEff();
256 new(this) AliITSTrackleterSPDEff(mr);
257 return *this;
258}
259//______________________________________________________________________
260AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
261 // Destructor
58e8dc31 262// from AliITSMultReconstructor
263 // delete arrays
264 for(Int_t i=0; i<300000; i++) {
265 delete [] fClustersLay1[i];
266 delete [] fClustersLay2[i];
267 delete [] fTracklets[i];
268 }
269 delete [] fClustersLay1;
270 delete [] fClustersLay2;
271 delete [] fTracklets;
272 delete [] fAssociationFlag;
273//
275a301c 274 // delete histograms
275 DeleteHistos();
276
277 delete [] fAssociationFlag1;
278
279 delete [] fChipPredOnLay2;
280 delete [] fChipPredOnLay1;
281
a3b31967 282 delete [] fChipUpdatedInEvent;
283
275a301c 284 delete [] fPredictionPrimary;
285 delete [] fPredictionSecondary;
286 delete [] fClusterPrimary;
287 delete [] fClusterSecondary;
a3b31967 288 delete [] fSuccessPP;
289 delete [] fSuccessTT;
290 delete [] fSuccessS;
291 delete [] fSuccessP;
292 delete [] fFailureS;
293 delete [] fFailureP;
294 delete [] fRecons;
295 delete [] fNonRecons;
275a301c 296
297 // delete PlaneEff
298 delete fPlaneEffSPD;
03ee9629 299 fPlaneEffSPD=0;
300 if(fPlaneEffBkg) {
301 delete fPlaneEffBkg;
302 fPlaneEffBkg=0;
303
304 }
275a301c 305}
306//____________________________________________________________________
307void
03ee9629 308AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef, Bool_t lbkg) {
275a301c 309 //
58e8dc31 310 // - you have to take care of the following, before of using Reconstruct
311 // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
312 // and convert the cluster coordinates to theta, phi (seen from the
313 // interaction vertex).
314 // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
315 // - Find the extrapolation/interpolation point.
275a301c 316 // - Find the chip corresponding to that
317 // - Check if there is a cluster near that point
318 //
275a301c 319 // reset counters
03ee9629 320 if(lbkg && !GetLightBkgStudyInParallel()) {
321 AliError("You asked for lightBackground in the Reconstruction without proper call to SetLightBkgStudyInParallel(1)");
322 return;
323 }
324 AliITSPlaneEffSPD *pe;
325 if(lbkg) {
326 pe=fPlaneEffBkg;
327 } else {
328 pe=fPlaneEffSPD;
329 }
275a301c 330 fNTracklets = 0;
58e8dc31 331 // retrieve the vertex position
332 Float_t vtx[3];
333 vtx[0]=(Float_t)GetX();
334 vtx[1]=(Float_t)GetY();
335 vtx[2]=(Float_t)GetZ();
0fce916f 336 // to study residual background (i.e. contribution from TT' to measured efficiency)
03ee9629 337 if(fReflectClusterAroundZAxisForLayer0 && !lbkg) ReflectClusterAroundZAxisForLayer(0);
338 if(fReflectClusterAroundZAxisForLayer1 && !lbkg) ReflectClusterAroundZAxisForLayer(1);
0fce916f 339 //
03ee9629 340 if(fMC && !pStack && !lbkg) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
341 if(fMC && !tRef && !lbkg) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
275a301c 342 Bool_t found;
343 Int_t nfTraPred1=0; Int_t ntTraPred1=0;
344 Int_t nfTraPred2=0; Int_t ntTraPred2=0;
345 Int_t nfClu1=0; Int_t ntClu1=0;
346 Int_t nfClu2=0; Int_t ntClu2=0;
347
a3b31967 348 // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
349 // for this new event)
350 for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
275a301c 351
352 // find the tracklets
353 AliDebug(1,"Looking for tracklets... ");
354 AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
355
356 //###########################################################
357 // Loop on layer 1 : finding theta, phi and z
358 UInt_t key;
359 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
360 Float_t x = fClustersLay1[iC1][0] - vtx[0];
361 Float_t y = fClustersLay1[iC1][1] - vtx[1];
362 Float_t z = fClustersLay1[iC1][2] - vtx[2];
363
364 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
365
366 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
367 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
368 fClustersLay1[iC1][2] = z; // Store z
369
370 // find the Radius and the chip corresponding to the extrapolation point
371
372 found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
373 if (!found) {
374 AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
375 key=999999; // also some other actions should be taken if not Found
376 }
377 nfTraPred2+=(Int_t)found; // this for debugging purpose
378 ntTraPred2++; // to check efficiency of the method FindChip
379 fChipPredOnLay2[iC1] = key;
380 fAssociationFlag1[iC1] = kFALSE;
381
03ee9629 382 if (fHistOn && !lbkg) {
275a301c 383 Float_t eta=fClustersLay1[iC1][0];
384 eta= TMath::Tan(eta/2.);
385 eta=-TMath::Log(eta);
17d531c2 386 fhetaClustersLay1->Fill(eta);
275a301c 387 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
17d531c2 388 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
275a301c 389 }
390 }
275a301c 391 // Loop on layer 2 : finding theta, phi and r
392 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
393 Float_t x = fClustersLay2[iC2][0] - vtx[0];
394 Float_t y = fClustersLay2[iC2][1] - vtx[1];
395 Float_t z = fClustersLay2[iC2][2] - vtx[2];
396
397 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
398
399 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
400 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi (done properly in the range [0,2pi])
401 fClustersLay2[iC2][2] = z; // Store z
402
403 // find the Radius and the chip corresponding to the extrapolation point
404
405 found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
406 if (!found) {
407 AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
408 key=999999;
409 }
410 nfTraPred1+=(Int_t)found; // this for debugging purpose
411 ntTraPred1++; // to check efficiency of the method FindChip
412 fChipPredOnLay1[iC2] = key;
413 fAssociationFlag[iC2] = kFALSE;
414
03ee9629 415 if (fHistOn && !lbkg) {
275a301c 416 Float_t eta=fClustersLay2[iC2][0];
417 eta= TMath::Tan(eta/2.);
418 eta=-TMath::Log(eta);
419 fhetaClustersLay2->Fill(eta);
420 fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
17d531c2 421 fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
275a301c 422 }
423 }
424
425 //###########################################################
426
427 // First part : Extrapolation to Layer 2
428
429 // Loop on layer 1
430 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
431
a3b31967 432 // here the control to check whether the efficiency of the chip traversed by this tracklet
433 // prediction has already been updated in this event using another tracklet prediction
434 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
435
275a301c 436 // reset of variables for multiple candidates
437 Int_t iC2WithBestDist = 0; // reset
438 Float_t distmin = 100.; // just to put a huge number!
439 Float_t dPhimin = 0.; // Used for histograms only!
440 Float_t dThetamin = 0.; // Used for histograms only!
441 Float_t dZetamin = 0.; // Used for histograms only!
442
443 // in any case, if MC has been required, store statistics of primaries and secondaries
a3b31967 444 Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
03ee9629 445 if (fMC && !lbkg) {
275a301c 446 Int_t lab1=(Int_t)fClustersLay1[iC1][3];
447 Int_t lab2=(Int_t)fClustersLay1[iC1][4];
448 Int_t lab3=(Int_t)fClustersLay1[iC1][5];
449 // do it always as a function of the chip number used to built the prediction
450 found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
451 if (!found) {AliWarning(
452 Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
453 else {
454 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
455 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
456 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
457 { // this cluster is from a primary particle
458 fClusterPrimary[key]++;
a3b31967 459 primary=kTRUE;
275a301c 460 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
461 } else { // this cluster is from a secondary particle
462 fClusterSecondary[key]++;
a3b31967 463 secondary=kTRUE;
275a301c 464 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
465 }
466 }
467 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
468 // (in case the prediction is reliable)
469 if( fChipPredOnLay2[iC1]<1200) {
470 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
471 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
472 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
473 else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
a3b31967 474 if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
475 (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
476 (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
477 else fNonRecons[fChipPredOnLay2[iC1]]++;
275a301c 478 }
479 }
480
481 // Loop on layer 2
482 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
483
484 // The following excludes double associations
485 if (!fAssociationFlag[iC2]) {
486
487 // find the difference in angles
488 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
489 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
490 // take into account boundary condition
491 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
492
493 // find the difference in z (between linear projection from layer 1
494 // and the actual point: Dzeta= z1/r1*r2 -z2)
495 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
496 Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
497
03ee9629 498 if (fHistOn && !lbkg) {
275a301c 499 fhClustersDPhiAll->Fill(dPhi);
500 fhClustersDThetaAll->Fill(dTheta);
501 fhClustersDZetaAll->Fill(dZeta);
502 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
503 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
504 }
505
506 // make "elliptical" cut in Phi and Zeta!
7284b2b2 507 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
508 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
275a301c 509
510 if (d>1) continue;
511
512 //look for the minimum distance: the minimum is in iC2WithBestDist
513 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
514 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
515 dPhimin = dPhi;
516 dThetamin = dTheta;
517 dZetamin = dZeta;
518 iC2WithBestDist = iC2;
519 }
520 }
521 } // end of loop over clusters in layer 2
522
523 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
524
03ee9629 525 if (fHistOn && !lbkg) {
275a301c 526 fhClustersDPhiAcc->Fill(dPhimin);
527 fhClustersDThetaAcc->Fill(dThetamin);
528 fhClustersDZetaAcc->Fill(dZetamin);
529 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
530 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
531 }
532
533 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
534 // flag the association
535
536 // store the tracklet
537
538 // use the theta from the clusters in the first layer
539 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
540 // use the phi from the clusters in the first layer
541 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
542 // Store the difference between phi1 and phi2
543 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
544
545 // find labels
546 Int_t label1 = 0;
547 Int_t label2 = 0;
548 while (label2 < 3)
549 {
550 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
551 break;
552 label1++;
553 if (label1 == 3)
554 {
555 label1 = 0;
556 label2++;
557 }
558 }
559
560 if (label2 < 3)
561 {
562 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
563 }
564 else
565 {
566 fTracklets[fNTracklets][3] = -2;
567 }
568
03ee9629 569 if (fHistOn && !lbkg) {
275a301c 570 Float_t eta=fTracklets[fNTracklets][0];
571 eta= TMath::Tan(eta/2.);
572 eta=-TMath::Log(eta);
573 fhetaTracklets->Fill(eta);
574 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
575 }
576
577// Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
578 found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
579 if(!found){
580 AliWarning(
581 Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
582 key=999999;
583 }
584 nfClu2+=(Int_t)found; // this for debugging purpose
585 ntClu2++; // to check efficiency of the method FindChip
586 if(key<1200) { // the Chip has been found
03ee9629 587 if(fMC && !lbkg) { // this part only for MC
275a301c 588 // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
589 // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
590 // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
a3b31967 591 if (label2 < 3) {
592 fSuccessTT[key]++;
593 if(primary) fSuccessPP[key]++;
594 }
275a301c 595 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
596 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
597 }
598
599 if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
a3b31967 600 // OK, success
03ee9629 601 pe->UpDatePlaneEff(kTRUE,key); // success
a3b31967 602 fChipUpdatedInEvent[key]=kTRUE;
03ee9629 603 if(fMC && !lbkg) {
a3b31967 604 if(primary) fSuccessP[key]++;
605 if(secondary) fSuccessS[key]++;
606 }
275a301c 607 }
608 else {
03ee9629 609 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
a3b31967 610 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
03ee9629 611 if(fMC && !lbkg) {
a3b31967 612 if(primary) fSuccessP[key]++;
613 if(secondary) fSuccessS[key]++;
614 }
275a301c 615 }
616 }
617
618 fNTracklets++;
619
620 } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
a3b31967 621 else if (fChipPredOnLay2[iC1]<1200) {
03ee9629 622 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
a3b31967 623 fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
03ee9629 624 if(fMC && !lbkg) {
a3b31967 625 if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
626 if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
627 }
628 }
275a301c 629 } // end of loop over clusters in layer 1
630
631 fNTracklets1=fNTracklets;
632
633//###################################################################
634
635 // Second part : Interpolation to Layer 1
636
637 // Loop on layer 2
638 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
639
a3b31967 640 // here the control to check whether the efficiency of the chip traversed by this tracklet
641 // prediction has already been updated in this event using another tracklet prediction
642 if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
643
275a301c 644 // reset of variables for multiple candidates
645 Int_t iC1WithBestDist = 0; // reset
646 Float_t distmin = 100.; // just to put a huge number!
647 Float_t dPhimin = 0.; // Used for histograms only!
648 Float_t dThetamin = 0.; // Used for histograms only!
649 Float_t dZetamin = 0.; // Used for histograms only!
650
651 // in any case, if MC has been required, store statistics of primaries and secondaries
a3b31967 652 Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
03ee9629 653 if (fMC && !lbkg) {
275a301c 654 Int_t lab1=(Int_t)fClustersLay2[iC2][3];
655 Int_t lab2=(Int_t)fClustersLay2[iC2][4];
656 Int_t lab3=(Int_t)fClustersLay2[iC2][5];
657 // do it always as a function of the chip number used to built the prediction
658 found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
659 if (!found) {AliWarning(
660 Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
661 else {
662 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
663 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
664 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
665 { // this cluster is from a primary particle
666 fClusterPrimary[key]++;
a3b31967 667 primary=kTRUE;
275a301c 668 if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
669 } else { // this cluster is from a secondary particle
670 fClusterSecondary[key]++;
a3b31967 671 secondary=kTRUE;
275a301c 672 if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
673 }
674 }
675 // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
676 // (in case the prediction is reliable)
677 if( fChipPredOnLay1[iC2]<1200) {
678 if((lab1 != -2 && PrimaryTrackChecker(lab1,pStack) ) ||
679 (lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
680 (lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
681 else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
a3b31967 682 if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
683 (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
684 (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
685 else fNonRecons[fChipPredOnLay1[iC2]]++;
275a301c 686 }
687 }
688
689 // Loop on layer 1
690 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
691
692 // The following excludes double associations
693 if (!fAssociationFlag1[iC1]) {
694
695 // find the difference in angles
696 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
697 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
698 // take into account boundary condition
699 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
700
701
702 // find the difference in z (between linear projection from layer 2
703 // and the actual point: Dzeta= z2/r2*r1 -z1)
704 Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
705 Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
706
707
03ee9629 708 if (fHistOn && !lbkg) {
275a301c 709 fhClustersDPhiInterpAll->Fill(dPhi);
710 fhClustersDThetaInterpAll->Fill(dTheta);
711 fhClustersDZetaInterpAll->Fill(dZeta);
712 fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
713 fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
714 }
715 // make "elliptical" cut in Phi and Zeta!
716 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
717 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
718
719 if (d>1) continue;
720
721 //look for the minimum distance: the minimum is in iC1WithBestDist
722 if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
723 distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
724 dPhimin = dPhi;
725 dThetamin = dTheta;
726 dZetamin = dZeta;
727 iC1WithBestDist = iC1;
728 }
729 }
730 } // end of loop over clusters in layer 1
731
a3b31967 732 if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
275a301c 733
03ee9629 734 if (fHistOn && !lbkg) {
275a301c 735 fhClustersDPhiInterpAcc->Fill(dPhimin);
736 fhClustersDThetaInterpAcc->Fill(dThetamin);
737 fhClustersDZetaInterpAcc->Fill(dZetamin);
738 fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
739 fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
740 }
741
742 if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
743 // flag the association
744
745 // store the tracklet
746
747 // use the theta from the clusters in the first layer
748 fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
749 // use the phi from the clusters in the first layer
750 fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
751 // Store the difference between phi1 and phi2
752 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
753
754 // find labels
755 Int_t label1 = 0;
756 Int_t label2 = 0;
757 while (label2 < 3)
758 {
759 if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
760 break;
761 label1++;
762 if (label1 == 3)
763 {
764 label1 = 0;
765 label2++;
766 }
767 }
768
769 if (label2 < 3)
770 {
771 fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
772 }
773 else
774 {
775 fTracklets[fNTracklets][3] = -2;
776 }
777
778// Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
779 found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
780 if(!found){
781 AliWarning(
782 Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
783 key=999999;
784 }
785 nfClu1+=(Int_t)found; // this for debugging purpose
786 ntClu1++; // to check efficiency of the method FindChip
787 if(key<1200) {
03ee9629 788 if(fMC && !lbkg) { // this part only for MC
275a301c 789 // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
790 // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
791 // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
a3b31967 792 if (label2 < 3) { // same label
793 fSuccessTT[key]++;
794 if(primary) fSuccessPP[key]++;
795 }
275a301c 796 if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
797 if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
798 }
799
800 if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
a3b31967 801 // OK, success
03ee9629 802 pe->UpDatePlaneEff(kTRUE,key); // success
a3b31967 803 fChipUpdatedInEvent[key]=kTRUE;
03ee9629 804 if(fMC && !lbkg) {
a3b31967 805 if(primary) fSuccessP[key]++;
806 if(secondary) fSuccessS[key]++;
807 }
275a301c 808 } else {
03ee9629 809 pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
a3b31967 810 fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
03ee9629 811 if(fMC && !lbkg) {
a3b31967 812 if(primary) fSuccessP[key]++;
813 if(secondary) fSuccessS[key]++;
814 }
275a301c 815 }
816 }
817
818 fNTracklets++;
819
820 } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
a3b31967 821 else if (fChipPredOnLay1[iC2]<1200) {
03ee9629 822 pe->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
a3b31967 823 fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
03ee9629 824 if(fMC && !lbkg) {
a3b31967 825 if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
826 if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
827 }
828 }
275a301c 829 } // end of loop over clusters in layer 2
830
831 AliDebug(1,Form("%d tracklets found", fNTracklets));
832 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
833 AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
834 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
835 AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
836}
837//____________________________________________________________________
838Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
839 Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
840//
841// Input: a) layer number in the range [0,1]
842// b) vtx[3]: actual vertex
843// c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx
844// d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx
845// e) phiVtx /
846// Output: Unique key to locate a chip
847// return: kTRUE if succesfull
848
849 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
850 Double_t r=GetRLayer(layer);
851 //AliInfo(Form("Radius on layer %d is %f cm",layer,r));
852
853 // set phiVtx in the range [0,2pi]
854 if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
855
856 Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
857 // of the intersection of the tracklet with the pixel layer.
858 if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
859 else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
860 AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
861 Double_t vtxy[2]={vtx[0],vtx[1]};
862 if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
863 // this method gives you two interceptions
864 if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
865 // set phiAbs in the range [0,2pi]
866 if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
867 // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
868 // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
869 // taking the closest one to phiVtx
870 AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
871 } else phiAbs=phiVtx;
872 Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
873
874 // now you need to locate the chip within the idet detector,
875 // starting from the local coordinates in such a detector
876
877 Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
878 Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
879 if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
880
881 key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
882 return kTRUE;
883}
884//______________________________________________________________________________
885Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
a3b31967 886//
887// Return the average radius of a layer from Geometry
888//
275a301c 889 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
890 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
891
892 Double_t xyz[3], &x=xyz[0], &y=xyz[1];
893 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
894 Double_t r=TMath::Sqrt(x*x + y*y);
895
896 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
897 r += TMath::Sqrt(x*x + y*y);
898 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
899 r += TMath::Sqrt(x*x + y*y);
900 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
901 r += TMath::Sqrt(x*x + y*y);
902 r*=0.25;
903 return r;
904}
905//______________________________________________________________________________
906Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
907 Float_t &xloc, Float_t &zloc) {
908 // this method transform Global Cilindrical coordinates into local (i.e. module)
909 // cartesian coordinates
910 //
911 //Compute Cartesian Global Coordinate
912 Double_t xyzGlob[3],xyzLoc[3];
913 xyzGlob[2]=z;
914 xyzGlob[0]=r*TMath::Cos(phi);
915 xyzGlob[1]=r*TMath::Sin(phi);
916
917 xloc=0.;
918 zloc=0.;
919
920 if(idet<0) return kFALSE;
921
922 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
923 Int_t lad = Int_t(idet/ndet) + 1;
924 Int_t det = idet - (lad-1)*ndet + 1;
925
926 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
927
928 xloc = (Float_t)xyzLoc[0];
929 zloc = (Float_t)xyzLoc[2];
930
931return kTRUE;
932}
933//______________________________________________________________________________
934Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
935 //--------------------------------------------------------------------
a3b31967 936 // This function finds the detector crossed by the track
937 // Input: layer in range [0,1]
938 // phi in ALICE absolute reference system
939 // z " " " " "
275a301c 940 //--------------------------------------------------------------------
941 if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
942 Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
943 Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
944 Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
945
946 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
947 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
948 Double_t phiOffset=TMath::ATan2(y,x);
949 Double_t zOffset=z2;
950
951 Double_t dphi;
952 if (zOffset<0) // old geometry
953 dphi = -(phi-phiOffset);
954 else // new geometry
955 dphi = phi-phiOffset;
956
957 if (dphi < 0) dphi += 2*TMath::Pi();
958 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
959 Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
960 if (np>=nladders) np-=nladders;
961 if (np<0) np+=nladders;
962
963 Double_t dz=zOffset-z;
964 Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
965 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
966 if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;}
967 if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
968
969 return np*ndetectors + nz;
970}
971//____________________________________________________________
972Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
973// this method find the intersection in xy between a tracklet (straight line) and
974// a circonference (r=R), using polar coordinates.
975/*
976Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
977 - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
978 - R: radius of the circle
979Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
980
981Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
982r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
983In the same system, the equation of a semi-line is: phi=phiVtx;
984Hence you get one interception only: P=(r,phiVtx)
a3b31967 985Finally you want P in the ABSOLUTE ALICE reference system.
275a301c 986*/
987Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
988Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
989Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
990Double_t cC=rO*rO-R*R;
991Double_t dDelta=bB*bB-4*cC;
992if(dDelta<0) return kFALSE;
993Double_t r1,r2;
994r1=(-bB-TMath::Sqrt(dDelta))/2;
995r2=(-bB+TMath::Sqrt(dDelta))/2;
996if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
997Double_t r=TMath::Max(r1,r2); // take the positive
998Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
999Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
1000pvtx[0]=r*TMath::Cos(phiVtx);
1001pvtx[1]=r*TMath::Sin(phiVtx);
1002pP[0]=vtx[0]+pvtx[0];
1003pP[1]=vtx[1]+pvtx[1];
1004phi=TMath::ATan2(pP[1],pP[0]);
1005return kTRUE;
1006}
1007//___________________________________________________________
1008Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
a3b31967 1009//
1010// simple method to reduce all angles (in rad)
1011// in range [0,2pi[
1012//
1013//
275a301c 1014while(angle >=2*TMath::Pi() || angle<0) {
1015 if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
1016 if(angle < 0) angle+=2*TMath::Pi();
1017}
1018return kTRUE;
1019}
1020//___________________________________________________________
1021Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
a3b31967 1022//
1023// This method check if a particle is primary; i.e.
1024// it comes from the main vertex and it is a "stable" particle, according to
1025// AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
1026// a stable particle: it has no effect on this analysis).
1027// This method can be called only for MC events, where Kinematics is available.
58e8dc31 1028// if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
a3b31967 1029// returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
1030// The latter (see below) try to verify if a primary particle is also "detectable".
1031//
275a301c 1032if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
1033if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
1034if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
1035// return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
1036 if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
1037 // like below: as in the correction for Multiplicity (i.e. by hand in macro)
1038 TParticle* part = stack->Particle(ipart);
1039 TParticle* part0 = stack->Particle(0); // first primary
1040 TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
1041 if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
1042 part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
1043
1044 if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
1045 TParticlePDG* pdgPart = part->GetPDG();
1046 if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
1047
1048 Double_t distx = part->Vx() - part0->Vx();
1049 Double_t disty = part->Vy() - part0->Vy();
1050 Double_t distz = part->Vz() - part0->Vz();
1051 Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
1052
1053 if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
1054 part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
1055 return kFALSE; }// primary if within 500 microns from true Vertex
1056
a3b31967 1057 if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
275a301c 1058 return kTRUE;
1059}
1060//_____________________________________________________________________________________________
1061Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
a3b31967 1062//
1063// This private method can be applied on MC particles (if stack is available),
1064// provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
1065//
1066// It define "detectable" a primary particle according to the following criteria:
1067//
1068// - if no decay products can be found in the stack (note that this does not
1069// means it is stable, since a particle is stored in stack if it has at least 1 hit in a
1070// sensitive detector)
1071// - if it has at least one decay daughter produced outside or just on the outer pixel layer
1072// - if the last decay particle is an electron (or a muon) which is not produced in-between
1073// the two pixel layers (this is likely to be a kink).
275a301c 1074if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
1075if(!stack) {AliError("null pointer to MC stack"); return 0;}
1076if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
1077
1078TParticle* part = stack->Particle(ipart);
1079//TParticle* part0 = stack->Particle(0); // first primary
1080
1081 Int_t nret=0;
1082 TParticle* dau = 0;
1083 Int_t nDau = 0;
a3b31967 1084 Int_t pdgDau;
1085 Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
1086 // its real fate ! But you have to take it !
1087 if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
275a301c 1088 Int_t lastDau = part->GetLastDaughter();
1089 nDau = lastDau - firstDau + 1;
a3b31967 1090 Double_t distMax=0.;
1091 Int_t jmax=0;
1092 for(Int_t j=firstDau; j<=lastDau; j++) {
1093 dau = stack->Particle(j);
1094 Double_t distx = dau->Vx();
1095 Double_t disty = dau->Vy();
1096 //Double_t distz = dau->Vz();
1097 Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
1098 if(distR<distMax) continue; // considere only the daughter produced at largest radius
1099 distMax=distR;
1100 jmax=j;
1101 }
1102 dau = stack->Particle(jmax);
1103 pdgDau=dau->GetPdgCode();
1104 if (pdgDau == 11 || pdgDau == 13 ) {
1105 if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
1106 else nret =0; // delta-ray emission in material (keep it)
1107 }
1108 else {// not ele or muon
1109 if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
275a301c 1110 }
a3b31967 1111return nret;
275a301c 1112}
1113//_________________________________________________________________
1114void AliITSTrackleterSPDEff::InitPredictionMC() {
a3b31967 1115//
1116// this method allocate memory for the MC related informations
1117// all the counters are set to 0
1118//
1119//
275a301c 1120if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
1121fPredictionPrimary = new Int_t[1200];
1122fPredictionSecondary = new Int_t[1200];
1123fClusterPrimary = new Int_t[1200];
1124fClusterSecondary = new Int_t[1200];
a3b31967 1125fSuccessPP = new Int_t[1200];
1126fSuccessTT = new Int_t[1200];
1127fSuccessS = new Int_t[1200];
1128fSuccessP = new Int_t[1200];
1129fFailureS = new Int_t[1200];
1130fFailureP = new Int_t[1200];
1131fRecons = new Int_t[1200];
1132fNonRecons = new Int_t[1200];
275a301c 1133for(Int_t i=0; i<1200; i++) {
1134 fPredictionPrimary[i]=0;
1135 fPredictionSecondary[i]=0;
1136 fPredictionSecondary[i]=0;
1137 fClusterSecondary[i]=0;
a3b31967 1138 fSuccessPP[i]=0;
1139 fSuccessTT[i]=0;
1140 fSuccessS[i]=0;
1141 fSuccessP[i]=0;
1142 fFailureS[i]=0;
1143 fFailureP[i]=0;
1144 fRecons[i]=0;
1145 fNonRecons[i]=0;
275a301c 1146}
1147return;
1148}
c6a05d92 1149//_________________________________________________________________
1150void AliITSTrackleterSPDEff::DeletePredictionMC() {
1151//
1152// this method deallocate memory for the MC related informations
1153// all the counters are set to 0
1154//
1155//
1156if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
1157if(fPredictionPrimary) {
1158 delete fPredictionPrimary; fPredictionPrimary=0;
1159}
1160if(fPredictionSecondary) {
1161 delete fPredictionSecondary; fPredictionSecondary=0;
1162}
1163if(fClusterPrimary) {
1164 delete fClusterPrimary; fClusterPrimary=0;
1165}
1166if(fClusterSecondary) {
1167 delete fClusterSecondary; fClusterSecondary=0;
1168}
1169if(fSuccessPP) {
1170 delete fSuccessPP; fSuccessPP=0;
1171}
1172if(fSuccessTT) {
1173 delete fSuccessTT; fSuccessTT=0;
1174}
1175if(fSuccessS) {
1176 delete fSuccessS; fSuccessS=0;
1177}
1178if(fSuccessP) {
1179 delete fSuccessP; fSuccessP=0;
1180}
1181if(fFailureS) {
1182 delete fFailureS; fFailureS=0;
1183}
1184if(fFailureP) {
1185 delete fFailureP; fFailureP=0;
1186}
1187if(fRecons) {
1188 delete fRecons; fRecons=0;
1189}
1190if(fNonRecons) {
1191 delete fNonRecons; fNonRecons=0;
1192}
1193return;
1194}
275a301c 1195//______________________________________________________________________
1196Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
84161aec 1197//
1198// This method return the Data menmber fPredictionPrimary [1200].
1199// You can call it only for MC events.
1200// fPredictionPrimary[key] contains the number of tracklet predictions on the
1201// given chip key built using a cluster on the other layer produced (at least)
1202// from a primary particle.
1203// Key refers to the chip crossed by the prediction
1204//
1205//
275a301c 1206if (!fMC) {CallWarningMC(); return 0;}
1207if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1208return fPredictionPrimary[(Int_t)key];
1209}
1210//______________________________________________________________________
1211Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
84161aec 1212//
1213// This method return the Data menmber fPredictionSecondary [1200].
1214// You can call it only for MC events.
1215// fPredictionSecondary[key] contains the number of tracklet predictions on the
1216// given chip key built using a cluster on the other layer produced (only)
1217// from a secondary particle
1218// Key refers to the chip crossed by the prediction
1219//
1220//
275a301c 1221if (!fMC) {CallWarningMC(); return 0;}
1222if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1223return fPredictionSecondary[(Int_t)key];
1224}
1225//______________________________________________________________________
1226Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
84161aec 1227//
1228// This method return the Data menmber fClusterPrimary [1200].
1229// You can call it only for MC events.
1230// fClusterPrimary[key] contains the number of tracklet predictions
1231// built using a cluster on that layer produced (only)
1232// from a primary particle
1233// Key refers to the chip used to build the prediction
1234//
1235//
275a301c 1236if (!fMC) {CallWarningMC(); return 0;}
1237if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1238return fClusterPrimary[(Int_t)key];
1239}
1240//______________________________________________________________________
1241Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
84161aec 1242//
1243// This method return the Data menmber fClusterSecondary [1200].
1244// You can call it only for MC events.
1245// fClusterSecondary[key] contains the number of tracklet predictions
1246// built using a cluster on that layer produced (only)
1247// from a secondary particle
1248// Key refers to the chip used to build the prediction
1249//
275a301c 1250if (!fMC) {CallWarningMC(); return 0;}
1251if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1252return fClusterSecondary[(Int_t)key];
1253}
1254//______________________________________________________________________
a3b31967 1255Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
1256//
1257// This method return the Data menmber fSuccessPP [1200].
1258// You can call it only for MC events.
1259// fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
1260// with a cluster on the other layer) built by using the same primary particle
1261// the unique chip key refers to the chip which get updated its efficiency
1262//
1263if (!fMC) {CallWarningMC(); return 0;}
1264if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1265return fSuccessPP[(Int_t)key];
1266}
1267//______________________________________________________________________
1268Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
1269//
1270// This method return the Data menmber fSuccessTT [1200].
1271// You can call it only for MC events.
1272// fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
1273// with a cluster on the other layer) built by using the same particle (whatever)
1274// the unique chip key refers to the chip which get updated its efficiency
1275//
1276if (!fMC) {CallWarningMC(); return 0;}
1277if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1278return fSuccessTT[(Int_t)key];
1279}
1280//______________________________________________________________________
1281Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
1282//
1283// This method return the Data menmber fSuccessS [1200].
1284// You can call it only for MC events.
1285// fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
1286// with a cluster on the other layer) built by using a secondary particle
1287// the unique chip key refers to the chip which get updated its efficiency
1288//
1289if (!fMC) {CallWarningMC(); return 0;}
1290if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1291return fSuccessS[(Int_t)key];
1292}
1293//______________________________________________________________________
1294Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
1295//
1296// This method return the Data menmber fSuccessP [1200].
1297// You can call it only for MC events.
1298// fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
1299// with a cluster on the other layer) built by using a primary particle
1300// the unique chip key refers to the chip which get updated its efficiency
1301//
1302if (!fMC) {CallWarningMC(); return 0;}
1303if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1304return fSuccessP[(Int_t)key];
1305}
1306//______________________________________________________________________
1307Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
1308//
1309// This method return the Data menmber fFailureS [1200].
1310// You can call it only for MC events.
1311// fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
1312// with a cluster on the other layer) built by using a secondary particle
1313// the unique chip key refers to the chip which get updated its efficiency
1314//
1315if (!fMC) {CallWarningMC(); return 0;}
1316if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1317return fFailureS[(Int_t)key];
1318}
1319//______________________________________________________________________
1320Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
1321//
1322// This method return the Data menmber fFailureP [1200].
1323// You can call it only for MC events.
1324// fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
1325// with a cluster on the other layer) built by using a primary particle
1326// the unique chip key refers to the chip which get updated its efficiency
1327//
1328if (!fMC) {CallWarningMC(); return 0;}
1329if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1330return fFailureP[(Int_t)key];
1331}
1332//_____________________________________________________________________
1333Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
1334//
1335// This method return the Data menmber fRecons [1200].
1336// You can call it only for MC events.
1337// fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
1338// has an hit in the detector)
1339// the unique chip key refers to the chip where fall the prediction
1340//
1341if (!fMC) {CallWarningMC(); return 0;}
1342if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1343return fRecons[(Int_t)key];
1344}
1345//_____________________________________________________________________
1346Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
1347//
1348// This method return the Data menmber fNonRecons [1200].
1349// You can call it only for MC events.
1350// fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
1351// has not any hit in the detector)
1352// the unique chip key refers to the chip where fall the prediction
1353//
1354if (!fMC) {CallWarningMC(); return 0;}
1355if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
1356return fNonRecons[(Int_t)key];
1357}
1358//______________________________________________________________________
275a301c 1359void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
1360 // Print out some class data values in Ascii Form to output stream
1361 // Inputs:
1362 // ostream *os Output stream where Ascii data is to be writen
1363 // Outputs:
1364 // none.
1365 // Return:
1366 // none.
7284b2b2 1367 *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindowL2 <<" "<< fZetaWindowL2
a3b31967 1368 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
0fce916f 1369 << " " << fUpdateOncePerEventPlaneEff << " " << fReflectClusterAroundZAxisForLayer0
1370 << " " << fReflectClusterAroundZAxisForLayer1;
275a301c 1371 *os << " " << fMC;
1372 if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
1373 *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
1374 << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle
1375 << " " << fUseOnlyStableParticle ;
1376 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ;
1377 for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
1378 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
1379 for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
a3b31967 1380 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
1381 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
1382 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
1383 for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
1384 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
1385 for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
1386 for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
1387 for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
275a301c 1388 return;
1389}
1390//______________________________________________________________________
1391void AliITSTrackleterSPDEff::ReadAscii(istream *is){
1392 // Read in some class data values in Ascii Form to output stream
1393 // Inputs:
1394 // istream *is Input stream where Ascii data is to be read in from
1395 // Outputs:
1396 // none.
1397 // Return:
1398 // none.
1399
c6a05d92 1400 Bool_t tmp= fMC;
7284b2b2 1401 *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindowL2 >> fZetaWindowL2
a3b31967 1402 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
0fce916f 1403 >> fUpdateOncePerEventPlaneEff >> fReflectClusterAroundZAxisForLayer0
1404 >> fReflectClusterAroundZAxisForLayer1;
c6a05d92 1405 //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
275a301c 1406 *is >> fMC;
c6a05d92 1407 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1408 else {
1409 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1410 *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
1411 >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
1412 >> fUseOnlyStableParticle;
1413 for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
1414 for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
1415 for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
1416 for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
1417 for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
1418 for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
1419 for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
1420 for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
1421 for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
1422 for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
1423 for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
1424 for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
1425 }
275a301c 1426 return;
1427}
1428//______________________________________________________________________
1429ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
1430 // Standard output streaming function
1431 // Inputs:
1432 // ostream &os output steam
1433 // AliITSTrackleterSPDEff &s class to be streamed.
1434 // Output:
1435 // none.
1436 // Return:
1437 // ostream &os The stream pointer
1438
1439 s.PrintAscii(&os);
1440 return os;
1441}
1442//______________________________________________________________________
1443istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
1444 // Standard inputput streaming function
1445 // Inputs:
1446 // istream &is input steam
1447 // AliITSTrackleterSPDEff &s class to be streamed.
1448 // Output:
1449 // none.
1450 // Return:
1451 // ostream &os The stream pointer
1452
1453 //printf("prova %d \n", (Int_t)s.GetMC());
1454 s.ReadAscii(&is);
1455 return is;
1456}
1457//______________________________________________________________________
1458void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
84161aec 1459//
58e8dc31 1460// This Method write into an either asci or root file
84161aec 1461// the used cuts and the statistics of the MC related quantities
1462// The method SetMC() has to be called before
1463// Input TString filename: name of file for output (it deletes already existing
1464// file)
1465// Output: none
1466//
1467//
c6a05d92 1468 //if(!fMC) {CallWarningMC(); return;}
1469 if (!filename.Contains(".root")) {
1470 ofstream out(filename.Data(),ios::out | ios::binary);
1471 out << *this;
1472 out.close();
1473 return;
1474 }
1475 else {
1476 TFile* mcfile = TFile::Open(filename, "RECREATE");
1477 TH1F* cuts = new TH1F("cuts", "list of cuts", 10, 0, 10); // TH1I containing cuts
1478 cuts->SetBinContent(1,fPhiWindowL1);
1479 cuts->SetBinContent(2,fZetaWindowL1);
7284b2b2 1480 cuts->SetBinContent(3,fPhiWindowL2);
1481 cuts->SetBinContent(4,fZetaWindowL2);
c6a05d92 1482 cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
1483 cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
1484 cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
1485 cuts->SetBinContent(8,fReflectClusterAroundZAxisForLayer0);
1486 cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer1);
1487 cuts->SetBinContent(10,fMC);
1488 cuts->Write();
1489 delete cuts;
1490 if(!fMC) {AliInfo("Writing only cuts, no MC info");}
1491 else {
1492 TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
1493 mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
1494 mc0->SetBinContent(2,fUseOnlySecondaryForPred);
1495 mc0->SetBinContent(3,fUseOnlySameParticle);
1496 mc0->SetBinContent(4,fUseOnlyDifferentParticle);
1497 mc0->SetBinContent(5,fUseOnlyStableParticle);
1498 mc0->Write();
1499 delete mc0;
1500 TH1I *mc1;
1501 mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
1502 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
1503 mc1->Write();
eb6a1ac8 1504 delete mc1;
c6a05d92 1505 mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
1506 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
1507 mc1->Write();
eb6a1ac8 1508 delete mc1;
c6a05d92 1509 mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
1510 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
1511 mc1->Write();
eb6a1ac8 1512 delete mc1;
c6a05d92 1513 mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
1514 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
1515 mc1->Write();
eb6a1ac8 1516 delete mc1;
c6a05d92 1517 mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
1518 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
1519 mc1->Write();
eb6a1ac8 1520 delete mc1;
c6a05d92 1521 mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
1522 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
1523 mc1->Write();
eb6a1ac8 1524 delete mc1;
c6a05d92 1525 mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
1526 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
1527 mc1->Write();
eb6a1ac8 1528 delete mc1;
c6a05d92 1529 mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
1530 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
1531 mc1->Write();
eb6a1ac8 1532 delete mc1;
c6a05d92 1533 mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
1534 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
1535 mc1->Write();
eb6a1ac8 1536 delete mc1;
c6a05d92 1537 mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
1538 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
1539 mc1->Write();
eb6a1ac8 1540 delete mc1;
c6a05d92 1541 mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
1542 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
1543 mc1->Write();
eb6a1ac8 1544 delete mc1;
c6a05d92 1545 mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
1546 for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
1547 mc1->Write();
1548 delete mc1;
1549 }
1550 mcfile->Close();
1551 }
275a301c 1552return;
1553}
1554//____________________________________________________________________
1555void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
84161aec 1556//
1557// This Method read from an asci file (do not know why binary does not work)
1558// the cuts to be used and the statistics of the MC related quantities
1559// Input TString filename: name of input file for output
1560// The method SetMC() has to be called before
1561// Output: none
1562//
1563//
c6a05d92 1564 //if(!fMC) {CallWarningMC(); return;}
275a301c 1565 if( gSystem->AccessPathName( filename.Data() ) ) {
1566 AliError( Form( "file (%s) not found", filename.Data() ) );
1567 return;
1568 }
1569
c6a05d92 1570 if (!filename.Contains(".root")) {
1571 ifstream in(filename.Data(),ios::in | ios::binary);
1572 in >> *this;
1573 in.close();
1574 return;
1575 }
1576 else {
1577 Bool_t tmp= fMC;
1578 TFile *mcfile = TFile::Open(filename);
1579 TH1F *cuts = (TH1F*)mcfile->Get("cuts");
1580 fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
1581 fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
7284b2b2 1582 fPhiWindowL2=(Float_t)cuts->GetBinContent(3);
1583 fZetaWindowL2=(Float_t)cuts->GetBinContent(4);
c6a05d92 1584 fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
1585 fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
1586 fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
1587 fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(8);
1588 fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(9);
1589 fMC=(Bool_t)cuts->GetBinContent(10);
1590 if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
1591 else { // only if file with MC predictions
1592 if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
1593 TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
1594 fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
1595 fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
1596 fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
1597 fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
1598 fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
1599 TH1I *mc1;
1600 mc1 =(TH1I*)mcfile->Get("mc1");
1601 for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1602 mc1 =(TH1I*)mcfile->Get("mc2");
1603 for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1604 mc1 =(TH1I*)mcfile->Get("mc3");
1605 for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1606 mc1 =(TH1I*)mcfile->Get("mc4");
1607 for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
1608 mc1 =(TH1I*)mcfile->Get("mc5");
1609 for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1610 mc1 =(TH1I*)mcfile->Get("mc6");
1611 for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
1612 mc1 =(TH1I*)mcfile->Get("mc7");
1613 for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1614 mc1 =(TH1I*)mcfile->Get("mc8");
1615 for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1616 mc1 =(TH1I*)mcfile->Get("mc9");
1617 for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
1618 mc1 =(TH1I*)mcfile->Get("mc10");
1619 for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
1620 mc1 =(TH1I*)mcfile->Get("mc11");
1621 for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1622 mc1 =(TH1I*)mcfile->Get("mc12");
1623 for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
1624 }
1625 mcfile->Close();
1626 }
275a301c 1627 return;
1628}
1629//____________________________________________________________________
1630Bool_t AliITSTrackleterSPDEff::SaveHists() {
84161aec 1631 // This (private) method save the histograms on the output file
275a301c 1632 // (only if fHistOn is TRUE).
84161aec 1633 // Also the histograms from the base class are saved through the
1634 // AliITSMultReconstructor::SaveHists() call
275a301c 1635
1636 if (!GetHistOn()) return kFALSE;
1637
58e8dc31 1638// AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
1639 fhClustersDPhiAll->Write();
1640 fhClustersDThetaAll->Write();
1641 fhClustersDZetaAll->Write();
1642 fhDPhiVsDThetaAll->Write();
1643 fhDPhiVsDZetaAll->Write();
1644
1645 fhClustersDPhiAcc->Write();
1646 fhClustersDThetaAcc->Write();
1647 fhClustersDZetaAcc->Write();
1648 fhDPhiVsDThetaAcc->Write();
1649 fhDPhiVsDZetaAcc->Write();
1650
1651 fhetaTracklets->Write();
1652 fhphiTracklets->Write();
1653 fhetaClustersLay1->Write();
1654 fhphiClustersLay1->Write();
275a301c 1655
1656 fhClustersDPhiInterpAll->Write();
1657 fhClustersDThetaInterpAll->Write();
1658 fhClustersDZetaInterpAll->Write();
1659 fhDPhiVsDThetaInterpAll->Write();
1660 fhDPhiVsDZetaInterpAll->Write();
1661
1662 fhClustersDPhiInterpAcc->Write();
1663 fhClustersDThetaInterpAcc->Write();
1664 fhClustersDZetaInterpAcc->Write();
1665 fhDPhiVsDThetaInterpAcc->Write();
1666 fhDPhiVsDZetaInterpAcc->Write();
1667
1668 fhetaClustersLay2->Write();
1669 fhphiClustersLay2->Write();
17d531c2 1670 fhClustersInChip->Write();
1671 for (Int_t nhist=0;nhist<80;nhist++){
1672 fhClustersInModuleLay1[nhist]->Write();
1673 }
1674 for (Int_t nhist=0;nhist<160;nhist++){
1675 fhClustersInModuleLay2[nhist]->Write();
1676 }
275a301c 1677 return kTRUE;
1678}
1679//__________________________________________________________
1680Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1681 //
1682 // Saves the histograms into a tree and saves the trees into a file
84161aec 1683 // Also the histograms from the base class are saved
275a301c 1684 //
1685 if (!GetHistOn()) return kFALSE;
5bd7ec3a 1686 if (!strcmp(filename.Data(),"")) {
275a301c 1687 AliWarning("WriteHistosToFile: null output filename!");
1688 return kFALSE;
1689 }
1690 TFile *hFile=new TFile(filename.Data(),option,
1691 "The File containing the histos for SPD efficiency studies with tracklets");
1692 if(!SaveHists()) return kFALSE;
1693 hFile->Write();
1694 hFile->Close();
1695 return kTRUE;
1696}
1697//____________________________________________________________
1698void AliITSTrackleterSPDEff::BookHistos() {
84161aec 1699//
1700// This method books addtitional histograms
1701// w.r.t. those of the base class.
1702// In particular, the differences of cluster coordinate between the two SPD
1703// layers are computed in the interpolation phase
1704//
275a301c 1705 if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
58e8dc31 1706//
1707 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
1708 fhClustersDPhiAcc->SetDirectory(0);
1709 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
1710 fhClustersDThetaAcc->SetDirectory(0);
1711 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
1712 fhClustersDZetaAcc->SetDirectory(0);
1713
1714 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
1715 fhDPhiVsDZetaAcc->SetDirectory(0);
1716 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
1717 fhDPhiVsDThetaAcc->SetDirectory(0);
1718
1719 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
1720 fhClustersDPhiAll->SetDirectory(0);
1721 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
1722 fhClustersDThetaAll->SetDirectory(0);
1723 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
1724 fhClustersDZetaAll->SetDirectory(0);
1725
1726 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
1727 fhDPhiVsDZetaAll->SetDirectory(0);
1728 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
1729 fhDPhiVsDThetaAll->SetDirectory(0);
1730
1731 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
1732 fhetaTracklets->SetDirectory(0);
1733 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
1734 fhphiTracklets->SetDirectory(0);
1735 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
1736 fhetaClustersLay1->SetDirectory(0);
1737 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
1738 fhphiClustersLay1->SetDirectory(0);
1739//
275a301c 1740 fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
1741 fhClustersDPhiInterpAcc->SetDirectory(0);
1742 fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1743 fhClustersDThetaInterpAcc->SetDirectory(0);
1744 fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1745 fhClustersDZetaInterpAcc->SetDirectory(0);
1746
1747 fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1748 fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1749 fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1750 fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1751
1752 fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5);
1753 fhClustersDPhiInterpAll->SetDirectory(0);
1754 fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1755 fhClustersDThetaInterpAll->SetDirectory(0);
1756 fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1757 fhClustersDZetaInterpAll->SetDirectory(0);
1758
1759 fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1760 fhDPhiVsDZetaInterpAll->SetDirectory(0);
1761 fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1762 fhDPhiVsDThetaInterpAll->SetDirectory(0);
1763
1764 fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.);
1765 fhetaClustersLay2->SetDirectory(0);
1766 fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1767 fhphiClustersLay2->SetDirectory(0);
17d531c2 1768 fhClustersInChip = new TH1F("fhClustersInChip", "ClustersPerChip", 1200, -0.5, 1199.5);
1769 fhClustersInChip->SetDirectory(0);
1770// each chip is divided 8(z) x 4(y), i.e. in 32 squares, each containing 4 columns and 64 rows.
1771 Float_t bz[160]; const Float_t kconv = 1.0E-04; // converts microns to cm.
1772 for(Int_t i=0;i<160;i++) bz[i] = 425.0; // most are 425 microns except below
1773 bz[ 31] = bz[ 32] = 625.0; // first chip boundry
1774 bz[ 63] = bz[ 64] = 625.0; // first chip boundry
1775 bz[ 95] = bz[ 96] = 625.0; // first chip boundry
1776 bz[127] = bz[128] = 625.0; // first chip boundry
1777 Double_t xbins[41]; // each bin in x (Z loc coordinate) includes 4 columns
1778 //xbins[0]=0;
1779 Float_t xmn,xmx,zmn,zmx;
1780 if(!fPlaneEffSPD->GetBlockBoundaries(0,xmn,xmx,zmn,zmx)) AliWarning("Could not book histo properly");
1781 xbins[0]=(Double_t)zmn;
1782 for(Int_t i=0;i<40;i++) {
1783 xbins[i+1]=xbins[i] + (bz[4*i]+bz[4*i+1]+bz[4*i+2]+bz[4*i+3])*kconv;
1784 }
1785 TString histname="ClustersLay1_mod_",aux;
1786 fhClustersInModuleLay1 =new TH2F*[80];
1787 for (Int_t nhist=0;nhist<80;nhist++){
1788 aux=histname;
1789 aux+=nhist;
1790 //
1791 fhClustersInModuleLay1[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1792 fhClustersInModuleLay1[nhist]->SetName(aux.Data());
1793 fhClustersInModuleLay1[nhist]->SetTitle(aux.Data());
1794 fhClustersInModuleLay1[nhist]->SetDirectory(0);
1795 }
1796 histname="ClustersLay2_mod_";
1797 fhClustersInModuleLay2 =new TH2F*[160];
1798 for (Int_t nhist=0;nhist<160;nhist++){
1799 aux=histname;
1800 aux+=nhist;
1801 fhClustersInModuleLay2[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
1802 fhClustersInModuleLay2[nhist]->SetName(aux.Data());
1803 fhClustersInModuleLay2[nhist]->SetTitle(aux.Data());
1804 fhClustersInModuleLay2[nhist]->SetDirectory(0);
1805 }
1806//
275a301c 1807 return;
1808}
1809//____________________________________________________________
1810void AliITSTrackleterSPDEff::DeleteHistos() {
84161aec 1811//
1812// Private method to delete Histograms from memory
1813// it is called. e.g., by the destructor.
58e8dc31 1814//
1815// form AliITSMultReconstructor
1816 if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
1817 if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
1818 if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
1819 if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
1820 if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
1821 if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
1822 if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
1823 if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
1824 if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
1825 if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
1826 if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
1827 if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
1828 if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
1829 if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
84161aec 1830//
275a301c 1831 if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1832 if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1833 if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1834 if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1835 if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1836 if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1837 if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1838 if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1839 if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1840 if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1841 if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1842 if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
17d531c2 1843 if(fhClustersInChip) {delete fhClustersInChip; fhClustersInChip=0;}
1844 if(fhClustersInModuleLay1) {
1845 for (Int_t i=0; i<80; i++ ) delete fhClustersInModuleLay1[i];
1846 delete [] fhClustersInModuleLay1; fhClustersInModuleLay1=0;
1847 }
1848 if(fhClustersInModuleLay2) {
1849 for (Int_t i=0; i<160; i++ ) delete fhClustersInModuleLay2[i];
1850 delete [] fhClustersInModuleLay2; fhClustersInModuleLay2=0;
1851 }
275a301c 1852}
1853//_______________________________________________________________
a3b31967 1854Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
1855 Float_t* vtx, AliStack *stack, TTree *ref) {
1856// This (private) method can be used only for MC events, where both AliStack and the TrackReference
1857// are available.
1858// It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
1859// Input:
1860// - Int_t layer (either 0 or 1): layer which you want to chech if the tracklete can be
1861// reconstructed at
1862// - Int_t iC : cluster index used to build the tracklet prediction
1863// if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
1864// - Float_t* vtx: actual event vertex
1865// - stack: pointer to Stack
1866// - ref: pointer to TTRee of TrackReference
1867Bool_t ret=kFALSE; // returned value
1868Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
1869if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
1870if(!stack) {AliError("null pointer to MC stack"); return ret;}
1871if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
1872if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
1873if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
1874
1875AliTrackReference *tref=0x0;
1876Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
1877Int_t nentries = (Int_t)ref->GetEntries();
1878TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
1879TBranch *br = ref->GetBranch("TrackReferences");
1880br->SetAddress(&tcaRef);
1881for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
1882 br->GetEntry(itrack);
1883 Int_t nref=tcaRef->GetEntriesFast();
1884 if(nref>0) { //it is enough to look at the first one
1885 tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
1886 if(tref->GetTrack()==ipart) {imatch=itrack; break;}
1887 }
1888}
1889if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
1890br->GetEntry(imatch); // redundant, nevertheless ...
1891Int_t nref=tcaRef->GetEntriesFast();
1892for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
1893 tref=(AliTrackReference*)tcaRef->At(iref);
1894 if(tref->R()>10) continue; // not SPD ref
1895 if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
1896 if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
1897
1898// compute the proper quantities for this tref, as was done for fClustersLay1/2
1899 Float_t x = tref->X() - vtx[0];
1900 Float_t y = tref->Y() - vtx[1];
1901 Float_t z = tref->Z() - vtx[2];
1902
1903 Float_t r = TMath::Sqrt(x*x + y*y +z*z);
1904
1905 trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
1906 trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
1907 trefLayExtr[2] = z; // Store z
1908
1909 if(layer==1) { // try to see if it is reconstructable at the outer layer
1910// find the difference in angles
1911 Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
1912 // take into account boundary condition
1913 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1914
1915 // find the difference in z (between linear projection from layer 1
1916 // and the actual point: Dzeta= z1/r1*r2 -z2)
1917 Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1918 Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
1919
1920 // make "elliptical" cut in Phi and Zeta!
7284b2b2 1921 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
1922 dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
a3b31967 1923 if (d<1) {ret=kTRUE; break;}
1924 }
1925 if(layer==0) { // try to see if it is reconstructable at the inner layer
1926
1927 // find the difference in angles
1928 Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
1929 // take into account boundary condition
1930 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
1931
1932 // find the difference in z (between linear projection from layer 2
1933 // and the actual point: Dzeta= z2/r2*r1 -z1)
1934 Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
1935 Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
1936
1937 // make "elliptical" cut in Phi and Zeta!
1938 Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
1939 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
1940 if (d<1) {ret=kTRUE; break;};
1941 }
1942}
1943delete tcaRef;
1944return ret;
1945}
0fce916f 1946//_________________________________________________________________________
1947void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
1948//
1949// this method apply a rotation by 180 degree around the Z (beam) axis to all
1950// the RecPoints in a given layer to be used to build tracklets.
1951// **************** VERY IMPORTANT:: ***************
1952// It must be called just after LoadClusterArrays, since afterwards the datamember
1953// fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
1954// instead of Cartesian
1955//
1956if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
1957AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
1958if(ilayer==0) {
1959 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
1960 fClustersLay1[iC1][0]*=-1;
1961 fClustersLay1[iC1][1]*=-1;
1962 }
1963}
1964if(ilayer==1) {
1965 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
1966 fClustersLay2[iC2][0]*=-1;
1967 fClustersLay2[iC2][1]*=-1;
1968 }
1969}
1970return;
1971}
58e8dc31 1972//____________________________________________________________________________
1973Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *){
1974// This method is used to find the tracklets.
1975// It is called from AliReconstruction
1976// The vertex is supposed to be associated to the Tracker (i.e. to this) already
1977// The cluster is supposed to be associated to the Tracker already
1978// In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
1979//
1980 Int_t rc=1;
1981 AliRunLoader* runLoader = AliRunLoader::Instance();
1982 if (!runLoader) {
1983 Error("Clusters2Tracks", "no run loader found");
1984 return rc;
1985 }
1986 AliStack *pStack=0x0; TTree *tRefTree=0x0;
1987 if(GetMC()) {
1988 runLoader->LoadKinematics("read");
1989 runLoader->LoadTrackRefs("read");
1990 pStack= runLoader->Stack();
1991 tRefTree= runLoader->TreeTR();
1992 }
1993 Reconstruct(pStack,tRefTree);
03ee9629 1994
1995 if (GetLightBkgStudyInParallel()) {
1996 AliStack *dummy1=0x0; TTree *dummy2=0x0;
1997 ReflectClusterAroundZAxisForLayer(1);
1998 Reconstruct(dummy1,dummy2,kTRUE);
1999 }
58e8dc31 2000 return 0;
2001}
2002//____________________________________________________________________________
2003Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
2004//
2005// It is called from AliReconstruction
2006//
2007//
2008//
2009//
2010 Int_t rc=0;
2011 if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
2012 if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
03ee9629 2013 if(GetLightBkgStudyInParallel()) {
2014 TString name="AliITSPlaneEffSPDtrackletBkg.root";
2015 TFile* pefile = TFile::Open(name, "RECREATE");
2016 rc*=fPlaneEffBkg->Write();
2017 pefile->Close();
2018 }
58e8dc31 2019 return rc;
2020}
2021//____________________________________________________________________
2022void
2023AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
2024 // This method
2025 // - gets the clusters from the cluster tree
2026 // - convert them into global coordinates
2027 // - store them in the internal arrays
2028 // - count the number of cluster-fired chips
2029
2030 //AliDebug(1,"Loading clusters and cluster-fired chips ...");
2031
2032 fNClustersLay1 = 0;
2033 fNClustersLay2 = 0;
2034
2035 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
2036 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
2037
2038 itsClusterBranch->SetAddress(&itsClusters);
2039
2040 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
2041 Float_t cluGlo[3]={0.,0.,0.};
2042
2043 // loop over the its subdetectors
2044 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
2045
2046 if (!itsClusterTree->GetEvent(iIts))
2047 continue;
2048
2049 Int_t nClusters = itsClusters->GetEntriesFast();
2050
2051 // number of clusters in each chip of the current module
2052 Int_t layer = 0;
2053
2054 // loop over clusters
2055 while(nClusters--) {
2056 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
2057
2058 layer = cluster->GetLayer();
2059 if (layer>1) continue;
2060
2061 cluster->GetGlobalXYZ(cluGlo);
2062 Float_t x = cluGlo[0];
2063 Float_t y = cluGlo[1];
2064 Float_t z = cluGlo[2];
2065
2066 if (layer==0) {
2067 fClustersLay1[fNClustersLay1][0] = x;
2068 fClustersLay1[fNClustersLay1][1] = y;
2069 fClustersLay1[fNClustersLay1][2] = z;
2070
2071 for (Int_t i=0; i<3; i++)
2072 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
2073 fNClustersLay1++;
17d531c2 2074 if(fHistOn) {
2075 Int_t det=cluster->GetDetectorIndex();
2076 if(det<0 || det>79) {AliError("Cluster with det. index out of boundaries"); return;}
2077 fhClustersInModuleLay1[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2078 }
58e8dc31 2079 }
2080 if (layer==1) {
2081 fClustersLay2[fNClustersLay2][0] = x;
2082 fClustersLay2[fNClustersLay2][1] = y;
2083 fClustersLay2[fNClustersLay2][2] = z;
2084
2085 for (Int_t i=0; i<3; i++)
2086 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
2087 fNClustersLay2++;
17d531c2 2088 if(fHistOn) {
2089 Int_t det=cluster->GetDetectorIndex();
2090 if(det<0 || det>159) {AliError("Cluster with det. index out of boundaries"); return;}
2091 fhClustersInModuleLay2[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
2092 }
58e8dc31 2093 }
2094
2095 }// end of cluster loop
2096
2097 } // end of its "subdetector" loop
2098 if (itsClusters) {
2099 itsClusters->Delete();
2100 delete itsClusters;
2101 itsClusters = 0;
2102 }
2103 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
2104}
03ee9629 2105//_________________________________________________________________________
2106void
2107AliITSTrackleterSPDEff::SetLightBkgStudyInParallel(Bool_t b) {
2108// This method:
2109// - set Bool_t fLightBackgroundStudyInParallel = b
2110// a) if you set this kTRUE, then the estimation of the
2111// SPD efficiency is done as usual for data, but in
2112// parallel a light (i.e. without control histograms, etc.)
2113// evaluation of combinatorial background is performed
2114// with the usual ReflectClusterAroundZAxisForLayer method.
2115// b) if you set this kFALSE, then you would not have a second
2116// container for PlaneEfficiency statistics to be used for background
2117// (fPlaneEffBkg=0). If you want to have a full evaluation of the
2118// background (with all control histograms and additional data
2119// members referring to the background) then you have to call the
2120// method SetReflectClusterAroundZAxisForLayer(kTRUE) esplicitily
2121 fLightBkgStudyInParallel=b;
2122 if(fLightBkgStudyInParallel) {
2123 if(!fPlaneEffBkg) fPlaneEffBkg = new AliITSPlaneEffSPD();
2124 }
2125 else {
2126 delete fPlaneEffBkg;
2127 fPlaneEffBkg=0;
2128 }
2129}