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