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   // TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
 128   TBranch *branch;
 129   //if(fUseV2Clusters){
 130   //  branch = tR->GetBranch("ITSRecPoints");
 131   //  branch->SetAddress(&clusters);
 132   //}
 133   //else {
 134   branch = tR->GetBranch("ITSRecPoints");
 135     //}
 136
 137   Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
 138   Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
 139   if(fZCombc)delete fZCombc;
 140   fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
 141   if(fZCombf)delete fZCombf;
 142   fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
 143
 144   Int_t nrpL1 = 0;
 145   Int_t nrpL2 = 0;
 146   for(Int_t module= fFirstL1; module<=fLastL1;module++){
 147     if(module%4==0 || module%4==3)continue;
 148     //   cout<<"Procesing module "<<module<<" ";
 149     branch->GetEvent(module);
 150     //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
 151     //if(fUseV2Clusters){
 152     //  Clusters2RecPoints(clusters,module,itsRec);
 153     // }
 154     nrpL1+= itsRec->GetEntries();
 155     detTypeRec.ResetRecPoints();
 156   }
 157   for(Int_t module= fFirstL2; module<=fLastL2;module++){
 158     branch->GetEvent(module);
 159     //if(fUseV2Clusters){
 160     //  Clusters2RecPoints(clusters,module,itsRec);
 161     //}
 162     nrpL2+= itsRec->GetEntries();
 163     detTypeRec.ResetRecPoints();
 164   }
 165   if(nrpL1 == 0 || nrpL2 == 0){
 166     ResetHistograms();
 167     return fCurrentVertex;
 168   }
 169   Float_t *xc1 = new Float_t [nrpL1];
 170   Float_t *yc1 = new Float_t [nrpL1];
 171   Float_t *zc1 = new Float_t [nrpL1];
 172   Float_t *phi1 = new Float_t [nrpL1];
 173   Float_t *xc2 = new Float_t [nrpL2];
 174   Float_t *yc2 = new Float_t [nrpL2];
 175   Float_t *zc2 = new Float_t [nrpL2];
 176   Float_t *phi2 = new Float_t [nrpL2];
 177   Int_t ind = 0;
 178   for(Int_t module= fFirstL1; module<=fLastL1;module++){
 179     if(module%4==0 || module%4==3)continue;
 180     branch->GetEvent(module);
 181     //if(fUseV2Clusters){
 182     //  Clusters2RecPoints(clusters,module,itsRec);
 183     //}
 184     Int_t nrecp1 = itsRec->GetEntries();
 185     for(Int_t j=0;j<nrecp1;j++){
 186       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
 187       lc[0]=recp->GetDetLocalX();
 188       lc[2]=recp->GetDetLocalZ();
 189       geom->LtoG(module,lc,gc);
 190       gc[0]-=fX0;
 191       gc[1]-=fY0;
 192       xc1[ind]=gc[0];
 193       yc1[ind]=gc[1];
 194       zc1[ind]=gc[2];
 195       phi1[ind] = TMath::ATan2(gc[1],gc[0]);
 196       if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
 197       ind++;
 198     }
 199     detTypeRec.ResetRecPoints();
 200   }
 201   ind = 0;
 202   for(Int_t module= fFirstL2; module<=fLastL2;module++){
 203     branch->GetEvent(module);
 204     //if(fUseV2Clusters){
 205     //  Clusters2RecPoints(clusters,module,itsRec);
 206     //}
 207     Int_t nrecp2 = itsRec->GetEntries();
 208     for(Int_t j=0;j<nrecp2;j++){
 209       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
 210       lc[0]=recp->GetDetLocalX();
 211       lc[2]=recp->GetDetLocalZ();
 212       geom->LtoG(module,lc,gc);
 213       gc[0]-=fX0;
 214       gc[1]-=fY0;
 215       xc2[ind]=gc[0];
 216       yc2[ind]=gc[1];
 217       zc2[ind]=gc[2];
 218       phi2[ind] = TMath::ATan2(gc[1],gc[0]);
 219       if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
 220       ind++;
 221     }
 222     detTypeRec.ResetRecPoints();
 223   }
 224   for(Int_t i=0;i<nrpL1;i++){
 225     Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]);
 226     for(Int_t j=0;j<nrpL2;j++){
 227       Float_t diff = TMath::Abs(phi2[j]-phi1[i]);
 228       if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
 229       if(diff<fDiffPhiMax){
 230         Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]);
 231         Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1);
 232         fZCombf->Fill(zr0);
 233         fZCombc->Fill(zr0);
 234       }
 235     }
 236   }
 237   delete [] xc1;
 238   delete [] yc1;
 239   delete [] zc1;
 240   delete [] phi1;
 241   delete [] xc2;
 242   delete [] yc2;
 243   delete [] zc2;
 244   delete [] phi2;
 245   if(fZCombc->GetEntries()==0){
 246     Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
 247     ResetHistograms();
 248     return fCurrentVertex;
 249   }
 250   //  else {
 251   //    cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
 252   //  }
 253   Int_t bi = fZCombc->GetMaximumBin();
 254   Float_t centre = fZCombc->GetBinCenter(bi);
 255   Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 256   Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 257   Int_t niter = 0;
 258   Bool_t goon = kTRUE;
 259   Int_t num;
 260   while(goon){
 261     fZFound = 0.;
 262     fZsig = 0.;
 263     num=0;
 264     for(Int_t n=n1;n<=n2;n++){
 265       fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
 266       num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
 267       fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
 268     }
 269     if(num<2){
 270       fZsig = 0.;
 271     }
 272     else {
 273       Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
 274       if(radi>0.)fZsig=TMath::Sqrt(radi);
 275       else fZsig=0.;
 276       fZFound/=num;
 277     }
 278     goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
 279     n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 280     n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 281     niter++;
 282     if(niter>=10){
 283       goon = kFALSE;
 284       Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
 285     }
 286   }
 287   //  cout<<"Numer of Iterations "<<niter<<endl<<endl;
 288   fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
 289   fCurrentVertex->SetTitle("vertexer: B");
 290   ResetHistograms();
 291   return fCurrentVertex;
 292 }
 293
 294 //_____________________________________________________________________
 295 void AliITSVertexerZ::ResetHistograms(){
 296   // delete TH1 data members
 297   if(fZCombc)delete fZCombc;
 298   if(fZCombf)delete fZCombf;
 299   fZCombc = 0;
 300   fZCombf = 0;
 301 }
 302
 303 //______________________________________________________________________
 304 void AliITSVertexerZ::FindVertices(){
 305   // computes the vertices of the events in the range FirstEvent - LastEvent
 306   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 307   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 308   itsLoader->ReloadRecPoints();
 309   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
 310     cout<<"Processing event "<<i<<endl;
 311     rl->GetEvent(i);
 312     FindVertexForCurrentEvent(i);
 313     if(fCurrentVertex){
 314       WriteCurrentVertex();
 315     }
 316     else {
 317       if(fDebug>0){
 318         cout<<"Vertex not found for event "<<i<<endl;
 319         cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
 320       }
 321     }
 322   }
 323 }
 324
 325 //________________________________________________________
 326 void AliITSVertexerZ::PrintStatus() const {
 327   // Print current status
 328   cout <<"=======================================================\n";
 329   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
 330   cout <<fLastL1<<endl;
 331   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
 332   cout <<fLastL2<<endl;
 333   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
 334   cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
 335   cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
 336   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
 337   cout <<" Debug flag: "<<fDebug<<endl;
 338   cout <<"First event to be processed "<<fFirstEvent;
 339   cout <<"\n Last event to be processed "<<fLastEvent<<endl;
 340
 341
 342   cout <<"=======================================================\n";
 343 }
 344