ITS/AliITSTrackleterSPDEff.cxx

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