ITS/AliITSTrackleterSPDEff.cxx

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