ITS/AliITSVertexerIons.cxx

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