ITS/AliITSVertexerIons.cxx

   1 #include "AliRunLoader.h"
   2 #include "AliITSDetTypeRec.h"
   3 #include "AliITSVertexerIons.h"
   4 #include "AliITSVertexerZ.h"
   5 #include "AliESDVertex.h"
   6 #include "AliITSgeomTGeo.h"
   7 #include "AliITSLoader.h"
   8 #include "AliITSRecPoint.h"
   9 #include "AliLog.h"
  10 #include <TTree.h>
  11 #include <TH1.h>
  12 #include <TF1.h>
  13
  14 //////////////////////////////////////////////////////////////////////
  15 // AliITSVertexerIons is a class for full 3D primary vertex         //
  16 // finding optimized for Ion-Ion interactions                       //
  17 //                                                                  //
  18 //                                                                  //
  19 //                                                                  //
  20 //                                                                  //
  21 // Written by Giuseppe Lo Re and Francesco Riggi                    //
  22 // Giuseppe.Lore@ct.infn.it                                         //
  23 // Franco.Riggi@ct.infn.it                                          //
  24 //                                                                  //
  25 // Release date: Mar 2004                                           //
  26 //                                                                  //
  27 //                                                                  //
  28 //////////////////////////////////////////////////////////////////////
  29
  30
  31 ClassImp(AliITSVertexerIons)
  32
  33 //______________________________________________________________________
  34   AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(),
  35 fNpThreshold(0),
  36 fMaxDeltaPhi(0),
  37 fMaxDeltaZ(0){
  38   // Default Constructor
  39
  40   //fITS = 0;
  41   SetNpThreshold();
  42   SetMaxDeltaPhi();
  43   SetMaxDeltaZ();
  44 }
  45
  46 //______________________________________________________________________
  47 AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn),
  48 fNpThreshold(0),
  49 fMaxDeltaPhi(0),
  50 fMaxDeltaZ(0) {
  51   // Standard constructor
  52
  53   //fITS = 0;
  54   SetNpThreshold();
  55   SetMaxDeltaPhi();
  56   SetMaxDeltaZ();
  57 }
  58 /*
  59 //______________________________________________________________________
  60 AliITSVertexerIons::AliITSVertexerIons(const AliITSVertexerIons &source):AliITSVertexer(source) {
  61   // Copy constructor
  62   // Copies are not allowed. The method is protected to avoid misuse.
  63   Error("AliITSVertexerIons","Copy constructor not allowed\n");
  64 }
  65 */
  66 //_________________________________________________________________________
  67 //AliITSVertexerIons& AliITSVertexerIons::operator=(const AliITSVertexerIons &/*source*/) {
  68   // Assignment operator
  69   // Assignment is not allowed. The method is protected to avoid misuse.
  70   //Error("= operator","Assignment operator not allowed\n");
  71   //return *this;
  72 //}
  73
  74
  75 //______________________________________________________________________
  76 AliITSVertexerIons::~AliITSVertexerIons() {
  77   // Default Destructor
  78   //fITS = 0;
  79 }
  80
  81 //______________________________________________________________________
  82 AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
  83 // Defines the AliESDVertex for the current event
  84
  85   fCurrentVertex = 0;
  86
  87   AliRunLoader *rl = AliRunLoader::GetRunLoader();
  88   AliITSLoader* itsloader = (AliITSLoader*)rl->GetLoader("ITSLoader");
  89   /*
  90   TDirectory * olddir = gDirectory;
  91   rl->CdGAFile();
  92   AliITSgeom* g2  = (AliITSgeom*)gDirectory->Get("AliITSgeom");
  93   olddir->cd();
  94   */
  95
  96   TTree *tr =  itsloader->TreeR();
  97   AliITSDetTypeRec detTypeRec;
  98
  99   detTypeRec.SetTreeAddressR(tr);
 100
 101   TClonesArray  *recpoints = detTypeRec.RecPoints();
 102   AliITSRecPoint *pnt;
 103
 104
 105   Int_t npoints=0;
 106   Int_t nopoints1=40000;
 107   Int_t nopoints2=40000;
 108   Float_t p[3];
 109
 110   Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
 111   z1=new Double_t[nopoints1];
 112   z2=new Double_t[nopoints2];
 113   y1=new Double_t[nopoints1];
 114   y2=new Double_t[nopoints2];
 115   x1=new Double_t[nopoints1];
 116   x2=new Double_t[nopoints2];
 117   phi1=new Double_t[nopoints1];
 118   phi2=new Double_t[nopoints2];
 119   r1=new Double_t[nopoints1];
 120   r2=new Double_t[nopoints2];
 121
 122   Double_t mxpiu = 0;
 123   Double_t mxmeno = 0;
 124   Double_t mypiu = 0;
 125   Double_t mymeno = 0;
 126   Double_t r=0;
 127
 128   Int_t np1=0, np2=0;
 129   for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) {
 130     detTypeRec.ResetRecPoints();
 131     tr->GetEvent(i);
 132     npoints = recpoints->GetEntries();
 133     for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
 134       pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
 135       /*
 136       l[0]=pnt->GetDetLocalX();
 137       l[1]=0;
 138       l[2]=pnt->GetDetLocalZ();
 139       g2->LtoG(i, l, p);
 140       */
 141       pnt->GetGlobalXYZ(p);
 142       r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
 143
 144       if(i<80 && TMath::Abs(p[2])<14.35)  {
 145                 y1[np1]=p[1];
 146                 x1[np1]=p[0];
 147                 z1[np1]=p[2];
 148                 if(p[0]>0) mxpiu++;
 149                 if(p[0]<0) mxmeno++;
 150                 if(p[1]>0) mypiu++;
 151                 if(p[1]<0) mymeno++;
 152                 np1++;
 153       }
 154       if(i>=80 && TMath::Abs(p[2])<14.35) {
 155                 y2[np2]=p[1];
 156                 x2[np2]=p[0];
 157                 z2[np2]=p[2];
 158                 np2++;
 159       }
 160     }
 161   }
 162
 163   if(np1<fNpThreshold) {
 164     Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerZ with default parameters...\n");
 165     Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
 166     AliITSVertexerZ *dovert = new AliITSVertexerZ("default");
 167     fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
 168     delete dovert;
 169     return fCurrentVertex;
 170   }
 171
 172   if(!np1 || !np2) {
 173     Error("FindVertexForCurrentEvent","No points in the pixer layers");
 174     return fCurrentVertex;
 175   }
 176   AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2));
 177
 178   Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
 179   Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
 180
 181   Double_t x0 = 2.86*asparx;
 182   Double_t y0 = 2.86*aspary;
 183
 184   AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0));
 185
 186   for(Int_t i=0;i<np1;i++) {
 187     x1[i]-=x0;
 188     y1[i]-=y0;
 189     PhiFunc(x1[i],y1[i],phi1[i]);
 190     r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
 191   }
 192   for(Int_t i=0;i<np2;i++) {
 193     x2[i]-=x0;
 194     y2[i]-=y0;
 195     PhiFunc(x2[i],y2[i],phi2[i]);
 196     r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
 197   }
 198
 199   Int_t nbinxy=400;
 200   Int_t nbinz=400;
 201   TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
 202   TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
 203   TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
 204
 205   Double_t dphi;
 206   for(Int_t j=0; j<np1; j++) {
 207     for(Int_t k=0; k<np2; k++) {
 208       dphi=TMath::Abs(phi2[k]-phi1[j]);
 209       if(dphi>180) dphi = 360-dphi;
 210       if(dphi>fMaxDeltaPhi) continue;
 211       hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
 212      }
 213   }
 214
 215   // Fitting ...
 216   Double_t max;
 217   Int_t binMax;
 218   Double_t maxcenter;
 219
 220   max = hzv->GetMaximum();
 221   binMax=hzv->GetMaximumBin();
 222   maxcenter=hzv->GetBinCenter(binMax);
 223   Double_t dxy=1.5;
 224
 225   TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy);
 226   fz->SetParameter(0,max);
 227   fz->SetParameter(1,maxcenter);
 228   fz->SetParameter(2,0.1);
 229   fz->SetLineColor(kRed);
 230   hzv->Fit("fz","RQ0");
 231
 232   if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) {
 233     Warning("FindVertexForCurrentEvent","AliITSVertexerIonsGeneral finder is not reliable for events with x and y vertex coordinates too far from the origin of the reference system. Their values will be set to 0.\n");
 234     Double_t position[3]={0,0,fz->GetParameter(1)};
 235     Double_t resolution[3]={0,0,0};
 236     Double_t snr[3]={0,0,0};
 237     Char_t name[30];
 238     AliDebug(1,Form("Vertex found for event %d",evnumber));
 239     sprintf(name,"Vertex");
 240     fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
 241     return fCurrentVertex;
 242   }
 243
 244   for(Int_t j=0; j<np1; j++) {
 245     for(Int_t k=0; k<np2; k++) {
 246       if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
 247       if(y2[k]==y1[j]) continue;
 248       hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
 249       if(x2[k]==x1[j]) continue;
 250       hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
 251     }
 252   }
 253
 254   TF1 *fx = new TF1 ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",x0-dxy,x0+dxy);
 255   fx->SetParameter(0,100);
 256   Double_t dist=0.3;
 257   Double_t xapprox=FindMaxAround(x0,hxv,dist);
 258   AliDebug(1,Form("xapprox = %f",xapprox));
 259   fx->SetParameter(1,xapprox);
 260   Double_t difcentroid=0.07;
 261   fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid);
 262   fx->SetParameter(2,0.1);
 263   fx->SetLineColor(kRed);
 264   hxv->Fit("fx","RQW0");
 265
 266   TF1 *fy = new TF1 ("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",y0-dxy,y0+dxy);
 267   fy->SetParameter(0,100);
 268   Double_t yapprox=FindMaxAround(y0,hyv,dist);
 269   AliDebug(1,Form("yapprox = %f",yapprox));
 270   fy->SetParameter(1,yapprox);
 271   fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid);
 272   fy->SetParameter(2,0.1);
 273   fy->SetLineColor(kRed);
 274   hyv->Fit("fy","RQW0");
 275
 276   delete [] z1;
 277   delete [] z2;
 278   delete [] y1;
 279   delete [] y2;
 280   delete [] x1;
 281   delete [] x2;
 282   delete [] r1;
 283   delete [] r2;
 284   delete [] phi1;
 285   delete [] phi2;
 286   delete hxv;
 287   delete hyv;
 288   delete hzv;
 289
 290   Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
 291   Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
 292   Double_t snr[3]={0,0,0};
 293
 294   Char_t name[30];
 295   AliDebug(1,Form("Vertex found for event %d",evnumber));
 296   sprintf(name,"Vertex");
 297   fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
 298
 299   return fCurrentVertex;
 300 }
 301
 302 //______________________________________________________________________
 303 void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
 304   // Method for the computation of the phi angle given x and y. The result is in degrees.
 305   if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
 306   if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
 307   if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
 308   if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;;
 309 }
 310
 311 //______________________________________________________________________
 312 void AliITSVertexerIons::FindVertices(){
 313   // computes the vertices of the events in the range FirstEvent - LastEvent
 314   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 315   AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 316   itsloader->LoadRecPoints("read");
 317   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
 318     rl->GetEvent(i);
 319     FindVertexForCurrentEvent(i);
 320     if(fCurrentVertex){
 321       WriteCurrentVertex();
 322     }
 323     else {
 324       AliDebug(1,Form("Vertex not found for event %d",i));
 325     }
 326   }
 327 }
 328
 329 //________________________________________________________
 330 void AliITSVertexerIons::PrintStatus() const {
 331   // Print current status
 332   cout <<"=======================================================\n";
 333   cout<<"First event to be processed "<<fFirstEvent;
 334   cout<<"\n Last event to be processed "<<fLastEvent<<endl;
 335   if(fCurrentVertex)fCurrentVertex->PrintStatus();
 336 }
 337
 338 //________________________________________________________
 339 Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
 340   // It finds the maximum of h within point-distance and distance+point.
 341   Int_t maxcontent=0;
 342   Int_t maxbin=0;
 343   for(Int_t i=0;i<h->GetNbinsX();i++) {
 344     Int_t content=(Int_t)h->GetBinContent(i);
 345     Double_t center=(Double_t)h->GetBinCenter(i);
 346     if(TMath::Abs(center-point)>distance) continue;
 347     if(content>maxcontent) {maxcontent=content;maxbin=i;}
 348   }
 349   Double_t max=h->GetBinCenter(maxbin);
 350   return max;
 351 }
 352
 353