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