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