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