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