ITS/AliITSTrackleterSPDEff.cxx

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