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