ITS/AliITSVertexerZ.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2003, 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 #include <AliITSVertexerZ.h>
  16 #include <TString.h>
  17 #include "AliITSLoader.h"
  18 #include<TBranch.h>
  19 #include<TClonesArray.h>
  20 #include<TH1.h>
  21 #include<TTree.h>
  22 #include <AliITSgeom.h>
  23 #include "AliITSDetTypeRec.h"
  24 #include <AliITSRecPoint.h>
  25
  26 /////////////////////////////////////////////////////////////////
  27 // this class implements a fast method to determine
  28 // the Z coordinate of the primary vertex
  29 // for p-p collisions it seems to give comparable or better results
  30 // with respect to what obtained with AliITSVertexerPPZ
  31 // It can be used successfully with Pb-Pb collisions
  32 ////////////////////////////////////////////////////////////////
  33
  34 ClassImp(AliITSVertexerZ)
  35
  36
  37
  38 //______________________________________________________________________
  39 AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer() {
  40   // Default constructor
  41   SetDiffPhiMax(0);
  42   fX0 = 0.;
  43   fY0 = 0.;
  44   SetFirstLayerModules(0);
  45   SetSecondLayerModules(0);
  46   fZFound = 0;
  47   fZsig = 0.;
  48   fZCombc = 0;
  49   fZCombf = 0;
  50   SetLowLimit(0.);
  51   SetHighLimit(0.);
  52   SetBinWidthCoarse(0.);
  53   SetBinWidthFine(0.);
  54   SetTolerance(0.);
  55 }
  56
  57 //______________________________________________________________________
  58 AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
  59   // Standard Constructor
  60   SetDiffPhiMax();
  61   fX0 = x0;
  62   fY0 = y0;
  63   SetFirstLayerModules();
  64   SetSecondLayerModules();
  65   fZFound = 0;
  66   fZsig = 0.;
  67   fZCombc = 0;
  68   fZCombf = 0;
  69   SetLowLimit();
  70   SetHighLimit();
  71   SetBinWidthCoarse();
  72   SetBinWidthFine();
  73   SetTolerance();
  74
  75 }
  76
  77 //______________________________________________________________________
  78 AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr) {
  79   // Copy constructor
  80   // Copies are not allowed. The method is protected to avoid misuse.
  81   Error("AliITSVertexerZ","Copy constructor not allowed\n");
  82 }
  83
  84 //______________________________________________________________________
  85 AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ& /* vtxr */){
  86   // Assignment operator
  87   // Assignment is not allowed. The method is protected to avoid misuse.
  88   Error("= operator","Assignment operator not allowed\n");
  89   return *this;
  90 }
  91
  92
  93 //______________________________________________________________________
  94 AliITSVertexerZ::~AliITSVertexerZ() {
  95   // Default Destructor
  96   //fITS = 0;
  97   if(fZCombc)delete fZCombc;
  98   if(fZCombf)delete fZCombf;
  99 }
 100
 101 //______________________________________________________________________
 102 AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
 103   // Defines the AliESDVertex for the current event
 104
 105   fCurrentVertex = 0;
 106   AliRunLoader *rl =AliRunLoader::GetRunLoader();
 107   AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
 108   TDirectory * olddir = gDirectory;
 109   rl->CdGAFile();
 110   AliITSgeom* geom = (AliITSgeom*)gDirectory->Get("AliITSgeom");
 111   olddir->cd();
 112
 113   itsLoader->LoadRecPoints();
 114   rl->GetEvent(evnumber);
 115
 116   AliITSDetTypeRec detTypeRec;
 117
 118   TTree *tR = itsLoader->TreeR();
 119   detTypeRec.SetTreeAddressR(tR);
 120   TClonesArray *itsRec  = 0;
 121   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
 122   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
 123   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
 124   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
 125
 126   itsRec = detTypeRec.RecPoints();
 127   TBranch *branch;
 128   branch = tR->GetBranch("ITSRecPoints");
 129
 130   Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
 131   Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
 132   if(fZCombc)delete fZCombc;
 133   fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
 134   if(fZCombf)delete fZCombf;
 135   fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
 136
 137   Int_t nrpL1 = 0;
 138   Int_t nrpL2 = 0;
 139   for(Int_t module= fFirstL1; module<=fLastL1;module++){
 140     if(module%4==0 || module%4==3)continue;
 141     branch->GetEvent(module);
 142     nrpL1+= itsRec->GetEntries();
 143     detTypeRec.ResetRecPoints();
 144   }
 145   for(Int_t module= fFirstL2; module<=fLastL2;module++){
 146     branch->GetEvent(module);
 147     nrpL2+= itsRec->GetEntries();
 148     detTypeRec.ResetRecPoints();
 149   }
 150   if(nrpL1 == 0 || nrpL2 == 0){
 151     ResetHistograms();
 152     return fCurrentVertex;
 153   }
 154   Float_t *xc1 = new Float_t [nrpL1];
 155   Float_t *yc1 = new Float_t [nrpL1];
 156   Float_t *zc1 = new Float_t [nrpL1];
 157   Float_t *phi1 = new Float_t [nrpL1];
 158   Float_t *xc2 = new Float_t [nrpL2];
 159   Float_t *yc2 = new Float_t [nrpL2];
 160   Float_t *zc2 = new Float_t [nrpL2];
 161   Float_t *phi2 = new Float_t [nrpL2];
 162   Int_t ind = 0;
 163   for(Int_t module= fFirstL1; module<=fLastL1;module++){
 164     if(module%4==0 || module%4==3)continue;
 165     branch->GetEvent(module);
 166     Int_t nrecp1 = itsRec->GetEntries();
 167     for(Int_t j=0;j<nrecp1;j++){
 168       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
 169       lc[0]=recp->GetDetLocalX();
 170       lc[2]=recp->GetDetLocalZ();
 171       geom->LtoG(module,lc,gc);
 172       gc[0]-=fX0;
 173       gc[1]-=fY0;
 174       xc1[ind]=gc[0];
 175       yc1[ind]=gc[1];
 176       zc1[ind]=gc[2];
 177       phi1[ind] = TMath::ATan2(gc[1],gc[0]);
 178       if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
 179       ind++;
 180     }
 181     detTypeRec.ResetRecPoints();
 182   }
 183   ind = 0;
 184   for(Int_t module= fFirstL2; module<=fLastL2;module++){
 185     branch->GetEvent(module);
 186     Int_t nrecp2 = itsRec->GetEntries();
 187     for(Int_t j=0;j<nrecp2;j++){
 188       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
 189       lc[0]=recp->GetDetLocalX();
 190       lc[2]=recp->GetDetLocalZ();
 191       geom->LtoG(module,lc,gc);
 192       gc[0]-=fX0;
 193       gc[1]-=fY0;
 194       xc2[ind]=gc[0];
 195       yc2[ind]=gc[1];
 196       zc2[ind]=gc[2];
 197       phi2[ind] = TMath::ATan2(gc[1],gc[0]);
 198       if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
 199       ind++;
 200     }
 201     detTypeRec.ResetRecPoints();
 202   }
 203   for(Int_t i=0;i<nrpL1;i++){
 204     Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]);
 205     for(Int_t j=0;j<nrpL2;j++){
 206       Float_t diff = TMath::Abs(phi2[j]-phi1[i]);
 207       if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
 208       if(diff<fDiffPhiMax){
 209         Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]);
 210         Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1);
 211         fZCombf->Fill(zr0);
 212         fZCombc->Fill(zr0);
 213       }
 214     }
 215   }
 216   delete [] xc1;
 217   delete [] yc1;
 218   delete [] zc1;
 219   delete [] phi1;
 220   delete [] xc2;
 221   delete [] yc2;
 222   delete [] zc2;
 223   delete [] phi2;
 224   if(fZCombc->GetEntries()==0){
 225     Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
 226     ResetHistograms();
 227     return fCurrentVertex;
 228   }
 229   //  else {
 230   //    cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
 231   //  }
 232   Int_t bi = fZCombc->GetMaximumBin();
 233   Float_t centre = fZCombc->GetBinCenter(bi);
 234   Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 235   Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 236   Int_t niter = 0;
 237   Bool_t goon = kTRUE;
 238   Int_t num;
 239   while(goon){
 240     fZFound = 0.;
 241     fZsig = 0.;
 242     num=0;
 243     for(Int_t n=n1;n<=n2;n++){
 244       fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
 245       num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
 246       fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
 247     }
 248     if(num<2){
 249       fZsig = 0.;
 250     }
 251     else {
 252       Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
 253       if(radi>0.)fZsig=TMath::Sqrt(radi);
 254       else fZsig=0.;
 255       fZFound/=num;
 256     }
 257     goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
 258     n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 259     n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 260     niter++;
 261     if(niter>=10){
 262       goon = kFALSE;
 263       Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
 264     }
 265   }
 266   //  cout<<"Numer of Iterations "<<niter<<endl<<endl;
 267   fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
 268   fCurrentVertex->SetTitle("vertexer: B");
 269   ResetHistograms();
 270   return fCurrentVertex;
 271 }
 272
 273 //_____________________________________________________________________
 274 void AliITSVertexerZ::ResetHistograms(){
 275   // delete TH1 data members
 276   if(fZCombc)delete fZCombc;
 277   if(fZCombf)delete fZCombf;
 278   fZCombc = 0;
 279   fZCombf = 0;
 280 }
 281
 282 //______________________________________________________________________
 283 void AliITSVertexerZ::FindVertices(){
 284   // computes the vertices of the events in the range FirstEvent - LastEvent
 285   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 286   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 287   itsLoader->ReloadRecPoints();
 288   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
 289     cout<<"Processing event "<<i<<endl;
 290     rl->GetEvent(i);
 291     FindVertexForCurrentEvent(i);
 292     if(fCurrentVertex){
 293       WriteCurrentVertex();
 294     }
 295     else {
 296       if(fDebug>0){
 297         cout<<"Vertex not found for event "<<i<<endl;
 298         cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
 299       }
 300     }
 301   }
 302 }
 303
 304 //________________________________________________________
 305 void AliITSVertexerZ::PrintStatus() const {
 306   // Print current status
 307   cout <<"=======================================================\n";
 308   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
 309   cout <<fLastL1<<endl;
 310   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
 311   cout <<fLastL2<<endl;
 312   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
 313   cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
 314   cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
 315   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
 316   cout <<" Debug flag: "<<fDebug<<endl;
 317   cout <<"First event to be processed "<<fFirstEvent;
 318   cout <<"\n Last event to be processed "<<fLastEvent<<endl;
 319
 320
 321   cout <<"=======================================================\n";
 322 }
 323