ITS/AliITSMultReconstructor.cxx

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