ITS/AliITSTrackleterSPDEff.cxx

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