ITS/AliITSTrackleterSPDEff.cxx

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