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