ITS/AliITSVertexerIons.cxx

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