ITS/AliITSVertexerIons.cxx

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