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() {
  48   // Default Destructor
  49   //fITS = 0;
  50 }
  51
  52 //______________________________________________________________________
  53 AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(TTree *itsClusterTree){
  54 // Defines the AliESDVertex for the current event
  55
  56   fCurrentVertex = 0;
  57
  58   AliITSDetTypeRec detTypeRec;
  59
  60   detTypeRec.SetTreeAddressR(itsClusterTree);
  61
  62   TClonesArray  *recpoints = detTypeRec.RecPoints();
  63   AliITSRecPoint *pnt;
  64
  65
  66   Int_t npoints=0;
  67   Int_t nopoints1=40000;
  68   Int_t nopoints2=40000;
  69   Float_t p[3];
  70
  71   Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
  72   z1=new Double_t[nopoints1];
  73   z2=new Double_t[nopoints2];
  74   y1=new Double_t[nopoints1];
  75   y2=new Double_t[nopoints2];
  76   x1=new Double_t[nopoints1];
  77   x2=new Double_t[nopoints2];
  78   phi1=new Double_t[nopoints1];
  79   phi2=new Double_t[nopoints2];
  80   r1=new Double_t[nopoints1];
  81   r2=new Double_t[nopoints2];
  82
  83   Double_t mxpiu = 0;
  84   Double_t mxmeno = 0;
  85   Double_t mypiu = 0;
  86   Double_t mymeno = 0;
  87   Double_t r=0;
  88
  89   Int_t np1=0, np2=0;
  90   for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) {
  91     detTypeRec.ResetRecPoints();
  92     itsClusterTree->GetEvent(i);
  93     npoints = recpoints->GetEntries();
  94     for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
  95       pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
  96       /*
  97       l[0]=pnt->GetDetLocalX();
  98       l[1]=0;
  99       l[2]=pnt->GetDetLocalZ();
 100       g2->LtoG(i, l, p);
 101       */
 102       pnt->GetGlobalXYZ(p);
 103       r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
 104
 105       if(i<80 && TMath::Abs(p[2])<14.35)  {
 106                 y1[np1]=p[1];
 107                 x1[np1]=p[0];
 108                 z1[np1]=p[2];
 109                 if(p[0]>0) mxpiu++;
 110                 if(p[0]<0) mxmeno++;
 111                 if(p[1]>0) mypiu++;
 112                 if(p[1]<0) mymeno++;
 113                 np1++;
 114       }
 115       if(i>=80 && TMath::Abs(p[2])<14.35) {
 116                 y2[np2]=p[1];
 117                 x2[np2]=p[0];
 118                 z2[np2]=p[2];
 119                 np2++;
 120       }
 121     }
 122   }
 123
 124   if(np1<fNpThreshold) {
 125     Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerZ with default parameters...\n");
 126     Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
 127     AliITSVertexerZ *dovert = new AliITSVertexerZ();
 128     fCurrentVertex =dovert->FindVertexForCurrentEvent(itsClusterTree);
 129     delete dovert;
 130     return fCurrentVertex;
 131   }
 132
 133   if(!np1 || !np2) {
 134     Error("FindVertexForCurrentEvent","No points in the pixer layers");
 135     return fCurrentVertex;
 136   }
 137   AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2));
 138
 139   Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
 140   Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
 141
 142   Double_t x0 = 2.86*asparx;
 143   Double_t y0 = 2.86*aspary;
 144
 145   AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0));
 146
 147   for(Int_t i=0;i<np1;i++) {
 148     x1[i]-=x0;
 149     y1[i]-=y0;
 150     PhiFunc(x1[i],y1[i],phi1[i]);
 151     r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
 152   }
 153   for(Int_t i=0;i<np2;i++) {
 154     x2[i]-=x0;
 155     y2[i]-=y0;
 156     PhiFunc(x2[i],y2[i],phi2[i]);
 157     r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
 158   }
 159
 160   Int_t nbinxy=400;
 161   Int_t nbinz=400;
 162   TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
 163   TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
 164   TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
 165
 166   Double_t dphi;
 167   for(Int_t j=0; j<np1; j++) {
 168     for(Int_t k=0; k<np2; k++) {
 169       dphi=TMath::Abs(phi2[k]-phi1[j]);
 170       if(dphi>180) dphi = 360-dphi;
 171       if(dphi>fMaxDeltaPhi) continue;
 172       hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
 173      }
 174   }
 175
 176   // Fitting ...
 177   Double_t max;
 178   Int_t binMax;
 179   Double_t maxcenter;
 180
 181   max = hzv->GetMaximum();
 182   binMax=hzv->GetMaximumBin();
 183   maxcenter=hzv->GetBinCenter(binMax);
 184   Double_t dxy=1.5;
 185
 186   TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy);
 187   fz->SetParameter(0,max);
 188   fz->SetParameter(1,maxcenter);
 189   fz->SetParameter(2,0.1);
 190   fz->SetLineColor(kRed);
 191   hzv->Fit("fz","RQ0");
 192
 193   if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) {
 194     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");
 195     Double_t position[3]={0,0,fz->GetParameter(1)};
 196     Double_t resolution[3]={0,0,0};
 197     Double_t snr[3]={0,0,0};
 198     Char_t name[30];
 199     sprintf(name,"Vertex");
 200     fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
 201     return fCurrentVertex;
 202   }
 203
 204   for(Int_t j=0; j<np1; j++) {
 205     for(Int_t k=0; k<np2; k++) {
 206       if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
 207       if(y2[k]==y1[j]) continue;
 208       hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
 209       if(x2[k]==x1[j]) continue;
 210       hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
 211     }
 212   }
 213
 214   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);
 215   fx->SetParameter(0,100);
 216   Double_t dist=0.3;
 217   Double_t xapprox=FindMaxAround(x0,hxv,dist);
 218   AliDebug(1,Form("xapprox = %f",xapprox));
 219   fx->SetParameter(1,xapprox);
 220   Double_t difcentroid=0.07;
 221   fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid);
 222   fx->SetParameter(2,0.1);
 223   fx->SetLineColor(kRed);
 224   hxv->Fit("fx","RQW0");
 225
 226   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);
 227   fy->SetParameter(0,100);
 228   Double_t yapprox=FindMaxAround(y0,hyv,dist);
 229   AliDebug(1,Form("yapprox = %f",yapprox));
 230   fy->SetParameter(1,yapprox);
 231   fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid);
 232   fy->SetParameter(2,0.1);
 233   fy->SetLineColor(kRed);
 234   hyv->Fit("fy","RQW0");
 235
 236   delete [] z1;
 237   delete [] z2;
 238   delete [] y1;
 239   delete [] y2;
 240   delete [] x1;
 241   delete [] x2;
 242   delete [] r1;
 243   delete [] r2;
 244   delete [] phi1;
 245   delete [] phi2;
 246   delete hxv;
 247   delete hyv;
 248   delete hzv;
 249
 250   Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
 251   Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
 252   Double_t snr[3]={0,0,0};
 253
 254   Char_t name[30];
 255   sprintf(name,"Vertex");
 256   fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
 257
 258   return fCurrentVertex;
 259 }
 260
 261 //______________________________________________________________________
 262 void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
 263   // Method for the computation of the phi angle given x and y. The result is in degrees.
 264   if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
 265   if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
 266   if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
 267   if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;;
 268 }
 269
 270 //________________________________________________________
 271 void AliITSVertexerIons::PrintStatus() const {
 272   // Print current status
 273   cout <<"=======================================================\n";
 274   if(fCurrentVertex)fCurrentVertex->PrintStatus();
 275 }
 276
 277 //________________________________________________________
 278 Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
 279   // It finds the maximum of h within point-distance and distance+point.
 280   Int_t maxcontent=0;
 281   Int_t maxbin=0;
 282   for(Int_t i=0;i<h->GetNbinsX();i++) {
 283     Int_t content=(Int_t)h->GetBinContent(i);
 284     Double_t center=(Double_t)h->GetBinCenter(i);
 285     if(TMath::Abs(center-point)>distance) continue;
 286     if(content>maxcontent) {maxcontent=content;maxbin=i;}
 287   }
 288   Double_t max=h->GetBinCenter(maxbin);
 289   return max;
 290 }
 291
 292