ITS/AliITSVertexerPPZ.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 <AliITSVertexerPPZ.h>
  16 #include <TArrayF.h>
  17 #include <TH1F.h>
  18 #include <TObjArray.h>
  19 #include <TTree.h>
  20 #include <TClonesArray.h>
  21 #include "AliITSDetTypeRec.h"
  22 #include "AliITSgeom.h"
  23 #include "AliITSLoader.h"
  24 #include "AliITSRecPoint.h"
  25 /////////////////////////////////////////////////////////////////////////
  26 //                                                                     //
  27 // This class is intended to compute the Z coordinate of the           //
  28 // primary vertex for p-p interactions, or in general when the         //
  29 // number of secondaries is too slow to use the class                  //
  30 // AliITSVertexerIons, which in turn is optimized for A-A collsions    //
  31 // Origin: masera@to.infn.it      9/12/2002                            //
  32 //                                                                     //
  33 /////////////////////////////////////////////////////////////////////////
  34 ClassImp(AliITSVertexerPPZ)
  35
  36
  37
  38 //______________________________________________________________________
  39 AliITSVertexerPPZ::AliITSVertexerPPZ():AliITSVertexer() {
  40   // Default Constructor
  41
  42   SetDiffPhiMax(0);
  43   fX0 = 0.;
  44   fY0 = 0.;
  45   SetFirstLayerModules(0);
  46   SetSecondLayerModules(0);
  47   fZFound = 0;
  48   fZsig = 0.;
  49   //fITS = 0;
  50   SetWindow(0);
  51 }
  52
  53 AliITSVertexerPPZ::AliITSVertexerPPZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
  54   // Standard constructor
  55   SetDiffPhiMax();
  56   fX0 = x0;
  57   fY0 = y0;
  58   SetFirstLayerModules();
  59   SetSecondLayerModules();
  60   fZFound = 0;
  61   fZsig = 0.;
  62   //fITS = 0;
  63   SetWindow();
  64 }
  65
  66 //______________________________________________________________________
  67 AliITSVertexerPPZ::~AliITSVertexerPPZ() {
  68   // Default Destructor
  69   //fITS = 0;
  70 }
  71
  72 //________________________________________________________
  73 void  AliITSVertexerPPZ::EvalZ(TH1F *hist,Int_t sepa, Int_t ncoinc, TArrayF *zval) {
  74
  75   //Evaluation of Z
  76   Float_t deltaVal = hist->GetBinWidth(1)*fWindow; // max window in Z for searching
  77   fZFound=0;
  78   fZsig=0;
  79   Int_t nN=0;
  80   Int_t nbinNotZero=0;
  81   Float_t totst = 0.;
  82   Float_t totst2 = 0.;
  83   Float_t curz = 0.;
  84   for(Int_t i=1;i<=sepa;i++){
  85     Float_t cont=hist->GetBinContent(i);
  86     if(cont!=0)curz = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i);
  87     totst+=cont;
  88     totst2+=cont*cont;
  89     nN++;
  90     if(cont!=0)nbinNotZero++;
  91   }
  92   if(nbinNotZero==0){fZFound=-100; fZsig=-100; return;}
  93   if(nbinNotZero==1){
  94     fZFound = curz;
  95     fZsig=0;
  96     fCurrentVertex = new AliESDVertex(fZFound,fZsig,nbinNotZero);
  97     return;
  98   }
  99   Float_t errsq = totst2/(nN-1)-totst*totst/nN/(nN-1);
 100   if(errsq>=0){
 101   totst2=TMath::Sqrt(errsq);
 102   }
 103   else {
 104     Error("EvalZ","Negative variance: %d - %12.7f - %12.7f",nN,totst2,totst);
 105     fZFound=-100;
 106     fZsig=-100;
 107     return;
 108   }
 109   totst /= nN;
 110   Float_t cut = totst+totst2*2.;
 111   if(fDebug>1)cout<<"totst, totst2, cut: "<<totst<<", "<<totst2<<", "<<cut<<endl;
 112   Float_t val1=hist->GetBinLowEdge(sepa);
 113   Float_t val2=hist->GetBinLowEdge(1);
 114   Float_t valm = 0.;
 115   Float_t zmax = 0.;
 116   for(Int_t i=1;i<=sepa;i++){
 117     Float_t cont=hist->GetBinContent(i);
 118     if(cont>valm){
 119       valm = cont;
 120       zmax = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(1);
 121     }
 122     if(cont>cut){
 123       curz=hist->GetBinLowEdge(i);
 124       if(curz<val1)val1=curz;
 125       if(curz>val2)val2=curz;
 126     }
 127   }
 128   val2+=hist->GetBinWidth(1);
 129   if((val2-val1)>deltaVal){
 130     val1 = zmax-deltaVal/2.;
 131     val2 = zmax+deltaVal/2.;
 132     if(fDebug>0)cout<<"val1 and val2 recomputed\n";
 133   }
 134   if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<" "<<val2-val1<<endl;
 135   fZFound=0;
 136   fZsig=0;
 137   nN=0;
 138   for(Int_t i=0; i<ncoinc; i++){
 139     Float_t z=zval->At(i);
 140     if(z<val1)continue;
 141     if(z>val2)continue;
 142
 143     fZFound+=z;
 144     fZsig+=z*z;
 145     /*   weights defined by the curvature
 146     Float_t wei = 1./curv->At(i);
 147     fZFound+=z*wei;
 148     fZsig+=wei;
 149     */
 150     nN++;
 151   }
 152   if(nN<1){fZFound=-110; fZsig=-110; return;}
 153   if(nN==1){
 154     fZsig = 0;
 155     fCurrentVertex = new AliESDVertex(fZFound,fZsig,nN);
 156     return;
 157   }
 158   errsq = (fZsig/(nN-1)-fZFound*fZFound/nN/(nN-1))/nN;
 159   if(errsq>=0.){
 160     fZsig=TMath::Sqrt(errsq);
 161   }
 162   else {
 163     Error("evalZ","Negative variance: %d - %12.7f %12.7f",nN,fZsig,fZFound);
 164     fZsig=0;
 165   }
 166   fZFound=fZFound/nN;
 167   /* weights defined by the curvature
 168   fZsig=1./fZsig;
 169   fZFound*=fZsig;
 170   fZsig = TMath::Sqrt(fZsig);
 171   */
 172   fCurrentVertex = new AliESDVertex(fZFound,fZsig,nN);
 173 }
 174
 175 //______________________________________________________________________
 176 AliESDVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
 177   // Defines the AliESDVertex for the current event
 178   fCurrentVertex = 0;
 179   fZFound = -999;
 180   fZsig = -999;
 181   AliRunLoader *rl =AliRunLoader::GetRunLoader();
 182   AliITSLoader* iTSloader = (AliITSLoader*)rl->GetLoader("ITSLoader");
 183   TDirectory * olddir = gDirectory;
 184   rl->CdGAFile();
 185   AliITSgeom* geom = (AliITSgeom*)gDirectory->Get("AliITSgeom");
 186   olddir->cd();
 187
 188   AliITSDetTypeRec detTypeRec;
 189
 190   if(!geom) {
 191     Error("FindVertexForCurrentEvent","ITS geometry is not defined");
 192     return fCurrentVertex;
 193   }
 194   TTree *tR=0;
 195   TClonesArray *itsRec  = 0;
 196   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
 197   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
 198   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
 199   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
 200
 201   tR = iTSloader->TreeR();
 202   if(!tR){
 203     Error("FindVertexForCurrentEvent","TreeR not found");
 204     return fCurrentVertex;
 205   }
 206   detTypeRec.SetTreeAddressR(tR);
 207   itsRec = detTypeRec.RecPoints();
 208   // missing
 209   // TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
 210   TBranch *branch;
 211   //if(fUseV2Clusters){
 212   //  branch = tR->GetBranch("ITSRecPoints");
 213   //  branch->SetAddress(&clusters);
 214   //}
 215   //else {
 216   branch = tR->GetBranch("ITSRecPoints");
 217   if(!branch){
 218     branch = tR->GetBranch("ITSRecPointsF");
 219   }
 220   //}
 221   if(!branch){
 222    Error("FindVertexForCurrentEvent","branch for ITS rec points not found");
 223    return fCurrentVertex;
 224   }
 225   Float_t zave=0;
 226   Float_t rmszav=0;
 227   Float_t zave2=0;
 228   Int_t firipixe=0;
 229   for(Int_t module= fFirstL1; module<=fLastL1;module++){
 230     branch->GetEvent(module);
 231     //if(fUseV2Clusters){
 232     //  Clusters2RecPoints(clusters,module,itsRec);
 233     //}
 234     Int_t nrecp1 = itsRec->GetEntries();
 235     for(Int_t i=0; i<nrecp1;i++){
 236       AliITSRecPoint *current = (AliITSRecPoint*)itsRec->At(i);
 237       lc[0]=current->GetDetLocalX();
 238       lc[2]=current->GetDetLocalZ();
 239       geom->LtoG(module,lc,gc);
 240       zave+=gc[2];
 241       zave2+=gc[2]*gc[2];
 242       firipixe++;
 243     }
 244     detTypeRec.ResetRecPoints();
 245   }
 246   if(firipixe>1){
 247     rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
 248     zave=zave/firipixe;
 249     if(fDebug>1)cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++\n Number of firing pixels: "<<firipixe<<endl;
 250   }
 251   else {
 252     fZFound = -200;
 253     fZsig = -200;
 254     Warning("FindVertexForCurrentEvent","No rec points on first layer for this event");
 255     return fCurrentVertex;
 256   }
 257   Float_t zlim1=zave-rmszav;
 258   Float_t zlim2=zave+rmszav;
 259   Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
 260   zlim2=zlim1 + sepa/10.;
 261   TH1F *zvdis = new TH1F("z_ev","zv distr",sepa,zlim1,zlim2);
 262   if(fDebug>0){
 263     cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
 264     cout<<"Bin width: "<<zvdis->GetBinWidth(1)*10000.<<" micron\n";
 265   }
 266   Int_t sizarr=100;
 267   TArrayF *zval = new TArrayF(sizarr);
 268   //  TArrayF *curv = new TArrayF(sizarr);
 269   Int_t ncoinc=0;
 270   for(Int_t module= fFirstL1; module<=fLastL1;module++){
 271     if(fDebug>0)cout<<"processing module   "<<module<<"                  \r";
 272     branch->GetEvent(module);
 273     //if(fUseV2Clusters){
 274     //  Clusters2RecPoints(clusters,module,itsRec);
 275     //}
 276     Int_t nrecp1 = itsRec->GetEntries();
 277     TObjArray *poiL1 = new TObjArray(nrecp1);
 278     for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(itsRec->At(i),i);
 279     //fITS->ResetRecPoints();
 280     detTypeRec.ResetRecPoints();
 281     for(Int_t i=0; i<nrecp1;i++){
 282       AliITSRecPoint *current = (AliITSRecPoint*)poiL1->At(i);
 283       lc[0]=current->GetDetLocalX();
 284       lc[2]=current->GetDetLocalZ();
 285       geom->LtoG(module,lc,gc);
 286       gc[0]-=fX0;
 287       gc[1]-=fY0;
 288       Float_t r1=TMath::Sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
 289       Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
 290       if(phi1<0)phi1=2*TMath::Pi()+phi1;
 291       if(fDebug>1)cout<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1<<"     \n";
 292       for(Int_t modul2=fFirstL2; modul2<=fLastL2; modul2++){
 293         branch->GetEvent(modul2);
 294         //if(fUseV2Clusters){
 295         //  Clusters2RecPoints(clusters,modul2,itsRec);
 296         //}
 297         Int_t nrecp2 = itsRec->GetEntries();
 298         for(Int_t j=0; j<nrecp2;j++){
 299           AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
 300           lc2[0]=recp->GetDetLocalX();
 301           lc2[2]=recp->GetDetLocalZ();
 302           geom->LtoG(modul2,lc2,gc2);
 303           gc2[0]-=fX0;
 304           gc2[1]-=fY0;
 305           Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
 306           Float_t zr0=(r2*gc[2]-r1*gc2[2])/(r2-r1);
 307           Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
 308           if(phi2<0)phi2=2.*TMath::Pi()+phi2;
 309           Float_t diff = TMath::Abs(phi2-phi1);
 310           if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
 311           if(zr0>zlim1 && zr0<zlim2){
 312             if(diff<fDiffPhiMax ){
 313               zvdis->Fill(zr0);
 314               zval->AddAt(zr0,ncoinc);
 315               /* uncomment these lines to use curvature as a weight
 316               Float_t cu = Curv(0.,0.,gc[0],gc[1],gc2[0],gc[1]);
 317               curv->AddAt(cu,ncoinc);
 318               */
 319               ncoinc++;
 320               if(ncoinc==(sizarr-1)){
 321                 sizarr+=100;
 322                 zval->Set(sizarr);
 323                 //uncomment next line to use curvature as weight
 324                 //              curv->Set(sizarr);
 325               }
 326             }
 327           }
 328         }
 329         detTypeRec.ResetRecPoints();
 330       }
 331     }
 332     delete poiL1;
 333   }         // loop on modules
 334   if(fDebug>0){
 335     cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
 336   }
 337   //  EvalZ(zvdis,sepa,ncoinc,zval,curv);
 338   EvalZ(zvdis,sepa,ncoinc,zval);
 339   delete zvdis;
 340   delete zval;
 341   //  delete curv;
 342   if(fCurrentVertex){
 343     char name[30];
 344     sprintf(name,"Vertex_%d",evnumber);
 345     //    fCurrentVertex->SetName(name);
 346     fCurrentVertex->SetTitle("vertexer: PPZ");
 347   }
 348
 349   return fCurrentVertex;
 350 }
 351
 352 //______________________________________________________________________
 353 void AliITSVertexerPPZ::FindVertices(){
 354   // computes the vertices of the events in the range FirstEvent - LastEvent
 355   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 356   AliITSLoader* iTSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 357   iTSloader->ReloadRecPoints();
 358   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
 359     cout<<"Processing event "<<i<<endl;
 360     rl->GetEvent(i);
 361     FindVertexForCurrentEvent(i);
 362     if(fCurrentVertex){
 363       WriteCurrentVertex();
 364     }
 365     else {
 366       if(fDebug>0){
 367         cout<<"Vertex not found for event "<<i<<endl;
 368         cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
 369       }
 370     }
 371   }
 372 }
 373
 374 //________________________________________________________
 375 void AliITSVertexerPPZ::PrintStatus() const {
 376   // Print current status
 377   cout <<"=======================================================\n";
 378   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
 379   cout <<fLastL1<<endl;
 380   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
 381   cout <<fLastL2<<endl;
 382   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
 383   cout <<" Window for Z search: "<<fWindow<<endl;
 384   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
 385   cout <<" Debug flag: "<<fDebug<<endl;
 386   cout<<"First event to be processed "<<fFirstEvent;
 387   cout<<"\n Last event to be processed "<<fLastEvent<<endl;
 388 }
 389
 390 //________________________________________________________
 391 Float_t AliITSVertexerPPZ::Curv(Double_t x1,Double_t y1,
 392                    Double_t x2,Double_t y2,
 393                    Double_t x3,Double_t y3)
 394 {
 395   //-----------------------------------------------------------------
 396   // Initial approximation of the track curvature (Y. Belikov) squared
 397   //-----------------------------------------------------------------
 398   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 399   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 400                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 401   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 402                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 403
 404   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 405
 406   Float_t val = (Float_t)(( -xr*yr/TMath::Sqrt(xr*xr+yr*yr))
 407                           *( -xr*yr/TMath::Sqrt(xr*xr+yr*yr)));
 408   return val;
 409 }