HLT/ITS/AliHLTITSVertexerZ.cxx

   1 // $Id$
   2
   3 /**************************************************************************
   4  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   5  *                                                                        *
   6  * Author: The ALICE Off-line Project.                                    *
   7  * Contributors are mentioned in the code where appropriate.              *
   8  *                                                                        *
   9  * Permission to use, copy, modify and distribute this software and its   *
  10  * documentation strictly for non-commercial purposes is hereby granted   *
  11  * without fee, provided that the above copyright notice appears in all   *
  12  * copies and that both the copyright notice and this permission notice   *
  13  * appear in the supporting documentation. The authors make no claims     *
  14  * about the suitability of this software for any purpose. It is          *
  15  * provided "as is" without express or implied warranty.                  *
  16  **************************************************************************/
  17 #include "AliHLTITSVertexerZ.h"
  18 #include <TTree.h>
  19 #include<TString.h>
  20 #include<TH1.h>
  21 #include<TMath.h>
  22 #include <AliRun.h>
  23 #include <AliITS.h>
  24 #include "AliITSLoader.h"
  25 #include <AliITSgeom.h>
  26 #include <AliITSRecPoint.h>
  27 #include <AliITSclusterV2.h>
  28
  29 //-------------------------------------------------------------------------
  30 //                Implementation of the HLT ITS vertexer class
  31 //
  32 //          Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
  33 //-------------------------------------------------------------------------
  34
  35 ClassImp(AliHLTITSVertexerZ)
  36
  37 AliHLTITSVertexerZ::AliHLTITSVertexerZ():
  38   AliITSVertexerZ(),
  39   fZCombf(0),
  40   fStepFine(0)
  41 {
  42   // Constructor in case that there is no runloader
  43   SetBinWidthFine();
  44 }
  45
  46 AliHLTITSVertexerZ::AliHLTITSVertexerZ(Float_t x0, Float_t y0):
  47   AliITSVertexerZ(x0,y0),
  48   fZCombf(0),
  49   fStepFine(0)
  50 {
  51   // Standard Constructor
  52   SetBinWidthFine();
  53 }
  54
  55 AliHLTITSVertexerZ::~AliHLTITSVertexerZ()
  56 {
  57   // Destructor
  58   if (fZCombf) delete fZCombf;
  59 }
  60
  61 //______________________________________________________________________
  62 AliESDVertex* AliHLTITSVertexerZ::FindVertexForCurrentEvent(AliITSgeom *geom,TTree *tR){
  63   // Defines the AliESDVertex for the current event
  64
  65   fCurrentVertex = 0;
  66
  67   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
  68   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
  69   //Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
  70   //Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
  71
  72   TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
  73   TBranch *branch;
  74   branch = tR->GetBranch("Clusters");
  75   branch->SetAddress(&clusters);
  76
  77   Int_t nrpL1 = 0;
  78   Int_t nrpL2 = 0;
  79   for(Int_t module= fFirstL1; module<=fLastL1;module++){
  80     if(module%4==0 || module%4==3)continue;
  81     //   cout<<"Procesing module "<<module<<" ";
  82     tR->GetEvent(module);
  83     //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
  84     nrpL1+= clusters->GetEntriesFast();
  85   }
  86   for(Int_t module= fFirstL2; module<=fLastL2;module++){
  87     tR->GetEvent(module);
  88     nrpL2+= clusters->GetEntriesFast();
  89   }
  90   if(nrpL1 == 0 || nrpL2 == 0){
  91     ResetHistograms();
  92     return fCurrentVertex;
  93   }
  94
  95   Int_t nPhiBins = (Int_t)(TMath::Pi()/fDiffPhiMax);
  96   Float_t phiBinSize = 2*TMath::Pi()/(Float_t)nPhiBins;
  97
  98   Int_t maxind1 = 2*nrpL1/nPhiBins;
  99   Float_t **zc1 = new Float_t *[nPhiBins];
 100   Float_t **phi1 = new Float_t *[nPhiBins];
 101   Float_t **r1 = new Float_t *[nPhiBins];
 102   Int_t *ind1 = new Int_t [nPhiBins];
 103   Int_t maxind2 = 2*nrpL2/nPhiBins;
 104   Float_t **zc2 = new Float_t *[nPhiBins];
 105   Float_t **phi2 = new Float_t *[nPhiBins];
 106   Float_t **r2 = new Float_t *[nPhiBins];
 107   Int_t *ind2 = new Int_t [nPhiBins];
 108   for(Int_t i=0;i<nPhiBins;i++) {
 109     zc1[i] = new Float_t [maxind1];
 110     phi1[i] = new Float_t [maxind1];
 111     r1[i] = new Float_t [maxind1];
 112     zc2[i] = new Float_t [maxind2];
 113     phi2[i] = new Float_t [maxind2];
 114     r2[i] = new Float_t [maxind2];
 115   }
 116
 117   Float_t yshift = 0;
 118   Float_t zshift[4] = {-10.708000, -3.536000, 3.536000, 10.708000};
 119
 120   yshift = 0.248499;
 121   memset(ind1,0,nPhiBins*sizeof(Int_t));
 122   for(Int_t module= fFirstL1; module<=fLastL1;module++){
 123     if(module%4==0 || module%4==3)continue;
 124     tR->GetEvent(module);
 125     Int_t nrecp1 = clusters->GetEntriesFast();
 126     for(Int_t j=0;j<nrecp1;j++){
 127       AliITSclusterV2 *recp = (AliITSclusterV2*)clusters->UncheckedAt(j);
 128       lc[0]=-recp->GetY()+yshift;
 129       lc[2]=-recp->GetZ()+zshift[module%4];
 130       geom->LtoG(module,lc,gc);
 131       gc[0]-=GetNominalPos()[0];
 132       gc[1]-=GetNominalPos()[1];
 133       Float_t xc1,yc1;
 134       xc1=gc[0];
 135       yc1=gc[1];
 136       Float_t phi = TMath::ATan2(gc[1],gc[0]);
 137       if(phi<0)phi=2*TMath::Pi()+phi;
 138       Int_t bin = (Int_t)(phi/phiBinSize);
 139       if(bin>=nPhiBins || bin<0) bin = 0;
 140       Int_t ind = ind1[bin];
 141       if(ind<maxind1) {
 142         phi1[bin][ind] = phi;
 143         zc1[bin][ind]=gc[2]/fStepFine;
 144         r1[bin][ind]=sqrt(xc1*xc1+yc1*yc1);
 145         ind1[bin]++;
 146       }
 147     }
 148     clusters->Delete();
 149   }
 150   yshift = 3.096207;
 151   memset(ind2,0,nPhiBins*sizeof(Int_t));
 152   for(Int_t module= fFirstL2; module<=fLastL2;module++){
 153     tR->GetEvent(module);
 154     Int_t nrecp2 = clusters->GetEntriesFast();
 155     for(Int_t j=0;j<nrecp2;j++){
 156       AliITSclusterV2 *recp = (AliITSclusterV2*)clusters->UncheckedAt(j);
 157       lc[0]=recp->GetY()+yshift;
 158       lc[2]=-recp->GetZ()+zshift[module%4];
 159       geom->LtoG(module,lc,gc);
 160       gc[0]-=GetNominalPos()[0];
 161       gc[1]-=GetNominalPos()[1];
 162       Float_t xc2,yc2;
 163       xc2=gc[0];
 164       yc2=gc[1];
 165       Float_t phi = TMath::ATan2(gc[1],gc[0]);
 166       if(phi<0)phi=2*TMath::Pi()+phi;
 167       Int_t bin = (Int_t)(phi/phiBinSize+0.5);
 168       if(bin>=nPhiBins || bin<0) bin = 0;
 169       Int_t ind = ind2[bin];
 170       if(ind<maxind2) {
 171         phi2[bin][ind] = phi;
 172         zc2[bin][ind]=gc[2]/fStepFine;
 173         r2[bin][ind]=sqrt(xc2*xc2+yc2*yc2);
 174         ind2[bin]++;
 175       }
 176     }
 177     clusters->Delete();
 178   }
 179   Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
 180   Float_t lowz = fLowLim/fStepFine;
 181   Int_t *harray = new Int_t[nbinfine];
 182   memset(harray,0,nbinfine*sizeof(Int_t));
 183   for(Int_t ibin=0;ibin<nPhiBins;ibin++) {
 184     Float_t *pphi1 = phi1[ibin];
 185     Float_t *pr1 = r1[ibin];
 186     Float_t *pzc1 = zc1[ibin];
 187     for(Int_t i=0;i<ind1[ibin];i++){
 188       for(Int_t jbin=ibin;jbin<=(ibin+1);jbin++) {
 189         Int_t jbin2 = jbin;
 190         if(jbin==nPhiBins) jbin2 = 0;
 191         Float_t *pphi2 = phi2[jbin2];
 192         Float_t *pr2 = r2[jbin2];
 193         Float_t *pzc2 = zc2[jbin2];
 194         for(Int_t j=0;j<ind2[jbin2];j++){
 195           Float_t diff = TMath::Abs(pphi2[j]-pphi1[i]);
 196           if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
 197           if(diff<fDiffPhiMax){
 198             Float_t zr0=(pr2[j]*pzc1[i]-pr1[i]*pzc2[j])/(pr2[j]-pr1[i]);
 199             Int_t bin = (Int_t)(zr0-lowz);
 200             if(bin>=0 && bin<nbinfine){
 201               harray[bin]++;
 202             }
 203           }
 204         }
 205       }
 206     }
 207   }
 208   for(Int_t i=0;i<nPhiBins;i++) {
 209     delete [] zc1[i];
 210     delete [] phi1[i];
 211     delete [] r1[i];
 212     delete [] zc2[i];
 213     delete [] phi2[i];
 214     delete [] r2[i];
 215   }
 216   delete [] zc1;
 217   delete [] phi1;
 218   delete [] r1;
 219   delete [] ind1;
 220   delete [] zc2;
 221   delete [] phi2;
 222   delete [] r2;
 223   delete [] ind2;
 224
 225   Int_t nbinratio = (Int_t)(fStepCoarse/fStepFine+0.5);
 226   Int_t nbincoarse = nbinfine/nbinratio;
 227
 228   if(fZCombc)delete fZCombc;
 229   fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
 230   if(fZCombf)delete fZCombf;
 231   fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
 232   Stat_t contents=0;
 233   Int_t counter=0;
 234   for(Int_t bin=0;bin<nbinfine;bin++) {
 235     fZCombf->SetBinContent(bin+1,(Stat_t)harray[bin]);
 236     fZCombf->SetBinError(bin+1,TMath::Sqrt((Stat_t)harray[bin]));
 237     contents+=(Stat_t)harray[bin];
 238     counter++;
 239     if(counter==nbinratio) {
 240       Int_t binc=bin/nbinratio;
 241       fZCombc->SetBinContent(binc+1,contents);
 242       fZCombc->SetBinError(binc+1,TMath::Sqrt(contents));
 243       contents=0;
 244       counter=0;
 245     }
 246   }
 247
 248   delete [] harray;
 249
 250   if(fZCombc->GetEntries()==0){
 251     Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
 252     ResetHistograms();
 253     return fCurrentVertex;
 254   }
 255   //  else {
 256   //    cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
 257   //  }
 258   Int_t bi = fZCombc->GetMaximumBin();
 259   Float_t centre = fZCombc->GetBinCenter(bi);
 260   Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 261   Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 262   Int_t niter = 0;
 263   Bool_t goon = kTRUE;
 264   Int_t num;
 265   while(goon){
 266     fZFound = 0.;
 267     fZsig = 0.;
 268     num=0;
 269     for(Int_t n=n1;n<=n2;n++){
 270       fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
 271       num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
 272       fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
 273     }
 274     if(num<2){
 275       fZsig = 0.;
 276     }
 277     else {
 278       Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
 279       if(radi>0.)fZsig=TMath::Sqrt(radi);
 280       else fZsig=0.;
 281       fZFound/=num;
 282     }
 283     goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
 284     n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 285     n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
 286     niter++;
 287     if(niter>=10){
 288       goon = kFALSE;
 289       Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
 290     }
 291   }
 292   //  cout<<"Numer of Iterations "<<niter<<endl<<endl;
 293   fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
 294   fCurrentVertex->SetTitle("vertexer: HLT");
 295   ResetHistograms();
 296   PrintStatus();
 297   return fCurrentVertex;
 298 }