ITS/AliITSMultReconstructor.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 // AliITSMultReconstructor - find clusters in the pixels (theta and
  21 // phi) and tracklets.
  22 //
  23 // These can be used to extract charged particles multiplcicity from the ITS.
  24 //
  25 // A tracklet consist of two ITS clusters, one in the first pixel
  26 // layer and one in the second. The clusters are associates if the
  27 // differencies in Phi (azimuth) and Zeta (longitudinal) are inside
  28 // a fiducial volume. In case of multiple candidates it is selected the
  29 // candidate with minimum distance in Phi.
  30 // The parameter AssociationChoice allows to control if two clusters
  31 // in layer 2 can be associated to the same cluster in layer 1 or not.
  32 // (TRUE means double associations exluded; default = TRUE)
  33 //
  34 // Two methods return the number of traklets and the number of clusters
  35 // in the first SPD layer (GetNTracklets GetNSingleClusters)
  36 //
  37 // -----------------------------------------------------------------
  38 //
  39 // NOTE: The cuts on phi and zeta depend on the interacting system (p-p
  40 //  or Pb-Pb). Please, check the file AliITSMultReconstructor.h and be
  41 //  sure that SetPhiWindow and SetZetaWindow are defined accordingly.
  42 //
  43 //  Author :  Tiziano Virgili
  44 //
  45 //  Recent updates (D. Elia, INFN Bari):
  46 //     - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE)
  47 //     - phi definition changed to ALICE convention (0,2*TMath::pi())
  48 //     - cluster coordinates taken with GetGlobalXYZ()
  49 //
  50 //  fGeometry removed
  51 //____________________________________________________________________
  52
  53 #include <TClonesArray.h>
  54 #include <TH1F.h>
  55 #include <TH2F.h>
  56 #include <TTree.h>
  57
  58 #include "AliITSMultReconstructor.h"
  59 #include "AliITSRecPoint.h"
  60 #include "AliITSgeom.h"
  61 #include "AliLog.h"
  62
  63 //____________________________________________________________________
  64 ClassImp(AliITSMultReconstructor)
  65
  66
  67 //____________________________________________________________________
  68 AliITSMultReconstructor::AliITSMultReconstructor():
  69 fClustersLay1(0),
  70 fClustersLay2(0),
  71 fTracklets(0),
  72 fSClusters(0),
  73 fAssociationFlag(0),
  74 fNClustersLay1(0),
  75 fNClustersLay2(0),
  76 fNTracklets(0),
  77 fNSingleCluster(0),
  78 fPhiWindow(0),
  79 fZetaWindow(0),
  80 fOnlyOneTrackletPerC2(0),
  81 fHistOn(0),
  82 fhClustersDPhiAcc(0),
  83 fhClustersDThetaAcc(0),
  84 fhClustersDZetaAcc(0),
  85 fhClustersDPhiAll(0),
  86 fhClustersDThetaAll(0),
  87 fhClustersDZetaAll(0),
  88 fhDPhiVsDThetaAll(0),
  89 fhDPhiVsDThetaAcc(0),
  90 fhDPhiVsDZetaAll(0),
  91 fhDPhiVsDZetaAcc(0),
  92 fhetaTracklets(0),
  93 fhphiTracklets(0),
  94 fhetaClustersLay1(0),
  95 fhphiClustersLay1(0){
  96   // Method to reconstruct the charged particles multiplicity with the
  97   // SPD (tracklets).
  98
  99
 100   SetHistOn();
 101   SetPhiWindow();
 102   SetZetaWindow();
 103   SetOnlyOneTrackletPerC2();
 104
 105   fClustersLay1       = new Float_t*[300000];
 106   fClustersLay2       = new Float_t*[300000];
 107   fTracklets          = new Float_t*[300000];
 108   fSClusters          = new Float_t*[300000];
 109   fAssociationFlag    = new Bool_t[300000];
 110
 111   for(Int_t i=0; i<300000; i++) {
 112     fClustersLay1[i]       = new Float_t[6];
 113     fClustersLay2[i]       = new Float_t[6];
 114     fTracklets[i]          = new Float_t[5];
 115     fSClusters[i]           = new Float_t[2];
 116     fAssociationFlag[i]    = kFALSE;
 117   }
 118
 119   // definition of histograms
 120   fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,0.,0.1);
 121   fhClustersDPhiAcc->SetDirectory(0);
 122   fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
 123   fhClustersDThetaAcc->SetDirectory(0);
 124   fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
 125   fhClustersDZetaAcc->SetDirectory(0);
 126
 127   fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
 128   fhDPhiVsDZetaAcc->SetDirectory(0);
 129   fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
 130   fhDPhiVsDThetaAcc->SetDirectory(0);
 131
 132   fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,0.0,0.5);
 133   fhClustersDPhiAll->SetDirectory(0);
 134   fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
 135   fhClustersDThetaAll->SetDirectory(0);
 136   fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
 137   fhClustersDZetaAll->SetDirectory(0);
 138
 139   fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
 140   fhDPhiVsDZetaAll->SetDirectory(0);
 141   fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
 142   fhDPhiVsDThetaAll->SetDirectory(0);
 143
 144   fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
 145   fhetaTracklets->SetDirectory(0);
 146   fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100, 0., 2*TMath::Pi());
 147   fhphiTracklets->SetDirectory(0);
 148   fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
 149   fhetaClustersLay1->SetDirectory(0);
 150   fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
 151   fhphiClustersLay1->SetDirectory(0);
 152 }
 153
 154 //______________________________________________________________________
 155 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
 156 fClustersLay1(mr.fClustersLay1),
 157 fClustersLay2(mr.fClustersLay2),
 158 fTracklets(mr.fTracklets),
 159 fSClusters(mr.fSClusters),
 160 fAssociationFlag(mr.fAssociationFlag),
 161 fNClustersLay1(mr.fNClustersLay1),
 162 fNClustersLay2(mr.fNClustersLay2),
 163 fNTracklets(mr.fNTracklets),
 164 fNSingleCluster(mr.fNSingleCluster),
 165 fPhiWindow(mr.fPhiWindow),
 166 fZetaWindow(mr.fZetaWindow),
 167 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
 168 fHistOn(mr.fHistOn),
 169 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
 170 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
 171 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
 172 fhClustersDPhiAll(mr.fhClustersDPhiAll),
 173 fhClustersDThetaAll(mr.fhClustersDThetaAll),
 174 fhClustersDZetaAll(mr.fhClustersDZetaAll),
 175 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
 176 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
 177 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
 178 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
 179 fhetaTracklets(mr.fhetaTracklets),
 180 fhphiTracklets(mr.fhphiTracklets),
 181 fhetaClustersLay1(mr.fhetaClustersLay1),
 182 fhphiClustersLay1(mr.fhphiClustersLay1) {
 183   // Copy constructor
 184
 185 }
 186
 187 //______________________________________________________________________
 188 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
 189   // Assignment operator
 190   this->~AliITSMultReconstructor();
 191   new(this) AliITSMultReconstructor(mr);
 192   return *this;
 193 }
 194
 195 //______________________________________________________________________
 196 AliITSMultReconstructor::~AliITSMultReconstructor(){
 197   // Destructor
 198
 199   // delete histograms
 200   delete fhClustersDPhiAcc;
 201   delete fhClustersDThetaAcc;
 202   delete fhClustersDZetaAcc;
 203   delete fhClustersDPhiAll;
 204   delete fhClustersDThetaAll;
 205   delete fhClustersDZetaAll;
 206   delete fhDPhiVsDThetaAll;
 207   delete fhDPhiVsDThetaAcc;
 208   delete fhDPhiVsDZetaAll;
 209   delete fhDPhiVsDZetaAcc;
 210   delete fhetaTracklets;
 211   delete fhphiTracklets;
 212   delete fhetaClustersLay1;
 213   delete fhphiClustersLay1;
 214
 215   // delete arrays
 216   for(Int_t i=0; i<300000; i++) {
 217     delete [] fClustersLay1[i];
 218     delete [] fClustersLay2[i];
 219     delete [] fTracklets[i];
 220     delete [] fSClusters[i];
 221   }
 222   delete [] fClustersLay1;
 223   delete [] fClustersLay2;
 224   delete [] fTracklets;
 225   delete [] fSClusters;
 226
 227   delete [] fAssociationFlag;
 228 }
 229
 230 //____________________________________________________________________
 231 void
 232 AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
 233   //
 234   // - calls LoadClusterArray that finds the position of the clusters
 235   //   (in global coord)
 236   // - convert the cluster coordinates to theta, phi (seen from the
 237   //   interaction vertex). The third coordinate is used for ....
 238   // - makes an array of tracklets
 239   //
 240   // After this method has been called, the clusters of the two layers
 241   // and the tracklets can be retrieved by calling the Get'er methods.
 242
 243   // reset counters
 244   fNClustersLay1 = 0;
 245   fNClustersLay2 = 0;
 246   fNTracklets = 0;
 247   fNSingleCluster = 0;
 248   // loading the clusters
 249   LoadClusterArrays(clusterTree);
 250
 251   // find the tracklets
 252   AliDebug(1,"Looking for tracklets... ");
 253
 254   //###########################################################
 255   // Loop on layer 1 : finding theta, phi and z
 256   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 257     Float_t x = fClustersLay1[iC1][0] - vtx[0];
 258     Float_t y = fClustersLay1[iC1][1] - vtx[1];
 259     Float_t z = fClustersLay1[iC1][2] - vtx[2];
 260
 261     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
 262                                TMath::Power(y,2) +
 263                                TMath::Power(z,2));
 264
 265     fClustersLay1[iC1][0] = TMath::ACos(z/r);                   // Store Theta
 266     fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 267     fClustersLay1[iC1][2] = z/r;                                // Store scaled z
 268     if (fHistOn) {
 269       Float_t eta=fClustersLay1[iC1][0];
 270       eta= TMath::Tan(eta/2.);
 271       eta=-TMath::Log(eta);
 272       fhetaClustersLay1->Fill(eta);
 273       fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
 274     }
 275   }
 276
 277   // Loop on layer 2 : finding theta, phi and r
 278   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 279     Float_t x = fClustersLay2[iC2][0] - vtx[0];
 280     Float_t y = fClustersLay2[iC2][1] - vtx[1];
 281     Float_t z = fClustersLay2[iC2][2] - vtx[2];
 282
 283     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
 284                                TMath::Power(y,2) +
 285                                TMath::Power(z,2));
 286
 287     fClustersLay2[iC2][0] = TMath::ACos(z/r);                   // Store Theta
 288     fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 289     fClustersLay2[iC2][2] = z;                                  // Store z
 290
 291  // this only needs to be initialized for the fNClustersLay2 first associations
 292     fAssociationFlag[iC2] = kFALSE;
 293   }
 294
 295   //###########################################################
 296   // Loop on layer 1
 297   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 298
 299     // reset of variables for multiple candidates
 300     Int_t  iC2WithBestDist = 0;     // reset
 301     Float_t distmin        = 100.;  // just to put a huge number!
 302     Float_t dPhimin        = 0.;  // Used for histograms only!
 303     Float_t dThetamin      = 0.;  // Used for histograms only!
 304     Float_t dZetamin       = 0.;  // Used for histograms only!
 305
 306     // Loop on layer 2
 307     for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 308
 309       // The following excludes double associations
 310       if (!fAssociationFlag[iC2]) {
 311
 312         // find the difference in angles
 313         Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
 314         Float_t dPhi   = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
 315         // take into account boundary condition
 316         if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
 317
 318         // find the difference in z (between linear projection from layer 1
 319         // and the actual point: Dzeta= z1/r1*r2 -z2)
 320         Float_t r2   = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
 321         Float_t dZeta  = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
 322
 323         if (fHistOn) {
 324           fhClustersDPhiAll->Fill(dPhi);
 325           fhClustersDThetaAll->Fill(dTheta);
 326           fhClustersDZetaAll->Fill(dZeta);
 327           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
 328           fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
 329         }
 330         // make "elliptical" cut in Phi and Zeta!
 331         Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
 332
 333         if (d>1) continue;
 334
 335         //look for the minimum distance: the minimum is in iC2WithBestDist
 336         if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
 337           distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
 338           dPhimin = dPhi;
 339           dThetamin = dTheta;
 340           dZetamin = dZeta;
 341           iC2WithBestDist = iC2;
 342         }
 343       }
 344     } // end of loop over clusters in layer 2
 345
 346     if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
 347
 348       if (fHistOn) {
 349         fhClustersDPhiAcc->Fill(dPhimin);
 350         fhClustersDThetaAcc->Fill(dThetamin);
 351         fhClustersDZetaAcc->Fill(dZetamin);
 352         fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
 353         fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
 354       }
 355
 356       if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
 357
 358       // store the tracklet
 359
 360       // use the theta from the clusters in the first layer
 361       fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
 362       // use the phi from the clusters in the first layer
 363       fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
 364       // Store the difference between phi1 and phi2
 365       fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
 366
 367       // find label
 368       // if equal label in both clusters found this label is assigned
 369       // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
 370       Int_t label1 = 0;
 371       Int_t label2 = 0;
 372       while (label2 < 3)
 373       {
 374         if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
 375           break;
 376         label1++;
 377         if (label1 == 3)
 378         {
 379           label1 = 0;
 380           label2++;
 381         }
 382       }
 383
 384       if (label2 < 3)
 385       {
 386         AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) fClustersLay1[iC1][3+label1], (Int_t) fClustersLay2[iC2WithBestDist][3+label2], fNTracklets));
 387         fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
 388         fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3+label2];
 389       }
 390       else
 391       {
 392         AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", (Int_t) fClustersLay1[iC1][3], (Int_t) fClustersLay1[iC1][4], (Int_t) fClustersLay1[iC1][5], (Int_t) fClustersLay2[iC2WithBestDist][3], (Int_t) fClustersLay2[iC2WithBestDist][4], (Int_t) fClustersLay2[iC2WithBestDist][5], fNTracklets));
 393         fTracklets[fNTracklets][3] = fClustersLay1[iC1][3];
 394         fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3];
 395       }
 396
 397       if (fHistOn) {
 398         Float_t eta=fTracklets[fNTracklets][0];
 399         eta= TMath::Tan(eta/2.);
 400         eta=-TMath::Log(eta);
 401         fhetaTracklets->Fill(eta);
 402         fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
 403       }
 404
 405       AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
 406       AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1,
 407                       iC2WithBestDist));
 408       fNTracklets++;
 409     }
 410
 411     // Delete the following else if you do not want to save Clusters!
 412
 413     else { // This means that the cluster has not been associated
 414
 415       // store the cluster
 416
 417       fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
 418       fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
 419       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
 420                       fNSingleCluster, iC1));
 421       fNSingleCluster++;
 422     }
 423
 424   } // end of loop over clusters in layer 1
 425
 426   AliDebug(1,Form("%d tracklets found", fNTracklets));
 427 }
 428
 429 //____________________________________________________________________
 430 void
 431 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
 432   // This method
 433   // - gets the clusters from the cluster tree
 434   // - convert them into global coordinates
 435   // - store them in the internal arrays
 436
 437   AliDebug(1,"Loading clusters ...");
 438
 439   fNClustersLay1 = 0;
 440   fNClustersLay2 = 0;
 441
 442   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
 443   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
 444
 445   itsClusterBranch->SetAddress(&itsClusters);
 446
 447   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
 448   Float_t cluGlo[3]={0.,0.,0.};
 449
 450   // loop over the its subdetectors
 451   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
 452
 453     if (!itsClusterTree->GetEvent(iIts))
 454       continue;
 455
 456     Int_t nClusters = itsClusters->GetEntriesFast();
 457
 458     // loop over clusters
 459     while(nClusters--) {
 460       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 461
 462       if (cluster->GetLayer()>1)
 463         continue;
 464
 465       cluster->GetGlobalXYZ(cluGlo);
 466       Float_t x = cluGlo[0];
 467       Float_t y = cluGlo[1];
 468       Float_t z = cluGlo[2];
 469
 470       if (cluster->GetLayer()==0) {
 471         fClustersLay1[fNClustersLay1][0] = x;
 472         fClustersLay1[fNClustersLay1][1] = y;
 473         fClustersLay1[fNClustersLay1][2] = z;
 474         for (Int_t i=0; i<3; i++)
 475                 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
 476         fNClustersLay1++;
 477       }
 478       if (cluster->GetLayer()==1) {
 479         fClustersLay2[fNClustersLay2][0] = x;
 480         fClustersLay2[fNClustersLay2][1] = y;
 481         fClustersLay2[fNClustersLay2][2] = z;
 482         for (Int_t i=0; i<3; i++)
 483                 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
 484         fNClustersLay2++;
 485       }
 486
 487     }// end of cluster loop
 488   } // end of its "subdetector" loop
 489
 490   if (itsClusters) {
 491     itsClusters->Delete();
 492     delete itsClusters;
 493     itsClusters = 0;
 494   }
 495   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
 496 }
 497 //____________________________________________________________________
 498 void
 499 AliITSMultReconstructor::SaveHists() {
 500   // This method save the histograms on the output file
 501   // (only if fHistOn is TRUE).
 502
 503   if (!fHistOn)
 504     return;
 505
 506   fhClustersDPhiAll->Write();
 507   fhClustersDThetaAll->Write();
 508   fhClustersDZetaAll->Write();
 509   fhDPhiVsDThetaAll->Write();
 510   fhDPhiVsDZetaAll->Write();
 511
 512   fhClustersDPhiAcc->Write();
 513   fhClustersDThetaAcc->Write();
 514   fhClustersDZetaAcc->Write();
 515   fhDPhiVsDThetaAcc->Write();
 516   fhDPhiVsDZetaAcc->Write();
 517
 518   fhetaTracklets->Write();
 519   fhphiTracklets->Write();
 520   fhetaClustersLay1->Write();
 521   fhphiClustersLay1->Write();
 522 }
 523
 524