ITS/AliITSTrackleterSPDEff.cxx

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