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 "AliITSgeom.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 the SPD inner
  29 // layer (SPD1). 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 SPD1 and define the fake vertex as
  33 // the mid-point of the straight line joining the two clusters.
  34 //   We reject the backgroung by requiring at least one clusters on SPD2
  35 // closer than fMaxDistOnSPD2 to the tracklet prolongation.
  36 //   We can reject (potentially pathological) events with the muon track
  37 // tangential to the SPD1 layer by the requiring the radial position of
  38 // the vertex to be smaller than fMaxVtxRadius.
  39 //   Due to background clusters, more than one vertex per event can
  40 // be found. We consider the first found.
  41 //   The number of contributors set in the AliESDVertex object is the
  42 // number of vertices found in the event; if this number is <1,
  43 // the procedure could not find a vertex position and by default
  44 // the vertex coordinates are set to (0,0,0) with large errors (100,100,100)
  45 // Number of contributors = 0  --> No SPD1 tracklets matching criteria
  46 // Number of contributors = -1 --> No SPD1 recpoints
  47 //
  48 // Origin: A.Dainese, andrea.dainese@lnl.infn.it
  49 //-------------------------------------------------------------------------
  50
  51 ClassImp(AliITSVertexerCosmics)
  52
  53 //-------------------------------------------------------------------------
  54 AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
  55 fFirstSPD1(0),
  56 fLastSPD1(0),
  57 fFirstSPD2(0),
  58 fLastSPD2(0),
  59 fMaxDistOnSPD2(0),
  60 fMaxVtxRadius(0),
  61 fMinDist2Vtxs(0)
  62 {
  63   // Default constructor
  64   SetSPD1Modules();
  65   SetSPD2Modules();
  66   SetMaxDistOnSPD2();
  67   SetMaxVtxRadius();
  68   SetMinDist2Vtxs();
  69 }
  70 //--------------------------------------------------------------------------
  71 AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(Int_t evnumber)
  72 {
  73   // Defines the AliESDVertex for the current event
  74
  75   fCurrentVertex = 0;
  76   AliRunLoader *rl =AliRunLoader::GetRunLoader();
  77   AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
  78   AliITSgeom* geom = itsLoader->GetITSgeom();
  79   itsLoader->LoadRecPoints();
  80   rl->GetEvent(evnumber);
  81
  82   TTree *rpTree = itsLoader->TreeR();
  83
  84   TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
  85   rpTree->SetBranchAddress("ITSRecPoints",&recpoints);
  86
  87   Double_t xclspd1[100],yclspd1[100],zclspd1[100],modclspd1[100];
  88   Int_t nclspd1stored=0;
  89   Double_t xclspd2[100],yclspd2[100],zclspd2[100],modclspd2[100];
  90   Int_t nclspd2stored=0;
  91   Int_t nrecpoints,nrecpointsSPD1=0;
  92
  93   Double_t gc[3]={0.,0.,0.};
  94   Double_t lc[3]={0.,0.,0.};
  95   Int_t lay,lad,det;
  96
  97   Double_t x[100],y[100],z[100],p1[3],p2[3],p3[3];
  98   Int_t nvtxs;
  99   Bool_t good,matchtospd2;
 100   Double_t xvtx,yvtx,zvtx,rvtx;
 101
 102   for(Int_t imodule=fFirstSPD1; imodule<fLastSPD2; imodule++) { // SPD
 103     rpTree->GetEvent(imodule);
 104     geom->GetModuleId(imodule,lay,lad,det);
 105     nrecpoints=recpoints->GetEntriesFast();
 106     if(imodule<fLastSPD1) nrecpointsSPD1 += nrecpoints;
 107     //printf("cosmics: module %d clusters %d\n",imodule,nrecpoints);
 108     for(Int_t irp=0; irp<nrecpoints; irp++) {
 109       AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
 110       // Local coordinates of this recpoint
 111       lc[0]=rp->GetDetLocalX();
 112       lc[2]=rp->GetDetLocalZ();
 113       geom->LtoG(imodule,lc,gc); // global coordinates
 114       if(lay==1) { // store SPD1 clusters
 115         xclspd1[nclspd1stored]=gc[0];
 116         yclspd1[nclspd1stored]=gc[1];
 117         zclspd1[nclspd1stored]=gc[2];
 118         modclspd1[nclspd1stored]=imodule;
 119         nclspd1stored++;
 120       }
 121       if(lay==2) { // store SPD2 clusters
 122         xclspd2[nclspd2stored]=gc[0];
 123         yclspd2[nclspd2stored]=gc[1];
 124         zclspd2[nclspd2stored]=gc[2];
 125         modclspd2[nclspd2stored]=imodule;
 126         nclspd2stored++;
 127       }
 128       if(nclspd1stored>100 || nclspd2stored>100)
 129         AliFatal("More than 100 clusters per layer in SPD");
 130     }// end clusters in a module
 131   }// end SPD modules for a given event
 132
 133   // build fake vertices
 134   nvtxs=0;
 135   // SPD1 - first cluster
 136   for(Int_t i1spd1=0; i1spd1<nclspd1stored; i1spd1++) {
 137     p1[0]=xclspd1[i1spd1]; p1[1]=yclspd1[i1spd1]; p1[2]=zclspd1[i1spd1];
 138     // SPD1 - second cluster
 139     for(Int_t i2spd1=i1spd1+1; i2spd1<nclspd1stored; i2spd1++) {
 140       if(modclspd1[i1spd1]==modclspd1[i2spd1]) continue;
 141       p2[0]=xclspd1[i2spd1]; p2[1]=yclspd1[i2spd1]; p2[2]=zclspd1[i2spd1];
 142       // look for point on SPD2
 143       AliStrLine spd1line(p1,p2,kTRUE);
 144       matchtospd2 = kFALSE;
 145       for(Int_t ispd2=0; ispd2<nclspd2stored; ispd2++) {
 146         p3[0]=xclspd2[ispd2]; p3[1]=yclspd2[ispd2]; p3[2]=zclspd2[ispd2];
 147         //printf(" %f\n",spd1line.GetDistFromPoint(p3));
 148         if(spd1line.GetDistFromPoint(p3)<fMaxDistOnSPD2)
 149           { matchtospd2 = kTRUE; break; }
 150       }
 151       if(!matchtospd2) continue;
 152       xvtx = 0.5*(xclspd1[i1spd1]+xclspd1[i2spd1]);
 153       yvtx = 0.5*(yclspd1[i1spd1]+yclspd1[i2spd1]);
 154       zvtx = 0.5*(zclspd1[i1spd1]+zclspd1[i2spd1]);
 155       rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
 156       if(rvtx>fMaxVtxRadius) continue;
 157       good = kTRUE;
 158       for(Int_t iv=0; iv<nvtxs; iv++) {
 159         if(TMath::Sqrt((xvtx- x[iv])*(xvtx- x[iv])+
 160                        (yvtx- y[iv])*(yvtx- y[iv])+
 161                        (zvtx- z[iv])*(zvtx- z[iv])) < fMinDist2Vtxs)
 162           good = kFALSE;
 163       }
 164       if(good) {
 165         x[nvtxs]=xvtx;
 166         y[nvtxs]=yvtx;
 167         z[nvtxs]=zvtx;
 168         nvtxs++;
 169       }
 170     } // SPD1 - second cluster
 171   } // SPD1 - first cluster
 172
 173
 174   Double_t pos[3]={0.,0.,0.};
 175   Double_t err[3]={100.,100.,100.};
 176   if(nvtxs) {
 177     pos[0]=x[0];
 178     pos[1]=y[0];
 179     pos[2]=z[0];
 180     err[0]=0.1;
 181     err[1]=0.1;
 182     err[2]=0.1;
 183   }
 184   if(!nrecpointsSPD1) nvtxs=-1;
 185   fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
 186   fCurrentVertex->SetNContributors(nvtxs);
 187   fCurrentVertex->SetTitle("cosmics fake vertex");
 188
 189   //if(nvtxs>0) fCurrentVertex->Print();
 190
 191   delete recpoints;
 192   itsLoader->UnloadRecPoints();
 193
 194   return fCurrentVertex;
 195 }
 196 //-------------------------------------------------------------------------
 197 void AliITSVertexerCosmics::FindVertices()
 198 {
 199   // computes the vertices of the events in the range FirstEvent - LastEvent
 200   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 201   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 202   itsLoader->ReloadRecPoints();
 203   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
 204     //  printf("Processing event %d\n",i);
 205     rl->GetEvent(i);
 206     FindVertexForCurrentEvent(i);
 207     if(fCurrentVertex){
 208       WriteCurrentVertex();
 209     }
 210   }
 211 }
 212 //-------------------------------------------------------------------------
 213 void AliITSVertexerCosmics::PrintStatus() const
 214 {
 215   // Print current status
 216   printf("=======================================================\n");
 217   printf(" fMaxDistOnSPD2: %f\n",fMaxDistOnSPD2);
 218   printf(" fMaxVtxRadius:  %f\n",fMaxVtxRadius);
 219   printf(" fMinDist2Vtxs:  %f\n",fMinDist2Vtxs);
 220   printf(" First layer first and last modules: %d, %d\n",fFirstSPD1,fLastSPD1);
 221   printf(" Second layer first and last modules: %d, %d\n",fFirstSPD2,fLastSPD2);
 222   printf("=======================================================\n");
 223 }
 224 //-------------------------------------------------------------------------