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