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