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