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 //_________________________________________________________________________
  17 //
  18 //        Implementation of the ITS-SPD trackleter class
  19 //
  20 // It retrieves clusters in the pixels (theta and phi) and finds tracklets.
  21 // These can be used to extract charged particle multiplicity from the ITS.
  22 //
  23 // A tracklet consists of two ITS clusters, one in the first pixel layer and
  24 // one in the second. The clusters are associated if the differences in
  25 // Phi (azimuth) and Theta (polar angle) are within fiducial windows.
  26 // In case of multiple candidates the candidate with minimum
  27 // distance is selected.
  28 //
  29 // Two methods return the number of tracklets and the number of unassociated
  30 // clusters (i.e. not used in any tracklet) in the first SPD layer
  31 // (GetNTracklets and GetNSingleClusters)
  32 //
  33 // The cuts on phi and theta depend on the interacting system (p-p or Pb-Pb)
  34 // and can be set via AliITSRecoParam class
  35 // (SetPhiWindow and SetThetaWindow)
  36 //
  37 // Origin: Tiziano Virgili
  38 //
  39 // Current support and development:
  40 //         Domenico Elia, Maria Nicassio (INFN Bari)
  41 //         Domenico.Elia@ba.infn.it, Maria.Nicassio@ba.infn.it
  42 //
  43 // Most recent updates:
  44 //     - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE)
  45 //     - phi definition changed to ALICE convention (0,2*TMath::pi())
  46 //     - cluster coordinates taken with GetGlobalXYZ()
  47 //     - fGeometry removed
  48 //     - number of fired chips on the two layers
  49 //     - option to cut duplicates in the overlaps
  50 //     - options and fiducial cuts via AliITSRecoParam
  51 //     - move from DeltaZeta to DeltaTheta cut
  52 //     - update to the new algorithm by Mariella and Jan Fiete
  53 //     - store also DeltaTheta in the ESD
  54 //     - less new and delete calls when creating the needed arrays
  55 //
  56 //     - RS: to decrease the number of new/deletes the clusters data are stored
  57 //           not in float[6] attached to float**, but in 1-D array.
  58 //     - RS: Clusters are sorted in Z in roder to have the same numbering as in the ITS reco
  59 //     - RS: Clusters used by ESDtrack are flagged, this information is passed to AliMulitiplicity object
  60 //           when storing the tracklets and single cluster info
  61 //_________________________________________________________________________
  62
  63 #include <TClonesArray.h>
  64 #include <TH1F.h>
  65 #include <TH2F.h>
  66 #include <TTree.h>
  67 #include <TBits.h>
  68 #include <TArrayI.h>
  69
  70 #include "AliITSMultReconstructor.h"
  71 #include "AliITSReconstructor.h"
  72 #include "AliITSsegmentationSPD.h"
  73 #include "AliITSRecPoint.h"
  74 #include "AliITSRecPointContainer.h"
  75 #include "AliITSgeom.h"
  76 #include "AliITSgeomTGeo.h"
  77 #include "AliITSDetTypeRec.h"
  78 #include "AliESDEvent.h"
  79 #include "AliESDVertex.h"
  80 #include "AliESDtrack.h"
  81 #include "AliMultiplicity.h"
  82 #include "AliLog.h"
  83 #include "TGeoGlobalMagField.h"
  84 #include "AliMagF.h"
  85
  86 //____________________________________________________________________
  87 ClassImp(AliITSMultReconstructor)
  88
  89
  90 //____________________________________________________________________
  91 AliITSMultReconstructor::AliITSMultReconstructor():
  92 fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fUsedClusLay1(0),fUsedClusLay2(0),
  93 fClustersLay1(0),
  94 fClustersLay2(0),
  95 fDetectorIndexClustersLay1(0),
  96 fDetectorIndexClustersLay2(0),
  97 fOverlapFlagClustersLay1(0),
  98 fOverlapFlagClustersLay2(0),
  99 fTracklets(0),
 100 fSClusters(0),
 101 fNClustersLay1(0),
 102 fNClustersLay2(0),
 103 fNTracklets(0),
 104 fNSingleCluster(0),
 105 fPhiWindow(0),
 106 fThetaWindow(0),
 107 fPhiShift(0),
 108 fRemoveClustersFromOverlaps(0),
 109 fPhiOverlapCut(0),
 110 fZetaOverlapCut(0),
 111 fHistOn(0),
 112 fhClustersDPhiAcc(0),
 113 fhClustersDThetaAcc(0),
 114 fhClustersDPhiAll(0),
 115 fhClustersDThetaAll(0),
 116 fhDPhiVsDThetaAll(0),
 117 fhDPhiVsDThetaAcc(0),
 118 fhetaTracklets(0),
 119 fhphiTracklets(0),
 120 fhetaClustersLay1(0),
 121 fhphiClustersLay1(0){
 122
 123   fNFiredChips[0] = 0;
 124   fNFiredChips[1] = 0;
 125   // Method to reconstruct the charged particles multiplicity with the
 126   // SPD (tracklets).
 127
 128   SetHistOn();
 129
 130   if(AliITSReconstructor::GetRecoParam()) {
 131     SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow());
 132     SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow());
 133     SetPhiShift(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiShift());
 134     SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps());
 135     SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut());
 136     SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut());
 137   } else {
 138     SetPhiWindow();
 139     SetThetaWindow();
 140     SetPhiShift();
 141     SetRemoveClustersFromOverlaps();
 142     SetPhiOverlapCut();
 143     SetZetaOverlapCut();
 144   }
 145
 146   fClustersLay1              = 0;
 147   fClustersLay2              = 0;
 148   fDetectorIndexClustersLay1 = 0;
 149   fDetectorIndexClustersLay2 = 0;
 150   fOverlapFlagClustersLay1   = 0;
 151   fOverlapFlagClustersLay2   = 0;
 152   fTracklets                 = 0;
 153   fSClusters                 = 0;
 154
 155   // definition of histograms
 156   Bool_t oldStatus = TH1::AddDirectoryStatus();
 157   TH1::AddDirectory(kFALSE);
 158
 159   fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,-0.1,0.1);
 160   fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
 161
 162   fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
 163
 164   fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,0.0,0.5);
 165   fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,0.0,0.5);
 166
 167   fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,0.,0.5,100,0.,0.5);
 168
 169   fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
 170   fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100, 0., 2*TMath::Pi());
 171   fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
 172   fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
 173
 174   TH1::AddDirectory(oldStatus);
 175 }
 176
 177 //______________________________________________________________________
 178 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) :
 179 AliTrackleter(mr),
 180 fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fUsedClusLay1(0),fUsedClusLay2(0),
 181 fClustersLay1(0),
 182 fClustersLay2(0),
 183 fDetectorIndexClustersLay1(0),
 184 fDetectorIndexClustersLay2(0),
 185 fOverlapFlagClustersLay1(0),
 186 fOverlapFlagClustersLay2(0),
 187 fTracklets(0),
 188 fSClusters(0),
 189 fNClustersLay1(0),
 190 fNClustersLay2(0),
 191 fNTracklets(0),
 192 fNSingleCluster(0),
 193 fPhiWindow(0),
 194 fThetaWindow(0),
 195 fPhiShift(0),
 196 fRemoveClustersFromOverlaps(0),
 197 fPhiOverlapCut(0),
 198 fZetaOverlapCut(0),
 199 fHistOn(0),
 200 fhClustersDPhiAcc(0),
 201 fhClustersDThetaAcc(0),
 202 fhClustersDPhiAll(0),
 203 fhClustersDThetaAll(0),
 204 fhDPhiVsDThetaAll(0),
 205 fhDPhiVsDThetaAcc(0),
 206 fhetaTracklets(0),
 207 fhphiTracklets(0),
 208 fhetaClustersLay1(0),
 209 fhphiClustersLay1(0)
 210  {
 211   // Copy constructor :!!! RS ATTENTION: old c-tor reassigned the pointers instead of creating a new copy -> would crash on delete
 212    AliError("May not use");
 213 }
 214
 215 //______________________________________________________________________
 216 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
 217   // Assignment operator
 218   if (this != &mr) {
 219     this->~AliITSMultReconstructor();
 220     new(this) AliITSMultReconstructor(mr);
 221   }
 222   return *this;
 223 }
 224
 225 //______________________________________________________________________
 226 AliITSMultReconstructor::~AliITSMultReconstructor(){
 227   // Destructor
 228
 229   // delete histograms
 230   delete fhClustersDPhiAcc;
 231   delete fhClustersDThetaAcc;
 232   delete fhClustersDPhiAll;
 233   delete fhClustersDThetaAll;
 234   delete fhDPhiVsDThetaAll;
 235   delete fhDPhiVsDThetaAcc;
 236   delete fhetaTracklets;
 237   delete fhphiTracklets;
 238   delete fhetaClustersLay1;
 239   delete fhphiClustersLay1;
 240   delete[] fUsedClusLay1;
 241   delete[] fUsedClusLay2;
 242   // delete arrays
 243   for(Int_t i=0; i<fNTracklets; i++)
 244     delete [] fTracklets[i];
 245
 246   for(Int_t i=0; i<fNSingleCluster; i++)
 247     delete [] fSClusters[i];
 248
 249   delete [] fClustersLay1;
 250   delete [] fClustersLay2;
 251   delete [] fDetectorIndexClustersLay1;
 252   delete [] fDetectorIndexClustersLay2;
 253   delete [] fOverlapFlagClustersLay1;
 254   delete [] fOverlapFlagClustersLay2;
 255   delete [] fTracklets;
 256   delete [] fSClusters;
 257 }
 258
 259 //____________________________________________________________________
 260 void AliITSMultReconstructor::Reconstruct(AliESDEvent* esd, TTree* treeRP)
 261 {
 262   //
 263   // - calls LoadClusterArrays that finds the position of the clusters
 264   //   (in global coord)
 265   // - convert the cluster coordinates to theta, phi (seen from the
 266   //   interaction vertex).
 267   // - makes an array of tracklets
 268   //
 269   // After this method has been called, the clusters of the two layers
 270   // and the tracklets can be retrieved by calling the Get'er methods.
 271
 272   if (!treeRP) { AliError(" Invalid ITS cluster tree !\n"); return; }
 273   if (!esd) {AliError("ESDEvent is not available, use old reconstructor"); return;}
 274   // reset counters
 275   if (fMult) delete fMult; fMult = 0;
 276   fNClustersLay1 = 0;
 277   fNClustersLay2 = 0;
 278   fNTracklets = 0;
 279   fNSingleCluster = 0;
 280   //
 281   fESDEvent = esd;
 282   fTreeRP = treeRP;
 283   //
 284   // >>>> RS: this part is equivalent to former AliITSVertexer::FindMultiplicity
 285   //
 286   // see if there is a SPD vertex
 287   Bool_t isVtxOK=kTRUE, isCosmics=kFALSE;
 288   AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
 289   if (!vtx && vtx->GetNContributors()<0) isVtxOK = kFALSE;
 290   if (vtx && strstr(vtx->GetTitle(),"cosmics")) {
 291     isVtxOK = kFALSE;
 292     isCosmics = kTRUE;
 293   }
 294   //
 295   if (!isVtxOK) {
 296     if (!isCosmics) {
 297       AliDebug(1,"Tracklets multiplicity not determined because the primary vertex was not found");
 298       AliDebug(1,"Just counting the number of cluster-fired chips on the SPD layers");
 299     }
 300     vtx = 0;
 301   }
 302   float vtxf[3] = {vtx->GetX(),vtx->GetY(),vtx->GetZ()};
 303   FindTracklets(vtxf);
 304   //
 305   CreateMultiplicityObject();
 306 }
 307
 308 //____________________________________________________________________
 309 void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
 310   //
 311   // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder
 312   //
 313   // - calls LoadClusterArray that finds the position of the clusters
 314   //   (in global coord)
 315   // - convert the cluster coordinates to theta, phi (seen from the
 316   //   interaction vertex).
 317   // - makes an array of tracklets
 318   //
 319   // After this method has been called, the clusters of the two layers
 320   // and the tracklets can be retrieved by calling the Get'er methods.
 321   if (fMult) delete fMult; fMult = 0;
 322   fNClustersLay1 = 0;
 323   fNClustersLay2 = 0;
 324   fNTracklets = 0;
 325   fNSingleCluster = 0;
 326   //
 327   if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
 328   //
 329   fESDEvent = 0;
 330   fTreeRP = clusterTree;
 331   //
 332   FindTracklets(vtx);
 333   //
 334 }
 335
 336 //____________________________________________________________________
 337 void AliITSMultReconstructor::FindTracklets(const Float_t *vtx)
 338 {
 339   // Find tracklets converging to vertex
 340   //
 341   LoadClusterArrays(fTreeRP);
 342   // flag clusters used by ESD tracks
 343   if (fESDEvent) ProcessESDTracks();
 344
 345   if (!vtx) return;
 346
 347   const Double_t pi = TMath::Pi();
 348
 349   // dPhi shift is field dependent
 350   // get average magnetic field
 351   Float_t bz = 0;
 352   AliMagF* field = 0;
 353   if (TGeoGlobalMagField::Instance()) field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
 354   if (!field)
 355   {
 356     AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor.")
 357   }
 358   else
 359     bz = TMath::Abs(field->SolenoidField());
 360
 361   const Double_t dPhiShift = fPhiShift / 5 * bz;
 362   AliDebug(1, Form("Using phi shift of %f", dPhiShift));
 363
 364   const Double_t dPhiWindow2 = fPhiWindow * fPhiWindow;
 365   const Double_t dThetaWindow2 = fThetaWindow * fThetaWindow;
 366
 367   Int_t* partners = new Int_t[fNClustersLay2];
 368   Float_t* minDists = new Float_t[fNClustersLay2];
 369   Int_t* associatedLay1 = new Int_t[fNClustersLay1];
 370   TArrayI** blacklist = new TArrayI*[fNClustersLay1];
 371
 372   for (Int_t i=0; i<fNClustersLay2; i++) {
 373     partners[i] = -1;
 374     minDists[i] = 2;
 375   }
 376   for (Int_t i=0; i<fNClustersLay1; i++)
 377     associatedLay1[i] = 0;
 378   for (Int_t i=0; i<fNClustersLay1; i++)
 379     blacklist[i] = 0;
 380
 381   // find the tracklets
 382   AliDebug(1,"Looking for tracklets... ");
 383
 384   //###########################################################
 385   // Loop on layer 1 : finding theta, phi and z
 386   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 387     float *clPar = GetClusterLayer1(iC1);
 388     Float_t x = clPar[kClTh] - vtx[0];
 389     Float_t y = clPar[kClPh] - vtx[1];
 390     Float_t z = clPar[kClZ]  - vtx[2];
 391
 392     Float_t r    = TMath::Sqrt(x*x + y*y + z*z);
 393
 394     clPar[kClTh] = TMath::ACos(z/r);                   // Store Theta
 395     clPar[kClPh] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 396
 397     if (fHistOn) {
 398       Float_t eta = clPar[kClTh];
 399       eta= TMath::Tan(eta/2.);
 400       eta=-TMath::Log(eta);
 401       fhetaClustersLay1->Fill(eta);
 402       fhphiClustersLay1->Fill(clPar[kClPh]);
 403     }
 404   }
 405
 406   // Loop on layer 2 : finding theta, phi and r
 407   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 408     float *clPar = GetClusterLayer2(iC2);
 409     Float_t x = clPar[kClTh] - vtx[0];
 410     Float_t y = clPar[kClPh] - vtx[1];
 411     Float_t z = clPar[kClZ]  - vtx[2];
 412
 413     Float_t r    = TMath::Sqrt(x*x + y*y + z*z);
 414
 415     clPar[kClTh] = TMath::ACos(z/r);                   // Store Theta
 416     clPar[kClPh] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 417   }
 418
 419   //###########################################################
 420   Int_t found = 1;
 421   while (found > 0) {
 422     found = 0;
 423
 424     // Step1: find all tracklets allowing double assocation
 425     // Loop on layer 1
 426     for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 427
 428       // already used or in the overlap ?
 429       if (associatedLay1[iC1] != 0 || fOverlapFlagClustersLay1[iC1]) continue;
 430
 431       found++;
 432
 433       // reset of variables for multiple candidates
 434       Int_t  iC2WithBestDist = -1;   // reset
 435       Double_t minDist       =  2;   // reset
 436       float* clPar1 = GetClusterLayer1(iC1);
 437
 438       // Loop on layer 2
 439       for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 440
 441         // in the overlap ?
 442         if (fOverlapFlagClustersLay2[iC2]) continue;
 443         float* clPar2 = GetClusterLayer2(iC2);
 444
 445         if (blacklist[iC1]) {
 446           Bool_t blacklisted = kFALSE;
 447           for (Int_t i=blacklist[iC1]->GetSize(); i--;) {
 448             if (blacklist[iC1]->At(i) == iC2) {
 449               blacklisted = kTRUE;
 450               break;
 451             }
 452           }
 453           if (blacklisted) continue;
 454         }
 455
 456         // find the difference in angles
 457         Double_t dTheta = TMath::Abs(clPar2[kClTh] - clPar1[kClTh]);
 458         Double_t dPhi   = TMath::Abs(clPar2[kClPh] - clPar1[kClPh]);
 459         // take into account boundary condition
 460         if (dPhi>pi) dPhi=2.*pi-dPhi;
 461
 462         if (fHistOn) {
 463           fhClustersDPhiAll->Fill(dPhi);
 464           fhClustersDThetaAll->Fill(dTheta);
 465           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
 466         }
 467
 468         dPhi -= dPhiShift;
 469
 470         // make "elliptical" cut in Phi and Theta!
 471         Float_t d = dPhi*dPhi/dPhiWindow2 + dTheta*dTheta/dThetaWindow2;
 472
 473         // look for the minimum distance: the minimum is in iC2WithBestDist
 474         if (d<1 && d<minDist) {
 475           minDist=d;
 476           iC2WithBestDist = iC2;
 477         }
 478       } // end of loop over clusters in layer 2
 479
 480       if (minDist<1) { // This means that a cluster in layer 2 was found that matches with iC1
 481
 482         if (minDists[iC2WithBestDist] > minDist) {
 483           Int_t oldPartner = partners[iC2WithBestDist];
 484           partners[iC2WithBestDist] = iC1;
 485           minDists[iC2WithBestDist] = minDist;
 486
 487           // mark as assigned
 488           associatedLay1[iC1] = 1;
 489
 490           if (oldPartner != -1) {
 491             // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its blacklist
 492             if (blacklist[oldPartner] == 0) {
 493               blacklist[oldPartner] = new TArrayI(1);
 494             } else blacklist[oldPartner]->Set(blacklist[oldPartner]->GetSize()+1);
 495
 496             blacklist[oldPartner]->AddAt(iC2WithBestDist, blacklist[oldPartner]->GetSize()-1);
 497
 498             // mark as free
 499             associatedLay1[oldPartner] = 0;
 500           }
 501         } else {
 502           // try again to find a cluster without considering iC2WithBestDist
 503           if (blacklist[iC1] == 0) {
 504             blacklist[iC1] = new TArrayI(1);
 505           }
 506           else
 507             blacklist[iC1]->Set(blacklist[iC1]->GetSize()+1);
 508
 509           blacklist[iC1]->AddAt(iC2WithBestDist, blacklist[iC1]->GetSize()-1);
 510         }
 511
 512       } else // cluster has no partner; remove
 513         associatedLay1[iC1] = 2;
 514     } // end of loop over clusters in layer 1
 515   }
 516
 517   // Step2: store tracklets; remove used clusters
 518   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 519
 520     if (partners[iC2] == -1) continue;
 521
 522     if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (partners[iC2],iC2);
 523
 524
 525     if (fOverlapFlagClustersLay1[partners[iC2]] || fOverlapFlagClustersLay2[iC2]) continue;
 526
 527     float* clPar2 = GetClusterLayer2(iC2);
 528     float* clPar1 = GetClusterLayer1(partners[iC2]);
 529
 530     Float_t* tracklet = fTracklets[fNTracklets] = new Float_t[kTrNPar]; // RS Add also the cluster id's
 531
 532     // use the theta from the clusters in the first layer
 533     tracklet[kTrTheta] = clPar1[kClTh];
 534     // use the phi from the clusters in the first layer
 535     tracklet[kTrPhi] = clPar1[kClPh];
 536     // store the difference between phi1 and phi2
 537     tracklet[kTrDPhi] = clPar1[kClPh] - clPar2[kClPh];
 538
 539     // define dphi in the range [0,pi] with proper sign (track charge correlated)
 540     if (tracklet[kTrDPhi] > TMath::Pi())   tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi();
 541     if (tracklet[kTrDPhi] < -TMath::Pi())  tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi();
 542
 543     // store the difference between theta1 and theta2
 544     tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh];
 545
 546     if (fHistOn) {
 547       fhClustersDPhiAcc->Fill(tracklet[kTrDPhi]);
 548       fhClustersDThetaAcc->Fill(tracklet[kTrDTheta]);
 549       fhDPhiVsDThetaAcc->Fill(tracklet[kTrDTheta],tracklet[kTrDPhi]);
 550     }
 551
 552     // find label
 553     // if equal label in both clusters found this label is assigned
 554     // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
 555     Int_t label1 = 0;
 556     Int_t label2 = 0;
 557     while (label2 < 3) {
 558       if ((Int_t) clPar1[kClMC0+label1] != -2 && (Int_t) clPar1[kClMC0+label1] == (Int_t) clPar2[kClMC0+label2])
 559         break;
 560       label1++;
 561       if (label1 == 3) {
 562         label1 = 0;
 563         label2++;
 564       }
 565     }
 566     if (label2 < 3) {
 567       AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) clPar1[kClMC0+label1], (Int_t) clPar1[kClMC0+label2], fNTracklets));
 568       tracklet[kTrLab1] = clPar1[kClMC0+label1];
 569       tracklet[kTrLab2] = clPar2[kClMC0+label2];
 570     } else {
 571       AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", (Int_t) clPar1[kClMC0], (Int_t) clPar1[kClMC1], (Int_t) clPar1[kClMC2], (Int_t) clPar2[kClMC0], (Int_t) clPar2[kClMC1], (Int_t) clPar2[kClMC2], fNTracklets));
 572       tracklet[kTrLab1] = clPar1[kClMC0];
 573       tracklet[kTrLab2] = clPar2[kClMC0];
 574     }
 575
 576     if (fHistOn) {
 577       Float_t eta = tracklet[kTrTheta];
 578       eta= TMath::Tan(eta/2.);
 579       eta=-TMath::Log(eta);
 580       fhetaTracklets->Fill(eta);
 581       fhphiTracklets->Fill(tracklet[kTrPhi]);
 582     }
 583     //
 584     tracklet[kClID1] = partners[iC2];
 585     tracklet[kClID2] = iC2;
 586     //
 587     AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
 588     AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", partners[iC2], iC2));
 589     fNTracklets++;
 590
 591     associatedLay1[partners[iC2]] = 1;
 592   }
 593
 594   // Delete the following else if you do not want to save Clusters!
 595   // store the cluster
 596   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 597
 598     float* clPar1 = GetClusterLayer1(iC1);
 599
 600     if (associatedLay1[iC1]==2||associatedLay1[iC1]==0) {
 601       fSClusters[fNSingleCluster] = new Float_t[kClNPar];
 602       fSClusters[fNSingleCluster][kSCTh] = clPar1[kClTh];
 603       fSClusters[fNSingleCluster][kSCPh] = clPar1[kClPh];
 604       fSClusters[fNSingleCluster][kSCID] = iC1;
 605       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
 606                 fNSingleCluster, iC1));
 607       fNSingleCluster++;
 608     }
 609   }
 610
 611   delete[] partners;
 612   delete[] minDists;
 613
 614   for (Int_t i=0; i<fNClustersLay1; i++)
 615     if (blacklist[i])
 616       delete blacklist[i];
 617   delete[] blacklist;
 618
 619   AliDebug(1,Form("%d tracklets found", fNTracklets));
 620 }
 621
 622 //____________________________________________________________________
 623 void AliITSMultReconstructor::CreateMultiplicityObject()
 624 {
 625   // create AliMultiplicity object and store it in the ESD event
 626   //
 627   TBits fastOrFiredMap,firedChipMap;
 628   if (fDetTypeRec) {
 629    fastOrFiredMap  = fDetTypeRec->GetFastOrFiredMap();
 630    firedChipMap    = fDetTypeRec->GetFiredChipMap(fTreeRP);
 631   }
 632   //
 633   fMult = new AliMultiplicity(fNTracklets,fNSingleCluster,fNFiredChips[0],fNFiredChips[1],fastOrFiredMap);
 634   fMult->SetFiredChipMap(firedChipMap);
 635   AliITSRecPointContainer* rcont = AliITSRecPointContainer::Instance();
 636   fMult->SetITSClusters(0,rcont->GetNClustersInLayer(1,fTreeRP));
 637   for(Int_t kk=2;kk<=6;kk++) fMult->SetITSClusters(kk-1,rcont->GetNClustersInLayerFast(kk));
 638   //
 639   for (int i=fNTracklets;i--;)  {
 640     float* tlInfo = fTracklets[i];
 641     fMult->SetTrackletData(i,tlInfo, fUsedClusLay1[int(tlInfo[kClID1])]|fUsedClusLay2[int(tlInfo[kClID2])]);
 642   }
 643   //
 644   for (int i=fNSingleCluster;i--;) {
 645     float* clInfo = fSClusters[i];
 646     fMult->SetSingleClusterData(i,clInfo,fUsedClusLay1[int(clInfo[kSCID])]);
 647   }
 648   fMult->CompactBits();
 649   //
 650 }
 651
 652
 653 //____________________________________________________________________
 654 void AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree)
 655 {
 656   // This method
 657   // - gets the clusters from the cluster tree
 658   // - convert them into global coordinates
 659   // - store them in the internal arrays
 660   // - count the number of cluster-fired chips
 661   //
 662   // RS: This method was strongly modified wrt original by Jan Fiete. In order to have the same numbering
 663   // of clusters as in the ITS reco I had to introduce sorting in Z
 664   // Also note that now the clusters data are stored not in float[6] attached to float**, but in 1-D array
 665
 666   AliDebug(1,"Loading clusters and cluster-fired chips ...");
 667
 668   fNClustersLay1 = 0;
 669   fNClustersLay2 = 0;
 670   fNFiredChips[0] = 0;
 671   fNFiredChips[1] = 0;
 672
 673   AliITSsegmentationSPD seg;
 674
 675   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 676   TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
 677   if(!rpcont->IsSPDActive()){
 678     AliWarning("No SPD rec points found, multiplicity not calculated");
 679     return;
 680   }
 681   //
 682   // count clusters
 683   // loop over the SPD subdetectors
 684   TObjArray clArr(100);
 685   for (int il=0;il<2;il++) {
 686     int nclLayer = 0;
 687     int detMin = AliITSgeomTGeo::GetModuleIndex(il+1,1,1);
 688     int detMax = AliITSgeomTGeo::GetModuleIndex(il+2,1,1);
 689     for (int idt=detMin;idt<detMax;idt++) {
 690       itsClusters=rpcont->UncheckedGetClusters(idt);
 691       int nClusters = itsClusters->GetEntriesFast();
 692       if (!nClusters) continue;
 693       Int_t nClustersInChip[5] = {0,0,0,0,0};
 694       while(nClusters--) {
 695         AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 696         if (!cluster) continue;
 697         clArr.AddAtAndExpand(cluster,nclLayer++);
 698         nClustersInChip[ seg.GetChipFromLocal(0,cluster->GetDetLocalZ()) ]++;
 699       }
 700       for(Int_t ifChip=5;ifChip--;) if (nClustersInChip[ifChip]) fNFiredChips[il]++;
 701     }
 702     // sort the clusters in Z (to have the same numbering as in ITS reco
 703     Float_t *z    = new Float_t[nclLayer];
 704     Int_t * index = new Int_t[nclLayer];
 705     for (int ic=0;ic<nclLayer;ic++) z[ic] = ((AliITSRecPoint*)clArr[ic])->GetZ();
 706     TMath::Sort(nclLayer,z,index,kFALSE);
 707     Float_t*   clustersLay              = new Float_t[nclLayer*kClNPar];
 708     Int_t*     detectorIndexClustersLay = new Int_t[nclLayer];
 709     Bool_t*    overlapFlagClustersLay   = new Bool_t[nclLayer];
 710     Char_t*    usedClusLay              = new Char_t[nclLayer];
 711     //
 712     for (int ic=0;ic<nclLayer;ic++) {
 713       AliITSRecPoint* cluster = (AliITSRecPoint*)clArr[index[ic]];
 714       float* clPar = &clustersLay[ic*kClNPar];
 715       //
 716       cluster->GetGlobalXYZ( clPar );
 717       detectorIndexClustersLay[ic] = cluster->GetDetectorIndex();
 718       overlapFlagClustersLay[ic]   = kFALSE;
 719       usedClusLay[ic]              = 0;
 720       for (Int_t i=3;i--;) clPar[kClMC0+i] = cluster->GetLabel(i);
 721     }
 722     clArr.Clear();
 723     delete[] z;
 724     delete[] index;
 725     //
 726     if (il==0) {
 727       fClustersLay1              = clustersLay;
 728       fOverlapFlagClustersLay1   = overlapFlagClustersLay;
 729       fDetectorIndexClustersLay1 = detectorIndexClustersLay;
 730       fUsedClusLay1              = usedClusLay;
 731       fNClustersLay1             = nclLayer;
 732     }
 733     else {
 734       fClustersLay2              = clustersLay;
 735       fOverlapFlagClustersLay2   = overlapFlagClustersLay;
 736       fDetectorIndexClustersLay2 = detectorIndexClustersLay;
 737       fUsedClusLay2              = usedClusLay;
 738       fNClustersLay2             = nclLayer;
 739     }
 740   }
 741   //
 742   // no double association allowed
 743   int nmaxT                  = TMath::Min(fNClustersLay1, fNClustersLay2);
 744   fTracklets                 = new Float_t*[nmaxT];
 745   fSClusters                 = new Float_t*[fNClustersLay1];
 746   for (Int_t i=nmaxT;i--;) fTracklets[i] = 0;
 747   //
 748   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
 749   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 750 }
 751 //____________________________________________________________________
 752 void
 753 AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
 754   // This method
 755   // - gets the clusters from the cluster tree
 756   // - counts the number of (cluster)fired chips
 757
 758   AliDebug(1,"Loading cluster-fired chips ...");
 759
 760   fNFiredChips[0] = 0;
 761   fNFiredChips[1] = 0;
 762
 763   AliITSsegmentationSPD seg;
 764   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 765   TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
 766   if(!rpcont->IsSPDActive()){
 767     AliWarning("No SPD rec points found, multiplicity not calculated");
 768     return;
 769   }
 770
 771   // loop over the its subdetectors
 772   Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1);
 773   for (Int_t iIts=0; iIts < nSPDmodules; iIts++) {
 774     itsClusters=rpcont->UncheckedGetClusters(iIts);
 775     Int_t nClusters = itsClusters->GetEntriesFast();
 776
 777     // number of clusters in each chip of the current module
 778     Int_t nClustersInChip[5] = {0,0,0,0,0};
 779     Int_t layer = 0;
 780
 781     // loop over clusters
 782     while(nClusters--) {
 783       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 784
 785       layer = cluster->GetLayer();
 786       if (layer>1) continue;
 787
 788       // find the chip for the current cluster
 789       Float_t locz = cluster->GetDetLocalZ();
 790       Int_t iChip = seg.GetChipFromLocal(0,locz);
 791       nClustersInChip[iChip]++;
 792
 793     }// end of cluster loop
 794
 795     // get number of fired chips in the current module
 796     for(Int_t ifChip=0; ifChip<5; ifChip++) {
 797       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
 798     }
 799
 800   } // end of its "subdetector" loop
 801
 802
 803   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 804 }
 805 //____________________________________________________________________
 806 void
 807 AliITSMultReconstructor::SaveHists() {
 808   // This method save the histograms on the output file
 809   // (only if fHistOn is TRUE).
 810
 811   if (!fHistOn)
 812     return;
 813
 814   fhClustersDPhiAll->Write();
 815   fhClustersDThetaAll->Write();
 816   fhDPhiVsDThetaAll->Write();
 817
 818   fhClustersDPhiAcc->Write();
 819   fhClustersDThetaAcc->Write();
 820   fhDPhiVsDThetaAcc->Write();
 821
 822   fhetaTracklets->Write();
 823   fhphiTracklets->Write();
 824   fhetaClustersLay1->Write();
 825   fhphiClustersLay1->Write();
 826 }
 827
 828 //____________________________________________________________________
 829 void
 830 AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
 831
 832   Float_t distClSameMod=0.;
 833   Float_t distClSameModMin=0.;
 834   Int_t   iClOverlap =0;
 835   Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
 836   Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
 837
 838   Float_t zproj1=0.;
 839   Float_t zproj2=0.;
 840   Float_t deZproj=0.;
 841   Float_t* clPar1  = GetClusterLayer1(iC1);
 842   Float_t* clPar2B = GetClusterLayer2(iC2WithBestDist);
 843   // Loop on inner layer clusters
 844   for (Int_t iiC1=0; iiC1<fNClustersLay1; iiC1++) {
 845     if (!fOverlapFlagClustersLay1[iiC1]) {
 846       // only for adjacent modules
 847       if ((TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==4)||
 848          (TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==76)) {
 849         Float_t *clPar11 = GetClusterLayer1(iiC1);
 850         Float_t dePhi=TMath::Abs(clPar11[kClPh]-clPar1[kClPh]);
 851         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 852
 853         zproj1=meanRadiusLay1/TMath::Tan(clPar1[kClTh]);
 854         zproj2=meanRadiusLay1/TMath::Tan(clPar11[kClTh]);
 855
 856         deZproj=TMath::Abs(zproj1-zproj2);
 857
 858         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 859         if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
 860
 861 //        if (distClSameMod<=1.) {
 862 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 863 //            distClSameModMin=distClSameMod;
 864 //            iClOverlap=iiC1;
 865 //          }
 866 //        }
 867
 868
 869       } // end adjacent modules
 870     }
 871   } // end Loop on inner layer clusters
 872
 873 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
 874
 875   distClSameMod=0.;
 876   distClSameModMin=0.;
 877   iClOverlap =0;
 878   // Loop on outer layer clusters
 879   for (Int_t iiC2=0; iiC2<fNClustersLay2; iiC2++) {
 880     if (!fOverlapFlagClustersLay2[iiC2]) {
 881       // only for adjacent modules
 882       Float_t *clPar2 = GetClusterLayer2(iiC2);
 883       if ((TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==4) ||
 884          (TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==156)) {
 885         Float_t dePhi=TMath::Abs(clPar2[kClPh]-clPar2B[kClPh]);
 886         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 887
 888         zproj1=meanRadiusLay2/TMath::Tan(clPar2B[kClTh]);
 889         zproj2=meanRadiusLay2/TMath::Tan(clPar2[kClTh]);
 890
 891         deZproj=TMath::Abs(zproj1-zproj2);
 892         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 893         if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE;
 894
 895 //        if (distClSameMod<=1.) {
 896 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 897 //            distClSameModMin=distClSameMod;
 898 //            iClOverlap=iiC2;
 899 //          }
 900 //        }
 901
 902       } // end adjacent modules
 903     }
 904   } // end Loop on outer layer clusters
 905
 906 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;
 907
 908 }
 909
 910 //____________________________________________________________________
 911 void AliITSMultReconstructor::ProcessESDTracks()
 912 {
 913   // Flag the clusters used by ESD tracks
 914   // Flag primary tracks to be used for multiplicity counting
 915   //
 916   if (!fESDEvent) return;
 917   AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexTracks();
 918   if (!vtx) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
 919   if (!vtx) {
 920     AliDebug(1,"No primary vertex: cannot flag primary tracks");
 921     return;
 922   }
 923   Int_t ntracks = fESDEvent->GetNumberOfTracks();
 924   for(Int_t itr=0; itr<ntracks; itr++) {
 925     AliESDtrack* track = fESDEvent->GetTrack(itr);
 926     if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx
 927     FlagTrackClusters(track);
 928     FlagIfPrimary(track,vtx);
 929   }
 930   //
 931 }
 932
 933 //____________________________________________________________________
 934 void AliITSMultReconstructor::FlagTrackClusters(const AliESDtrack* track)
 935 {
 936   // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation
 937   //
 938   Int_t idx[12];
 939   if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit
 940   //
 941   char mark = track->IsOn(AliESDtrack::kITSpureSA) ? kITSSAPBit : kITSTPCBit;
 942   char *uClus[2] = {fUsedClusLay1,fUsedClusLay2};
 943   for (int i=AliESDfriendTrack::kMaxITScluster;i--;) {
 944     // note: i>=6 is for extra clusters
 945     if (idx[i]<0) continue;
 946     int layID= (idx[i] & 0xf0000000) >> 28;
 947     if (layID>1) continue; // SPD only
 948     int clID = (idx[i] & 0x0fffffff);
 949     uClus[layID][clID] |= mark;
 950   }
 951   //
 952 }
 953
 954 //____________________________________________________________________
 955 void AliITSMultReconstructor::FlagIfPrimary(AliESDtrack* track, const AliVertex* vtx)
 956 {
 957   // RS: check if the track is primary and set the flag
 958   const double kPDCASPD1 = 0.1;
 959   const double kPDCASPD0 = 0.3;
 960   //
 961   double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? kPDCASPD1 : kPDCASPD0;
 962   // in principle, the track must already have been propagated to vertex
 963   /*
 964   Double_t dzRec[2]={0,0}, covdzRec[3];
 965   track->PropagateToDCA(vtx, fESDEvent->GetMagneticField(), 3.0, dzRec, covdzRec);
 966   */
 967   double dist = track->GetD(vtx->GetX(),vtx->GetY(),fESDEvent->GetMagneticField());
 968   if (TMath::Abs(dist*track->P())<cut) track->SetStatus(AliESDtrack::kMultPrimary);
 969 }