ITS/AliITSVertexerCosmics.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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 #include <TClonesArray.h>
  17 #include "AliLog.h"
  18 #include "AliESDVertex.h"
  19 #include "AliRunLoader.h"
  20 #include "AliITSLoader.h"
  21 #include "AliITSgeomTGeo.h"
  22 #include "AliITSRecPoint.h"
  23 #include "AliITSVertexerCosmics.h"
  24 #include "AliStrLine.h"
  25
  26 //------------------------------------------------------------------------
  27 // This class implements a method to construct a "fake" primary
  28 // vertex for cosmic events in which the muon crosses one of 5 inner
  29 // ITS layers. A fake primary vertex is needed for the reconstruction,
  30 // with e.g. AliITStrackerSA, of the two tracks produced by the muon
  31 // in the ITS.
  32 //   We build pairs of clusters on a given layer and define the fake vertex as
  33 // the mid-point of the straight line joining the two clusters.
  34 //   We use the innermost layer that has at least two clusters.
  35 //   We reject the background by requiring at least one cluster on the outer
  36 // layer, closer than fMaxDistOnOuterLayer to the tracklet prolongation.
  37 //   We can reject (potentially pathological) events with the muon track
  38 // tangential to the layer by the requiring the radial position of
  39 // the vertex to be smaller than fMaxVtxRadius.
  40 //   Due to background clusters, more than one vertex per event can
  41 // be found. We consider the first found.
  42 //   The errors on x,y,z of the vertex are calculated as errors on the mean
  43 // of clusters coordinates. Non-diag elements of vertex cov. mat. are set to 0.
  44 //   The number of contributors set in the AliESDVertex object is the
  45 // number of the layer on which the tracklet was built; if this number is -1,
  46 // the procedure could not find a vertex position and by default
  47 // the vertex coordinates are set to (0,0,0) with large errors (100,100,100)
  48 //
  49 // Origin: A.Dainese, andrea.dainese@lnl.infn.it
  50 //-------------------------------------------------------------------------
  51
  52 ClassImp(AliITSVertexerCosmics)
  53
  54 //-------------------------------------------------------------------------
  55 AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
  56 fMaxDistOnOuterLayer(0),
  57 fMinDist2Vtxs(0)
  58 {
  59   // Default constructor
  60   SetFirstLastModules(0,0,79);
  61   SetFirstLastModules(1,80,239);
  62   SetFirstLastModules(2,240,323);
  63   SetFirstLastModules(3,324,499);
  64   SetFirstLastModules(4,500,1247);
  65   SetFirstLastModules(5,1248,2197);
  66   SetMaxVtxRadius(0,3.5);
  67   SetMaxVtxRadius(1,6.5);
  68   SetMaxVtxRadius(2,14.5);
  69   SetMaxVtxRadius(3,23.5);
  70   SetMaxVtxRadius(4,37.5);
  71   SetMaxVtxRadius(5,42.5);
  72   SetMaxDistOnOuterLayer();
  73   SetMinDist2Vtxs();
  74 }
  75 //--------------------------------------------------------------------------
  76 AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(Int_t evnumber)
  77 {
  78   // Defines the AliESDVertex for the current event
  79
  80   fCurrentVertex = 0;
  81   AliRunLoader *rl =AliRunLoader::GetRunLoader();
  82   AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
  83   itsLoader->LoadRecPoints();
  84   rl->GetEvent(evnumber);
  85
  86   TTree *rpTree = itsLoader->TreeR();
  87
  88   TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
  89   rpTree->SetBranchAddress("ITSRecPoints",&recpoints);
  90
  91   Int_t lay,lad,det;
  92
  93   // Search for innermost layer with at least two clusters
  94   // on two different modules
  95   Int_t ilayer=0;
  96   while(ilayer<6) {
  97     Int_t nHitModules=0;
  98     for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer]; imodule++) {
  99       rpTree->GetEvent(imodule);
 100       AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
 101       lay -= 1;  // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
 102       if(lay!=ilayer) AliFatal("Layer mismatch!");
 103       if(recpoints->GetEntriesFast()>0) nHitModules++;
 104     }
 105     if(nHitModules>=2) break;
 106     ilayer++;
 107   }
 108   printf("Building tracklets on layer %d\n",ilayer);
 109
 110   const Int_t arrSize = 1000;
 111   Float_t xclInnLay[arrSize],yclInnLay[arrSize],zclInnLay[arrSize],modclInnLay[arrSize];
 112   Float_t e2xclInnLay[arrSize],e2yclInnLay[arrSize],e2zclInnLay[arrSize];
 113   Int_t nclInnLayStored=0;
 114   Float_t xclOutLay[arrSize],yclOutLay[arrSize],zclOutLay[arrSize],modclOutLay[arrSize];
 115   Int_t nclOutLayStored=0;
 116   Int_t nRecPoints,nRecPointsInnLay=0;
 117
 118   Float_t gc[3],gcov[5];
 119
 120   Float_t x[arrSize],y[arrSize],z[arrSize],e2x[arrSize],e2y[arrSize],e2z[arrSize];
 121   Double_t p1[3],p2[3],p3[3];
 122   Int_t nvtxs;
 123   Bool_t good,matchtoOutLay;
 124   Float_t xvtx,yvtx,zvtx,rvtx;
 125
 126   Double_t pos[3]={0.,0.,0.};
 127   Double_t err[3]={100.,100.,100.};
 128   Int_t ncontributors = -1;
 129
 130   // Collect clusters in the selected layer and the outer one
 131   for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer+1]; imodule++) {
 132     rpTree->GetEvent(imodule);
 133     AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
 134     lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
 135     nRecPoints=recpoints->GetEntriesFast();
 136     if(imodule<=fLast[ilayer]) nRecPointsInnLay += nRecPoints;
 137     //printf("cosmics: module %d clusters %d\n",imodule,nRecPoints);
 138     for(Int_t irp=0; irp<nRecPoints; irp++) {
 139       AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
 140       // Local coordinates of this recpoint
 141       rp->GetGlobalXYZ(gc);
 142       if(lay==ilayer) { // store InnLay clusters
 143         xclInnLay[nclInnLayStored]=gc[0];
 144         yclInnLay[nclInnLayStored]=gc[1];
 145         zclInnLay[nclInnLayStored]=gc[2];
 146         rp->GetGlobalCov(gcov);
 147         e2xclInnLay[nclInnLayStored]=gcov[0];
 148         e2yclInnLay[nclInnLayStored]=gcov[3];
 149         e2zclInnLay[nclInnLayStored]=gcov[5];
 150         modclInnLay[nclInnLayStored]=imodule;
 151         nclInnLayStored++;
 152       }
 153       if(lay==ilayer+1) { // store OutLay clusters
 154         xclOutLay[nclOutLayStored]=gc[0];
 155         yclOutLay[nclOutLayStored]=gc[1];
 156         zclOutLay[nclOutLayStored]=gc[2];
 157         modclOutLay[nclOutLayStored]=imodule;
 158         nclOutLayStored++;
 159       }
 160       if(nclInnLayStored>arrSize || nclOutLayStored>arrSize) {
 161         //AliFatal("More than arrSize clusters per layer");
 162         AliWarning("Too many clusters per layer");
 163         delete recpoints;
 164         itsLoader->UnloadRecPoints();
 165         fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
 166         fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
 167         fCurrentVertex->SetNContributors(ncontributors);
 168         return fCurrentVertex;
 169       }
 170     }// end clusters in a module
 171   }// end modules
 172
 173   // build fake vertices
 174   nvtxs=0;
 175   // InnLay - first cluster
 176   for(Int_t i1InnLay=0; i1InnLay<nclInnLayStored; i1InnLay++) {
 177     p1[0]=xclInnLay[i1InnLay];
 178     p1[1]=yclInnLay[i1InnLay];
 179     p1[2]=zclInnLay[i1InnLay];
 180     // InnLay - second cluster
 181     for(Int_t i2InnLay=i1InnLay+1; i2InnLay<nclInnLayStored; i2InnLay++) {
 182       if(modclInnLay[i1InnLay]==modclInnLay[i2InnLay]) continue;
 183       p2[0]=xclInnLay[i2InnLay];
 184       p2[1]=yclInnLay[i2InnLay];
 185       p2[2]=zclInnLay[i2InnLay];
 186       // look for point on OutLay
 187       AliStrLine InnLayline(p1,p2,kTRUE);
 188       matchtoOutLay = kFALSE;
 189       for(Int_t iOutLay=0; iOutLay<nclOutLayStored; iOutLay++) {
 190         p3[0]=xclOutLay[iOutLay];
 191         p3[1]=yclOutLay[iOutLay];
 192         p3[2]=zclOutLay[iOutLay];
 193         //printf(" %f\n",InnLayline.GetDistFromPoint(p3));
 194         if(InnLayline.GetDistFromPoint(p3)<fMaxDistOnOuterLayer)
 195           { matchtoOutLay = kTRUE; break; }
 196       }
 197       if(!matchtoOutLay) continue;
 198       xvtx = 0.5*(xclInnLay[i1InnLay]+xclInnLay[i2InnLay]);
 199       yvtx = 0.5*(yclInnLay[i1InnLay]+yclInnLay[i2InnLay]);
 200       zvtx = 0.5*(zclInnLay[i1InnLay]+zclInnLay[i2InnLay]);
 201       rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
 202       if(rvtx>fMaxVtxRadius[ilayer]) continue;
 203       good = kTRUE;
 204       for(Int_t iv=0; iv<nvtxs; iv++) {
 205         if(TMath::Sqrt((xvtx- x[iv])*(xvtx- x[iv])+
 206                        (yvtx- y[iv])*(yvtx- y[iv])+
 207                        (zvtx- z[iv])*(zvtx- z[iv])) < fMinDist2Vtxs)
 208           good = kFALSE;
 209       }
 210       if(good) {
 211         x[nvtxs]=xvtx;
 212         y[nvtxs]=yvtx;
 213         z[nvtxs]=zvtx;
 214         e2x[nvtxs]=0.25*(e2xclInnLay[i1InnLay]+e2xclInnLay[i2InnLay]);
 215         e2y[nvtxs]=0.25*(e2yclInnLay[i1InnLay]+e2yclInnLay[i2InnLay]);
 216         e2z[nvtxs]=0.25*(e2zclInnLay[i1InnLay]+e2zclInnLay[i2InnLay]);
 217         nvtxs++;
 218       }
 219     } // InnLay - second cluster
 220   } // InnLay - first cluster
 221
 222   if(nvtxs) {
 223     ncontributors = ilayer;
 224     pos[0]=x[0];
 225     pos[1]=y[0];
 226     pos[2]=z[0];
 227     err[0]=TMath::Sqrt(e2x[0]);
 228     err[1]=TMath::Sqrt(e2y[0]);
 229     err[2]=TMath::Sqrt(e2z[0]);
 230   }
 231
 232   fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
 233   fCurrentVertex->SetTitle("cosmics fake vertex");
 234   fCurrentVertex->SetNContributors(ncontributors);
 235   //fCurrentVertex->Print();
 236
 237   delete recpoints;
 238   itsLoader->UnloadRecPoints();
 239
 240   return fCurrentVertex;
 241 }
 242 //-------------------------------------------------------------------------
 243 void AliITSVertexerCosmics::FindVertices()
 244 {
 245   // computes the vertices of the events in the range FirstEvent - LastEvent
 246   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 247   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 248   itsLoader->ReloadRecPoints();
 249   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
 250     //  printf("Processing event %d\n",i);
 251     rl->GetEvent(i);
 252     FindVertexForCurrentEvent(i);
 253     if(fCurrentVertex){
 254       WriteCurrentVertex();
 255     }
 256   }
 257 }
 258 //-------------------------------------------------------------------------
 259 void AliITSVertexerCosmics::PrintStatus() const
 260 {
 261   // Print current status
 262   printf("=======================================================\n");
 263   printf(" fMaxDistOnOuterLayer: %f\n",fMaxDistOnOuterLayer);
 264   printf(" fMaxVtxRadius[0]:  %f\n",fMaxVtxRadius[0]);
 265   printf(" fMinDist2Vtxs:  %f\n",fMinDist2Vtxs);
 266   printf("=======================================================\n");
 267 }
 268 //-------------------------------------------------------------------------