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