ITS/AliITSVertexerIons.cxx

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