ITS/AliITSTrackleterSPDEff.cxx

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