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 /* $Id$ */
  17
  18 //____________________________________________________________________
  19 //
  20 // AliITSTrackleterSPDEff - find SPD chips efficiencies by using tracklets.
  21 //
  22 // This class has been developed from AliITSMultReconstructor (see
  23 // it for more details). It is the class for the Trackleter used to estimate
  24 // SPD plane efficiency.
  25 // The trackleter prediction is built using the vertex and 1 cluster.
  26 //
  27 //
  28 //  Author :  Giuseppe Eugenio Bruno, based on the skeleton of Reconstruct method  provided by Tiziano Virgili
  29 //  email:    giuseppe.bruno@ba.infn.it
  30 //
  31 //____________________________________________________________________
  32
  33 #include <TFile.h>
  34 #include <TParticle.h>
  35 #include <TSystem.h>
  36 #include <Riostream.h>
  37
  38 #include "AliITSMultReconstructor.h"
  39 #include "AliITSTrackleterSPDEff.h"
  40 #include "AliITSgeomTGeo.h"
  41 #include "AliLog.h"
  42 #include "AliITSPlaneEffSPD.h"
  43 #include "AliStack.h"
  44
  45 //____________________________________________________________________
  46 ClassImp(AliITSTrackleterSPDEff)
  47
  48
  49 //____________________________________________________________________
  50 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
  51 AliITSMultReconstructor(),
  52 fAssociationFlag1(0),
  53 fChipPredOnLay2(0),
  54 fChipPredOnLay1(0),
  55 fNTracklets1(0),
  56 fPhiWindowL1(0),
  57 fZetaWindowL1(0),
  58 fOnlyOneTrackletPerC1(0),
  59 fPlaneEffSPD(0),
  60 fMC(0),
  61 fUseOnlyPrimaryForPred(0),
  62 fUseOnlySecondaryForPred(0),
  63 fUseOnlySameParticle(0),
  64 fUseOnlyDifferentParticle(0),
  65 fUseOnlyStableParticle(0),
  66 fPredictionPrimary(0),
  67 fPredictionSecondary(0),
  68 fClusterPrimary(0),
  69 fClusterSecondary(0),
  70 fhClustersDPhiInterpAcc(0),
  71 fhClustersDThetaInterpAcc(0),
  72 fhClustersDZetaInterpAcc(0),
  73 fhClustersDPhiInterpAll(0),
  74 fhClustersDThetaInterpAll(0),
  75 fhClustersDZetaInterpAll(0),
  76 fhDPhiVsDThetaInterpAll(0),
  77 fhDPhiVsDThetaInterpAcc(0),
  78 fhDPhiVsDZetaInterpAll(0),
  79 fhDPhiVsDZetaInterpAcc(0),
  80 fhetaClustersLay2(0),
  81 fhphiClustersLay2(0)
  82 {
  83
  84   // Method to check the SPD chips efficiencies by using tracklets
  85
  86
  87   SetPhiWindowL1();
  88   SetZetaWindowL1();
  89   SetOnlyOneTrackletPerC1();
  90
  91   fAssociationFlag1   = new Bool_t[300000];
  92   fChipPredOnLay2     = new UInt_t[300000];
  93   fChipPredOnLay1     = new UInt_t[300000];
  94
  95   for(Int_t i=0; i<300000; i++) {
  96     fAssociationFlag1[i]   = kFALSE;
  97   }
  98
  99   if (GetHistOn()) BookHistos();
 100
 101   fPlaneEffSPD = new AliITSPlaneEffSPD();
 102 }
 103 //______________________________________________________________________
 104 AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) : AliITSMultReconstructor(mr),
 105 fAssociationFlag1(mr.fAssociationFlag1),
 106 fChipPredOnLay2(mr.fChipPredOnLay2),
 107 fChipPredOnLay1(mr.fChipPredOnLay1),
 108 fNTracklets1(mr.fNTracklets1),
 109 fPhiWindowL1(mr.fPhiWindowL1),
 110 fZetaWindowL1(mr.fZetaWindowL1),
 111 fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
 112 fPlaneEffSPD(mr.fPlaneEffSPD),
 113 fMC(mr.fMC),
 114 fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
 115 fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
 116 fUseOnlySameParticle(mr.fUseOnlySameParticle),
 117 fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle),
 118 fUseOnlyStableParticle(mr.fUseOnlyStableParticle),
 119 fPredictionPrimary(mr.fPredictionPrimary),
 120 fPredictionSecondary(mr.fPredictionSecondary),
 121 fClusterPrimary(mr.fClusterPrimary),
 122 fClusterSecondary(mr.fClusterSecondary),
 123 fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
 124 fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
 125 fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
 126 fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll),
 127 fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll),
 128 fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll),
 129 fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll),
 130 fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc),
 131 fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
 132 fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
 133 fhetaClustersLay2(mr.fhetaClustersLay2),
 134 fhphiClustersLay2(mr.fhphiClustersLay2)
 135 {
 136   // Copy constructor
 137 }
 138
 139 //______________________________________________________________________
 140 AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){
 141   // Assignment operator
 142   this->~AliITSTrackleterSPDEff();
 143   new(this) AliITSTrackleterSPDEff(mr);
 144   return *this;
 145 }
 146 //______________________________________________________________________
 147 AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
 148   // Destructor
 149
 150   // delete histograms
 151   DeleteHistos();
 152
 153   delete [] fAssociationFlag1;
 154
 155   delete [] fChipPredOnLay2;
 156   delete [] fChipPredOnLay1;
 157
 158   delete [] fPredictionPrimary;
 159   delete [] fPredictionSecondary;
 160   delete [] fClusterPrimary;
 161   delete [] fClusterSecondary;
 162
 163   // delete PlaneEff
 164   delete fPlaneEffSPD;
 165 }
 166 //____________________________________________________________________
 167 void
 168 AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/,AliStack *pStack) {
 169   //
 170   // - calls LoadClusterArray that finds the position of the clusters
 171   //   (in global coord)
 172   // - convert the cluster coordinates to theta, phi (seen from the
 173   //   interaction vertex). Find the extrapolation/interpolation point.
 174   // - Find the chip corresponding to that
 175   // - Check if there is a cluster near that point
 176   //
 177
 178   // reset counters
 179   fNClustersLay1 = 0;
 180   fNClustersLay2 = 0;
 181   fNTracklets = 0;
 182   fNSingleCluster = 0;
 183   // loading the clusters
 184   LoadClusterArrays(clusterTree);
 185   if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
 186   Bool_t found;
 187   Int_t nfTraPred1=0;  Int_t ntTraPred1=0;
 188   Int_t nfTraPred2=0;  Int_t ntTraPred2=0;
 189   Int_t nfClu1=0;      Int_t ntClu1=0;
 190   Int_t nfClu2=0;      Int_t ntClu2=0;
 191
 192
 193   // find the tracklets
 194   AliDebug(1,"Looking for tracklets... ");
 195   AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
 196
 197   //###########################################################
 198   // Loop on layer 1 : finding theta, phi and z
 199   UInt_t key;
 200   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 201     Float_t x = fClustersLay1[iC1][0] - vtx[0];
 202     Float_t y = fClustersLay1[iC1][1] - vtx[1];
 203     Float_t z = fClustersLay1[iC1][2] - vtx[2];
 204
 205     Float_t r    = TMath::Sqrt(x*x + y*y +z*z);
 206
 207     fClustersLay1[iC1][0] = TMath::ACos(z/r);                   // Store Theta
 208     fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 209     fClustersLay1[iC1][2] = z;                  // Store z
 210
 211     // find the Radius and the chip corresponding to the extrapolation point
 212
 213     found=FindChip(key, 1, vtx, fClustersLay1[iC1][0],fClustersLay1[iC1][1]);
 214     if (!found) {
 215       AliDebug(1,Form("Reconstruct: cannot find chip prediction on outer layer for cluster %d on the inner layer",iC1));
 216       key=999999;               // also some other actions should be taken if not Found
 217     }
 218     nfTraPred2+=(Int_t)found; // this for debugging purpose
 219     ntTraPred2++;             // to check efficiency of the method FindChip
 220     fChipPredOnLay2[iC1] = key;
 221     fAssociationFlag1[iC1] = kFALSE;
 222
 223     if (fHistOn) {
 224       Float_t eta=fClustersLay1[iC1][0];
 225       eta= TMath::Tan(eta/2.);
 226       eta=-TMath::Log(eta);
 227       fhetaClustersLay1->Fill(eta);
 228       fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
 229     }
 230   }
 231
 232   // Loop on layer 2 : finding theta, phi and r
 233   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 234     Float_t x = fClustersLay2[iC2][0] - vtx[0];
 235     Float_t y = fClustersLay2[iC2][1] - vtx[1];
 236     Float_t z = fClustersLay2[iC2][2] - vtx[2];
 237
 238     Float_t r    = TMath::Sqrt(x*x + y*y +z*z);
 239
 240     fClustersLay2[iC2][0] = TMath::ACos(z/r);                   // Store Theta
 241     fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi (done properly in the range [0,2pi])
 242     fClustersLay2[iC2][2] = z;                  // Store z
 243
 244     // find the Radius and the chip corresponding to the extrapolation point
 245
 246     found=FindChip(key, 0, vtx, fClustersLay2[iC2][0],fClustersLay2[iC2][1]);
 247     if (!found) {
 248       AliWarning(Form("Reconstruct: cannot find chip prediction on inner layer for cluster %d on the outer layer",iC2));
 249       key=999999;
 250     }
 251     nfTraPred1+=(Int_t)found; // this for debugging purpose
 252     ntTraPred1++;             // to check efficiency of the method FindChip
 253     fChipPredOnLay1[iC2] = key;
 254     fAssociationFlag[iC2] = kFALSE;
 255
 256     if (fHistOn) {
 257       Float_t eta=fClustersLay2[iC2][0];
 258       eta= TMath::Tan(eta/2.);
 259       eta=-TMath::Log(eta);
 260       fhetaClustersLay2->Fill(eta);
 261       fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
 262     }
 263   }
 264
 265   //###########################################################
 266
 267  // First part : Extrapolation to Layer 2
 268
 269   // Loop on layer 1
 270   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 271
 272     // reset of variables for multiple candidates
 273     Int_t  iC2WithBestDist = 0;     // reset
 274     Float_t distmin        = 100.;  // just to put a huge number!
 275     Float_t dPhimin        = 0.;  // Used for histograms only!
 276     Float_t dThetamin      = 0.;  // Used for histograms only!
 277     Float_t dZetamin       = 0.;  // Used for histograms only!
 278
 279     // in any case, if MC has been required, store statistics of primaries and secondaries
 280     if (fMC) {
 281        Int_t lab1=(Int_t)fClustersLay1[iC1][3];
 282        Int_t lab2=(Int_t)fClustersLay1[iC1][4];
 283        Int_t lab3=(Int_t)fClustersLay1[iC1][5];
 284        // do it always as a function of the chip number used to built the prediction
 285        found=FindChip(key,0,vtx,fClustersLay1[iC1][0],fClustersLay1[iC1][1],fClustersLay1[iC1][2]);
 286        if (!found) {AliWarning(
 287          Form("Reconstruct MC: cannot find chip on inner layer for cluster %d",iC1)); }
 288        else {
 289          if((lab1 != -2  &&  PrimaryTrackChecker(lab1,pStack) ) ||
 290             (lab2 != -2  &&  PrimaryTrackChecker(lab2,pStack) ) ||
 291             (lab3 != -2  &&  PrimaryTrackChecker(lab3,pStack)))
 292          { // this cluster is from a primary particle
 293            fClusterPrimary[key]++;
 294            if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
 295          } else { // this cluster is from a secondary particle
 296             fClusterSecondary[key]++;
 297             if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
 298          }
 299        }
 300        // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
 301        // (in case the prediction is reliable)
 302        if( fChipPredOnLay2[iC1]<1200) {
 303          if((lab1 != -2  &&  PrimaryTrackChecker(lab1,pStack) ) ||
 304             (lab2 != -2  &&  PrimaryTrackChecker(lab2,pStack) ) ||
 305             (lab3 != -2  &&  PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
 306          else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
 307        }
 308     }
 309
 310     // Loop on layer 2
 311     for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 312
 313       // The following excludes double associations
 314       if (!fAssociationFlag[iC2]) {
 315
 316         // find the difference in angles
 317         Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
 318         Float_t dPhi   = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
 319         // take into account boundary condition
 320         if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
 321
 322         // find the difference in z (between linear projection from layer 1
 323         // and the actual point: Dzeta= z1/r1*r2 -z2)
 324         Float_t r2    = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
 325         Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
 326
 327         if (fHistOn) {
 328           fhClustersDPhiAll->Fill(dPhi);
 329           fhClustersDThetaAll->Fill(dTheta);
 330           fhClustersDZetaAll->Fill(dZeta);
 331           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
 332           fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
 333         }
 334
 335         // make "elliptical" cut in Phi and Zeta!
 336         Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
 337                                 dZeta*dZeta/fZetaWindow/fZetaWindow);
 338
 339         if (d>1) continue;
 340
 341         //look for the minimum distance: the minimum is in iC2WithBestDist
 342         if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
 343           distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
 344           dPhimin = dPhi;
 345           dThetamin = dTheta;
 346           dZetamin = dZeta;
 347           iC2WithBestDist = iC2;
 348         }
 349       }
 350     } // end of loop over clusters in layer 2
 351
 352     if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
 353
 354       if (fHistOn) {
 355         fhClustersDPhiAcc->Fill(dPhimin);
 356         fhClustersDThetaAcc->Fill(dThetamin);
 357         fhClustersDZetaAcc->Fill(dZetamin);
 358         fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
 359         fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
 360       }
 361
 362       if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE;
 363        // flag the association
 364
 365       // store the tracklet
 366
 367       // use the theta from the clusters in the first layer
 368       fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
 369       // use the phi from the clusters in the first layer
 370       fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
 371       // Store the difference between phi1 and phi2
 372       fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
 373
 374       // find labels
 375       Int_t label1 = 0;
 376       Int_t label2 = 0;
 377       while (label2 < 3)
 378       {
 379         if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
 380           break;
 381         label1++;
 382         if (label1 == 3)
 383         {
 384           label1 = 0;
 385           label2++;
 386         }
 387       }
 388
 389       if (label2 < 3)
 390       {
 391         fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
 392       }
 393       else
 394       {
 395         fTracklets[fNTracklets][3] = -2;
 396       }
 397
 398       if (fHistOn) {
 399         Float_t eta=fTracklets[fNTracklets][0];
 400         eta= TMath::Tan(eta/2.);
 401         eta=-TMath::Log(eta);
 402         fhetaTracklets->Fill(eta);
 403         fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
 404       }
 405
 406 // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation)
 407       found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]);
 408       if(!found){
 409         AliWarning(
 410          Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist));
 411         key=999999;
 412       }
 413       nfClu2+=(Int_t)found; // this for debugging purpose
 414       ntClu2++;             // to check efficiency of the method FindChip
 415       if(key<1200) { // the Chip has been found
 416         if(fMC) { // this part only for MC
 417           // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
 418           // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
 419           // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
 420           if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
 421           if (fUseOnlySameParticle && label2 == 3) continue;      // different label (reject it)
 422         }
 423
 424         if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
 425           // OK, success
 426                 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
 427         }
 428         else {
 429                 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
 430                                                          // (might be in the tracking tollerance)
 431         }
 432       }
 433
 434       fNTracklets++;
 435
 436     } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
 437     else if (fChipPredOnLay2[iC1]<1200) fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
 438
 439   } // end of loop over clusters in layer 1
 440
 441     fNTracklets1=fNTracklets;
 442
 443 //###################################################################
 444
 445   // Second part : Interpolation to Layer 1
 446
 447   // Loop on layer 2
 448   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 449
 450     // reset of variables for multiple candidates
 451     Int_t  iC1WithBestDist = 0;     // reset
 452     Float_t distmin        = 100.;  // just to put a huge number!
 453     Float_t dPhimin        = 0.;  // Used for histograms only!
 454     Float_t dThetamin      = 0.;  // Used for histograms only!
 455     Float_t dZetamin       = 0.;  // Used for histograms only!
 456
 457     // in any case, if MC has been required, store statistics of primaries and secondaries
 458     if (fMC) {
 459        Int_t lab1=(Int_t)fClustersLay2[iC2][3];
 460        Int_t lab2=(Int_t)fClustersLay2[iC2][4];
 461        Int_t lab3=(Int_t)fClustersLay2[iC2][5];
 462        // do it always as a function of the chip number used to built the prediction
 463        found=FindChip(key,1,vtx,fClustersLay2[iC2][0],fClustersLay2[iC2][1],fClustersLay2[iC2][2]);
 464        if (!found) {AliWarning(
 465          Form("Reconstruct MC: cannot find chip on outer layer for cluster %d",iC2)); }
 466        else {
 467          if((lab1 != -2  &&  PrimaryTrackChecker(lab1,pStack) ) ||
 468             (lab2 != -2  &&  PrimaryTrackChecker(lab2,pStack) ) ||
 469             (lab3 != -2  &&  PrimaryTrackChecker(lab3,pStack)))
 470          {  // this cluster is from a primary particle
 471             fClusterPrimary[key]++;
 472             if(fUseOnlySecondaryForPred) continue; //  skip this tracklet built with a primary track
 473          } else { // this cluster is from a secondary particle
 474            fClusterSecondary[key]++;
 475            if(fUseOnlyPrimaryForPred) continue; //  skip this tracklet built with a secondary track
 476          }
 477        }
 478        // do it as a function of the chip number where you exspect the cluster (i.e. tracklet prediction)
 479        // (in case the prediction is reliable)
 480        if( fChipPredOnLay1[iC2]<1200) {
 481          if((lab1 != -2  &&  PrimaryTrackChecker(lab1,pStack) ) ||
 482             (lab2 != -2  &&  PrimaryTrackChecker(lab2,pStack) ) ||
 483             (lab3 != -2  &&  PrimaryTrackChecker(lab3,pStack)))   fPredictionPrimary[fChipPredOnLay1[iC2]]++;
 484          else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
 485        }
 486     }
 487
 488     // Loop on layer 1
 489     for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 490
 491       // The following excludes double associations
 492       if (!fAssociationFlag1[iC1]) {
 493
 494         // find the difference in angles
 495         Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
 496         Float_t dPhi   = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
 497         // take into account boundary condition
 498         if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
 499
 500
 501         // find the difference in z (between linear projection from layer 2
 502         // and the actual point: Dzeta= z2/r2*r1 -z1)
 503         Float_t r1    = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]);
 504         Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
 505
 506
 507         if (fHistOn) {
 508           fhClustersDPhiInterpAll->Fill(dPhi);
 509           fhClustersDThetaInterpAll->Fill(dTheta);
 510           fhClustersDZetaInterpAll->Fill(dZeta);
 511           fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi);
 512           fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi);
 513         }
 514         // make "elliptical" cut in Phi and Zeta!
 515         Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
 516                                 dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
 517
 518         if (d>1) continue;
 519
 520         //look for the minimum distance: the minimum is in iC1WithBestDist
 521         if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) {
 522           distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi));
 523           dPhimin = dPhi;
 524           dThetamin = dTheta;
 525           dZetamin = dZeta;
 526           iC1WithBestDist = iC1;
 527         }
 528       }
 529     } // end of loop over clusters in layer 1
 530
 531     if (distmin<100) { // This means that a cluster in layer 1 was found that mathes with iC2
 532
 533       if (fHistOn) {
 534         fhClustersDPhiInterpAcc->Fill(dPhimin);
 535         fhClustersDThetaInterpAcc->Fill(dThetamin);
 536         fhClustersDZetaInterpAcc->Fill(dZetamin);
 537         fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin);
 538         fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin);
 539       }
 540
 541       if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association
 542        // flag the association
 543
 544       // store the tracklet
 545
 546       // use the theta from the clusters in the first layer
 547       fTracklets[fNTracklets][0] = fClustersLay2[iC2][0];
 548       // use the phi from the clusters in the first layer
 549       fTracklets[fNTracklets][1] = fClustersLay2[iC2][1];
 550       // Store the difference between phi1 and phi2
 551       fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1];
 552
 553       // find labels
 554       Int_t label1 = 0;
 555       Int_t label2 = 0;
 556       while (label2 < 3)
 557       {
 558         if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2])
 559           break;
 560         label1++;
 561         if (label1 == 3)
 562         {
 563           label1 = 0;
 564           label2++;
 565         }
 566       }
 567
 568       if (label2 < 3)
 569       {
 570         fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1];
 571       }
 572       else
 573       {
 574         fTracklets[fNTracklets][3] = -2;
 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,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]);
 579       if(!found){
 580         AliWarning(
 581          Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist));
 582         key=999999;
 583       }
 584       nfClu1+=(Int_t)found; // this for debugging purpose
 585       ntClu1++;             // to check efficiency of the method FindChip
 586       if(key<1200) {
 587         if(fMC) { // this part only for MC
 588           // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
 589           // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
 590           // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
 591           if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
 592           if (fUseOnlySameParticle && label2 == 3) continue;      // different label (reject it)
 593         }
 594
 595         if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
 596           // OK, success
 597                 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
 598         } else {
 599                 fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
 600                                                          // (might be in the tracking tollerance)
 601         }
 602       }
 603
 604     fNTracklets++;
 605
 606     } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
 607     else if (fChipPredOnLay1[iC2]<1200) fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
 608
 609   } // end of loop over clusters in layer 2
 610
 611   AliDebug(1,Form("%d tracklets found", fNTracklets));
 612   AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1));
 613   AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2));
 614   AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1));
 615   AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
 616 }
 617 //____________________________________________________________________
 618 Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer,  Float_t* vtx,
 619                                   Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
 620 //
 621 // Input: a) layer number in the range [0,1]
 622 //        b) vtx[3]: actual vertex
 623 //        c) zVtx     \ z of the cluster (-999 for tracklet) computed with respect to vtx
 624 //        d) thetaVtx  > theta and phi of the cluster/tracklet computed with respect to vtx
 625 //        e) phiVtx   /
 626 // Output: Unique key to locate a chip
 627 // return: kTRUE if succesfull
 628
 629     if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;}
 630     Double_t r=GetRLayer(layer);
 631     //AliInfo(Form("Radius on layer %d  is %f cm",layer,r));
 632
 633   // set phiVtx in the range [0,2pi]
 634   if(!SetAngleRange02Pi(phiVtx)) return kFALSE ;
 635
 636   Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
 637                         // of the intersection of the tracklet with the pixel layer.
 638   if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
 639   else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
 640   AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
 641   Double_t vtxy[2]={vtx[0],vtx[1]};
 642   if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
 643     // this method gives you two interceptions
 644     if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE;
 645     // set phiAbs in the range [0,2pi]
 646     if(!SetAngleRange02Pi(phiAbs)) return kFALSE;
 647     // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close;
 648     // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by
 649     // taking the closest one to phiVtx
 650     AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs));
 651   } else phiAbs=phiVtx;
 652   Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number
 653
 654   // now you need to locate the chip within the idet detector,
 655   // starting from the local coordinates in such a detector
 656
 657   Float_t locx; // local Cartesian coordinate (to be determined) corresponding to
 658   Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) .
 659   if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE;
 660
 661   key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz);
 662   return kTRUE;
 663 }
 664 //______________________________________________________________________________
 665 Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
 666     if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
 667     Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
 668
 669     Double_t xyz[3], &x=xyz[0], &y=xyz[1];
 670     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
 671     Double_t r=TMath::Sqrt(x*x + y*y);
 672
 673     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
 674     r += TMath::Sqrt(x*x + y*y);
 675     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
 676     r += TMath::Sqrt(x*x + y*y);
 677     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
 678     r += TMath::Sqrt(x*x + y*y);
 679     r*=0.25;
 680     return r;
 681 }
 682 //______________________________________________________________________________
 683 Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z,
 684                            Float_t &xloc, Float_t &zloc) {
 685   // this method transform Global Cilindrical coordinates into local (i.e. module)
 686   // cartesian coordinates
 687   //
 688   //Compute Cartesian Global Coordinate
 689   Double_t xyzGlob[3],xyzLoc[3];
 690   xyzGlob[2]=z;
 691   xyzGlob[0]=r*TMath::Cos(phi);
 692   xyzGlob[1]=r*TMath::Sin(phi);
 693
 694   xloc=0.;
 695   zloc=0.;
 696
 697   if(idet<0)  return kFALSE;
 698
 699   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
 700   Int_t lad = Int_t(idet/ndet) + 1;
 701   Int_t det = idet - (lad-1)*ndet + 1;
 702
 703   AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
 704
 705   xloc = (Float_t)xyzLoc[0];
 706   zloc = (Float_t)xyzLoc[2];
 707
 708 return kTRUE;
 709 }
 710 //______________________________________________________________________________
 711 Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
 712   //--------------------------------------------------------------------
 713   //This function finds the detector crossed by the track
 714   //--------------------------------------------------------------------
 715     if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
 716     Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
 717     Int_t nladders=AliITSgeomTGeo::GetNLadders(i);
 718     Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i);
 719
 720     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2];
 721     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
 722     Double_t phiOffset=TMath::ATan2(y,x);
 723     Double_t zOffset=z2;
 724
 725   Double_t dphi;
 726   if (zOffset<0)            // old geometry
 727     dphi = -(phi-phiOffset);
 728   else                       // new geometry
 729     dphi = phi-phiOffset;
 730
 731   if      (dphi <  0) dphi += 2*TMath::Pi();
 732   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
 733   Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5);
 734   if (np>=nladders) np-=nladders;
 735   if (np<0)          np+=nladders;
 736
 737   Double_t dz=zOffset-z;
 738   Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5;
 739   Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
 740   if (nz>=ndetectors) {AliDebug(1,Form("too  long: nz =%d",nz)); return -1;}
 741   if (nz<0)           {AliDebug(1,Form("too short: nz =%d",nz)); return -1;}
 742
 743   return np*ndetectors + nz;
 744 }
 745 //____________________________________________________________
 746 Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) {
 747 // this method find the intersection in xy between a tracklet (straight line) and
 748 // a circonference (r=R), using polar coordinates.
 749 /*
 750 Input: - vtx[2]: actual vertex w.r.t. ALICE reference system
 751        - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx
 752        - R: radius of the circle
 753 Output: - phi : phi angle of the unique interception in the ALICE Global ref. system
 754
 755 Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx:
 756 r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
 757 In the same system, the equation of a semi-line is: phi=phiVtx;
 758 Hence you get one interception only: P=(r,phiVtx)
 759 Finally you want P in the ABSOLUTE ALICE system.
 760 */
 761 Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
 762 Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]);          // in the system with vtx[2] as origin
 763 Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO);
 764 Double_t cC=rO*rO-R*R;
 765 Double_t dDelta=bB*bB-4*cC;
 766 if(dDelta<0) return kFALSE;
 767 Double_t r1,r2;
 768 r1=(-bB-TMath::Sqrt(dDelta))/2;
 769 r2=(-bB+TMath::Sqrt(dDelta))/2;
 770 if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;}
 771 Double_t r=TMath::Max(r1,r2); // take the positive
 772 Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx
 773 Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin
 774 pvtx[0]=r*TMath::Cos(phiVtx);
 775 pvtx[1]=r*TMath::Sin(phiVtx);
 776 pP[0]=vtx[0]+pvtx[0];
 777 pP[1]=vtx[1]+pvtx[1];
 778 phi=TMath::ATan2(pP[1],pP[0]);
 779 return kTRUE;
 780 }
 781 //___________________________________________________________
 782 Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
 783 while(angle >=2*TMath::Pi() || angle<0) {
 784   if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
 785   if(angle < 0) angle+=2*TMath::Pi();
 786 }
 787 return kTRUE;
 788 }
 789 //___________________________________________________________
 790 Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
 791 if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
 792 if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
 793 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
 794 // return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz)
 795  if(!stack->IsPhysicalPrimary(ipart)) return kFALSE;
 796  // like below: as in the correction for Multiplicity (i.e. by hand in macro)
 797  TParticle* part = stack->Particle(ipart);
 798  TParticle* part0 = stack->Particle(0); // first primary
 799  TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary
 800  if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f",
 801             part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() ));
 802
 803  if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; }
 804  TParticlePDG* pdgPart = part->GetPDG();
 805  if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;}
 806
 807   Double_t distx = part->Vx() - part0->Vx();
 808   Double_t disty = part->Vy() - part0->Vy();
 809   Double_t distz = part->Vz() - part0->Vz();
 810   Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz);
 811
 812   if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ",
 813                                  part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
 814                       return kFALSE; }// primary if within 500 microns from true Vertex
 815
 816  if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)<2) return kFALSE;
 817  return kTRUE;
 818 }
 819 //_____________________________________________________________________________________________
 820 Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
 821 if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
 822 if(!stack) {AliError("null pointer to MC stack"); return 0;}
 823 if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
 824
 825 TParticle* part = stack->Particle(ipart);
 826 //TParticle* part0 = stack->Particle(0); // first primary
 827
 828   Int_t nret=0;
 829   TParticle* dau = 0;
 830   Int_t nDau = 0;
 831   Int_t firstDau = part->GetFirstDaughter();
 832   if (firstDau > 0) {
 833     Int_t lastDau = part->GetLastDaughter();
 834     nDau = lastDau - firstDau + 1;
 835     //printf("number of daugthers %d \n",nDau);
 836     if (nDau > 0) {
 837       //for(Int_t j=firstDau; j<=lastDau; j++)
 838       for(Int_t j=firstDau; j<=firstDau; j++)
 839                                               { // only first one
 840         dau = stack->Particle(j);
 841         Double_t distx = dau->Vx()-part->Vx();
 842         Double_t disty = dau->Vy()-part->Vy();
 843         Double_t distz = dau->Vz()-part->Vz();
 844         Double_t distR = TMath::Sqrt(distx*distx+disty*disty+distz*distz);
 845         if (distR > GetRLayer(0)+0.5)  nret=1;  // decay after first pixel layer
 846         if (distR > GetRLayer(1)+0.5)  nret=2;  // decay after second pixel layer
 847       }
 848     }
 849   } else nret = 3; // stable particle
 850 return nret;
 851 }
 852 //_________________________________________________________________
 853 void AliITSTrackleterSPDEff::InitPredictionMC() {
 854 if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
 855 fPredictionPrimary   = new Int_t[1200];
 856 fPredictionSecondary = new Int_t[1200];
 857 fClusterPrimary      = new Int_t[1200];
 858 fClusterSecondary    = new Int_t[1200];
 859 for(Int_t i=0; i<1200; i++) {
 860  fPredictionPrimary[i]=0;
 861  fPredictionSecondary[i]=0;
 862  fPredictionSecondary[i]=0;
 863  fClusterSecondary[i]=0;
 864 }
 865 return;
 866 }
 867 //______________________________________________________________________
 868 Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
 869 if (!fMC) {CallWarningMC(); return 0;}
 870 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
 871 return fPredictionPrimary[(Int_t)key];
 872 }
 873 //______________________________________________________________________
 874 Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
 875 if (!fMC) {CallWarningMC(); return 0;}
 876 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
 877 return fPredictionSecondary[(Int_t)key];
 878 }
 879 //______________________________________________________________________
 880 Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
 881 if (!fMC) {CallWarningMC(); return 0;}
 882 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
 883 return fClusterPrimary[(Int_t)key];
 884 }
 885 //______________________________________________________________________
 886 Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
 887 if (!fMC) {CallWarningMC(); return 0;}
 888 if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
 889 return fClusterSecondary[(Int_t)key];
 890 }
 891 //______________________________________________________________________
 892 void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
 893     // Print out some class data values in Ascii Form to output stream
 894     // Inputs:
 895     //   ostream *os   Output stream where Ascii data is to be writen
 896     // Outputs:
 897     //   none.
 898     // Return:
 899     //   none.
 900     *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow ;
 901     *os << " " << fMC;
 902     if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
 903     *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
 904         << " " << fUseOnlySameParticle   << " " << fUseOnlyDifferentParticle
 905         << " " << fUseOnlyStableParticle ;
 906     for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i)  ;
 907     for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
 908     for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
 909     for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
 910     return;
 911 }
 912 //______________________________________________________________________
 913 void AliITSTrackleterSPDEff::ReadAscii(istream *is){
 914     // Read in some class data values in Ascii Form to output stream
 915     // Inputs:
 916     //   istream *is   Input stream where Ascii data is to be read in from
 917     // Outputs:
 918     //   none.
 919     // Return:
 920     //   none.
 921
 922     *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow;
 923     *is >> fMC;
 924     if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
 925     *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
 926         >> fUseOnlySameParticle   >> fUseOnlyDifferentParticle
 927         >> fUseOnlyStableParticle;
 928     for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
 929     for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
 930     for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
 931     for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
 932     return;
 933 }
 934 //______________________________________________________________________
 935 ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){
 936     // Standard output streaming function
 937     // Inputs:
 938     //   ostream            &os  output steam
 939     //   AliITSTrackleterSPDEff &s class to be streamed.
 940     // Output:
 941     //   none.
 942     // Return:
 943     //   ostream &os  The stream pointer
 944
 945     s.PrintAscii(&os);
 946     return os;
 947 }
 948 //______________________________________________________________________
 949 istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){
 950     // Standard inputput streaming function
 951     // Inputs:
 952     //   istream            &is  input steam
 953     //   AliITSTrackleterSPDEff &s class to be streamed.
 954     // Output:
 955     //   none.
 956     // Return:
 957     //   ostream &os  The stream pointer
 958
 959     //printf("prova %d \n", (Int_t)s.GetMC());
 960     s.ReadAscii(&is);
 961     return is;
 962 }
 963 //______________________________________________________________________
 964 void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
 965  if(!fMC) {CallWarningMC(); return;}
 966  ofstream out(filename.Data(),ios::out | ios::binary);
 967  out << *this;
 968  out.close();
 969 return;
 970 }
 971 //____________________________________________________________________
 972 void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
 973  if(!fMC) {CallWarningMC(); return;}
 974  if( gSystem->AccessPathName( filename.Data() ) ) {
 975       AliError( Form( "file (%s) not found", filename.Data() ) );
 976       return;
 977    }
 978
 979  ifstream in(filename.Data(),ios::in | ios::binary);
 980  in >> *this;
 981  in.close();
 982  return;
 983 }
 984 //____________________________________________________________________
 985 Bool_t AliITSTrackleterSPDEff::SaveHists() {
 986   // This method save the histograms on the output file
 987   // (only if fHistOn is TRUE).
 988
 989   if (!GetHistOn()) return kFALSE;
 990
 991   AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
 992
 993   fhClustersDPhiInterpAll->Write();
 994   fhClustersDThetaInterpAll->Write();
 995   fhClustersDZetaInterpAll->Write();
 996   fhDPhiVsDThetaInterpAll->Write();
 997   fhDPhiVsDZetaInterpAll->Write();
 998
 999   fhClustersDPhiInterpAcc->Write();
1000   fhClustersDThetaInterpAcc->Write();
1001   fhClustersDZetaInterpAcc->Write();
1002   fhDPhiVsDThetaInterpAcc->Write();
1003   fhDPhiVsDZetaInterpAcc->Write();
1004
1005   fhetaClustersLay2->Write();
1006   fhphiClustersLay2->Write();
1007   return kTRUE;
1008 }
1009 //__________________________________________________________
1010 Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
1011   //
1012   // Saves the histograms into a tree and saves the trees into a file
1013   //
1014   if (!GetHistOn()) return kFALSE;
1015   if (filename.Data()=="") {
1016      AliWarning("WriteHistosToFile: null output filename!");
1017      return kFALSE;
1018   }
1019   TFile *hFile=new TFile(filename.Data(),option,
1020                          "The File containing the histos for SPD efficiency studies with tracklets");
1021   if(!SaveHists()) return kFALSE;
1022   hFile->Write();
1023   hFile->Close();
1024   return kTRUE;
1025 }
1026 //____________________________________________________________
1027 void AliITSTrackleterSPDEff::BookHistos() {
1028   if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
1029   fhClustersDPhiInterpAcc   = new TH1F("dphiaccInterp",  "dphi Interpolation phase",  100,0.,0.1);
1030   fhClustersDPhiInterpAcc->SetDirectory(0);
1031   fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
1032   fhClustersDThetaInterpAcc->SetDirectory(0);
1033   fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.);
1034   fhClustersDZetaInterpAcc->SetDirectory(0);
1035
1036   fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1);
1037   fhDPhiVsDZetaInterpAcc->SetDirectory(0);
1038   fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1);
1039   fhDPhiVsDThetaInterpAcc->SetDirectory(0);
1040
1041   fhClustersDPhiInterpAll   = new TH1F("dphiallInterp",  "dphi Interpolation phase",  100,0.0,0.5);
1042   fhClustersDPhiInterpAll->SetDirectory(0);
1043   fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5);
1044   fhClustersDThetaInterpAll->SetDirectory(0);
1045   fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.);
1046   fhClustersDZetaInterpAll->SetDirectory(0);
1047
1048   fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5);
1049   fhDPhiVsDZetaInterpAll->SetDirectory(0);
1050   fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5);
1051   fhDPhiVsDThetaInterpAll->SetDirectory(0);
1052
1053   fhetaClustersLay2  = new TH1F("etaClustersLay2",  "etaCl2",  100,-2.,2.);
1054   fhetaClustersLay2->SetDirectory(0);
1055   fhphiClustersLay2  = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
1056   fhphiClustersLay2->SetDirectory(0);
1057   return;
1058 }
1059 //____________________________________________________________
1060 void AliITSTrackleterSPDEff::DeleteHistos() {
1061     if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
1062     if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
1063     if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
1064     if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;}
1065     if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;}
1066     if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;}
1067     if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;}
1068     if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;}
1069     if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;}
1070     if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
1071     if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
1072     if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
1073 }
1074 //_______________________________________________________________