ITS/AliITSVertexerIons.cxx

   1 #include "AliITSDetTypeRec.h"
   2 #include "AliITSVertexerIons.h"
   3 #include "AliITSVertexerZ.h"
   4 #include "AliESDVertex.h"
   5 #include "AliITSgeomTGeo.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   /*
  89   TDirectory * olddir = gDirectory;
  90   rl->CdGAFile();
  91   AliITSgeom* g2  = (AliITSgeom*)gDirectory->Get("AliITSgeom");
  92   olddir->cd();
  93   */
  94
  95   TTree *tr =  itsloader->TreeR();
  96   AliITSDetTypeRec detTypeRec;
  97
  98   detTypeRec.SetTreeAddressR(tr);
  99
 100   TClonesArray  *recpoints = detTypeRec.RecPoints();
 101   AliITSRecPoint *pnt;
 102
 103
 104   Int_t npoints=0;
 105   Int_t nopoints1=40000;
 106   Int_t nopoints2=40000;
 107   Float_t p[3];
 108
 109   Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
 110   z1=new Double_t[nopoints1];
 111   z2=new Double_t[nopoints2];
 112   y1=new Double_t[nopoints1];
 113   y2=new Double_t[nopoints2];
 114   x1=new Double_t[nopoints1];
 115   x2=new Double_t[nopoints2];
 116   phi1=new Double_t[nopoints1];
 117   phi2=new Double_t[nopoints2];
 118   r1=new Double_t[nopoints1];
 119   r2=new Double_t[nopoints2];
 120
 121   Double_t mxpiu = 0;
 122   Double_t mxmeno = 0;
 123   Double_t mypiu = 0;
 124   Double_t mymeno = 0;
 125   Double_t r=0;
 126
 127   Int_t np1=0, np2=0;
 128   for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) {
 129     detTypeRec.ResetRecPoints();
 130     tr->GetEvent(i);
 131     npoints = recpoints->GetEntries();
 132     for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
 133       pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
 134       /*
 135       l[0]=pnt->GetDetLocalX();
 136       l[1]=0;
 137       l[2]=pnt->GetDetLocalZ();
 138       g2->LtoG(i, l, p);
 139       */
 140       pnt->GetGlobalXYZ(p);
 141       r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
 142
 143       if(i<80 && TMath::Abs(p[2])<14.35)  {
 144                 y1[np1]=p[1];
 145                 x1[np1]=p[0];
 146                 z1[np1]=p[2];
 147                 if(p[0]>0) mxpiu++;
 148                 if(p[0]<0) mxmeno++;
 149                 if(p[1]>0) mypiu++;
 150                 if(p[1]<0) mymeno++;
 151                 np1++;
 152       }
 153       if(i>=80 && TMath::Abs(p[2])<14.35) {
 154                 y2[np2]=p[1];
 155                 x2[np2]=p[0];
 156                 z2[np2]=p[2];
 157                 np2++;
 158       }
 159     }
 160   }
 161
 162   if(np1<fNpThreshold) {
 163     Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerZ with default parameters...\n");
 164     Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
 165     AliITSVertexerZ *dovert = new AliITSVertexerZ("default");
 166     fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
 167     delete dovert;
 168     return fCurrentVertex;
 169   }
 170
 171   if(!np1 || !np2) {
 172     Error("FindVertexForCurrentEvent","No points in the pixer layers");
 173     return fCurrentVertex;
 174   }
 175   AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2));
 176
 177   Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
 178   Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
 179
 180   Double_t x0 = 2.86*asparx;
 181   Double_t y0 = 2.86*aspary;
 182
 183   AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0));
 184
 185   for(Int_t i=0;i<np1;i++) {
 186     x1[i]-=x0;
 187     y1[i]-=y0;
 188     PhiFunc(x1[i],y1[i],phi1[i]);
 189     r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
 190   }
 191   for(Int_t i=0;i<np2;i++) {
 192     x2[i]-=x0;
 193     y2[i]-=y0;
 194     PhiFunc(x2[i],y2[i],phi2[i]);
 195     r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
 196   }
 197
 198   Int_t nbinxy=400;
 199   Int_t nbinz=400;
 200   TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
 201   TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
 202   TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
 203
 204   Double_t dphi;
 205   for(Int_t j=0; j<np1; j++) {
 206     for(Int_t k=0; k<np2; k++) {
 207       dphi=TMath::Abs(phi2[k]-phi1[j]);
 208       if(dphi>180) dphi = 360-dphi;
 209       if(dphi>fMaxDeltaPhi) continue;
 210       hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
 211      }
 212   }
 213
 214   // Fitting ...
 215   Double_t max;
 216   Int_t binMax;
 217   Double_t maxcenter;
 218
 219   max = hzv->GetMaximum();
 220   binMax=hzv->GetMaximumBin();
 221   maxcenter=hzv->GetBinCenter(binMax);
 222   Double_t dxy=1.5;
 223
 224   TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy);
 225   fz->SetParameter(0,max);
 226   fz->SetParameter(1,maxcenter);
 227   fz->SetParameter(2,0.1);
 228   fz->SetLineColor(kRed);
 229   hzv->Fit("fz","RQ0");
 230
 231   if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) {
 232     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");
 233     Double_t position[3]={0,0,fz->GetParameter(1)};
 234     Double_t resolution[3]={0,0,0};
 235     Double_t snr[3]={0,0,0};
 236     Char_t name[30];
 237     AliDebug(1,Form("Vertex found for event %d",evnumber));
 238     sprintf(name,"Vertex");
 239     fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
 240     return fCurrentVertex;
 241   }
 242
 243   for(Int_t j=0; j<np1; j++) {
 244     for(Int_t k=0; k<np2; k++) {
 245       if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
 246       if(y2[k]==y1[j]) continue;
 247       hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
 248       if(x2[k]==x1[j]) continue;
 249       hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
 250     }
 251   }
 252
 253   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);
 254   fx->SetParameter(0,100);
 255   Double_t dist=0.3;
 256   Double_t xapprox=FindMaxAround(x0,hxv,dist);
 257   AliDebug(1,Form("xapprox = %f",xapprox));
 258   fx->SetParameter(1,xapprox);
 259   Double_t difcentroid=0.07;
 260   fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid);
 261   fx->SetParameter(2,0.1);
 262   fx->SetLineColor(kRed);
 263   hxv->Fit("fx","RQW0");
 264
 265   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);
 266   fy->SetParameter(0,100);
 267   Double_t yapprox=FindMaxAround(y0,hyv,dist);
 268   AliDebug(1,Form("yapprox = %f",yapprox));
 269   fy->SetParameter(1,yapprox);
 270   fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid);
 271   fy->SetParameter(2,0.1);
 272   fy->SetLineColor(kRed);
 273   hyv->Fit("fy","RQW0");
 274
 275   delete [] z1;
 276   delete [] z2;
 277   delete [] y1;
 278   delete [] y2;
 279   delete [] x1;
 280   delete [] x2;
 281   delete [] r1;
 282   delete [] r2;
 283   delete [] phi1;
 284   delete [] phi2;
 285   delete hxv;
 286   delete hyv;
 287   delete hzv;
 288
 289   Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
 290   Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
 291   Double_t snr[3]={0,0,0};
 292
 293   Char_t name[30];
 294   AliDebug(1,Form("Vertex found for event %d",evnumber));
 295   sprintf(name,"Vertex");
 296   fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
 297
 298   return fCurrentVertex;
 299 }
 300
 301 //______________________________________________________________________
 302 void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
 303   // Method for the computation of the phi angle given x and y. The result is in degrees.
 304   if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
 305   if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
 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)+360;;
 308 }
 309
 310 //______________________________________________________________________
 311 void AliITSVertexerIons::FindVertices(){
 312   // computes the vertices of the events in the range FirstEvent - LastEvent
 313   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 314   AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 315   itsloader->LoadRecPoints("read");
 316   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
 317     rl->GetEvent(i);
 318     FindVertexForCurrentEvent(i);
 319     if(fCurrentVertex){
 320       WriteCurrentVertex();
 321     }
 322     else {
 323       AliDebug(1,Form("Vertex not found for event %d",i));
 324     }
 325   }
 326 }
 327
 328 //________________________________________________________
 329 void AliITSVertexerIons::PrintStatus() const {
 330   // Print current status
 331   cout <<"=======================================================\n";
 332   cout<<"First event to be processed "<<fFirstEvent;
 333   cout<<"\n Last event to be processed "<<fLastEvent<<endl;
 334   if(fCurrentVertex)fCurrentVertex->PrintStatus();
 335 }
 336
 337 //________________________________________________________
 338 Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
 339   // It finds the maximum of h within point-distance and distance+point.
 340   Int_t maxcontent=0;
 341   Int_t maxbin=0;
 342   for(Int_t i=0;i<h->GetNbinsX();i++) {
 343     Int_t content=(Int_t)h->GetBinContent(i);
 344     Double_t center=(Double_t)h->GetBinCenter(i);
 345     if(TMath::Abs(center-point)>distance) continue;
 346     if(content>maxcontent) {maxcontent=content;maxbin=i;}
 347   }
 348   Double_t max=h->GetBinCenter(maxbin);
 349   return max;
 350 }
 351
 352