ITS/DoVerticesSPD.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include "TFile.h"
   3 #include "TGeoManager.h"
   4 #include "TInterpreter.h"
   5 #include "TClassTable.h"
   6 #include "TNtuple.h"
   7 #include "TTree.h"
   8 #include "TArrayF.h"
   9 #include "TParticle.h"
  10 #include "AliESD.h"
  11 #include "AliGenPythiaEventHeader.h"
  12 #include "AliTracker.h"
  13 #include "AliHeader.h"
  14 #include "AliITSLoader.h"
  15 #include "AliVertexerTracks.h"
  16 #include "AliCDBManager.h"
  17 #include "AliGeomManager.h"
  18 #include "AliITSVertexer3D.h"
  19 #include "AliITSVertexerZ.h"
  20 #include "AliESDVertex.h"
  21 #include "AliRun.h"
  22 #include "AliRunLoader.h"
  23 #include "AliGenEventHeader.h"
  24 #include "AliStack.h"
  25 #endif
  26
  27 /*  $Id$    */
  28
  29 Bool_t DoVerticesSPD(Int_t optdebug=1){
  30
  31   TFile *fint = new TFile("VertexSPD.root","recreate");
  32   TNtuple *nt = new TNtuple("ntvert","vertices","xtrue:ytrue:ztrue:zZ:zdiffZ:zerrZ:ntrksZ:x3D:xdiff3D:xerr3D:y3D:ydiff3D:yerr3D:z3D:zdiff3D:zerr3D:ntrks3D:dndy:ntrklets:nrp1:ptyp:is3D");
  33
  34   if (gClassTable->GetID("AliRun") < 0) {
  35     gInterpreter->ExecuteMacro("loadlibs.C");
  36   }
  37   // Set OCDB if needed
  38   AliCDBManager* man = AliCDBManager::Instance();
  39   if (!man->IsDefaultStorageSet()) {
  40     printf("Setting a local default storage and run number 0\n");
  41     man->SetDefaultStorage("local://$ALICE_ROOT");
  42     man->SetRun(0);
  43   }
  44   else {
  45     printf("Using deafult storage \n");
  46   }
  47
  48   // retrives geometry
  49   if(!gGeoManager){
  50     AliGeomManager::LoadGeometry("geometry.root");
  51   }
  52
  53   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  54   if (!runLoader) {
  55     Error("DoVertices", "getting run loader from file %s failed",
  56           "galice.root");
  57     return kFALSE;
  58   }
  59   runLoader->LoadgAlice();
  60   gAlice = runLoader->GetAliRun();
  61   if (!gAlice) {
  62     Error("DoVertices", "no galice object found");
  63     return kFALSE;
  64   }
  65   runLoader->LoadKinematics();
  66   runLoader->LoadHeader();
  67   AliITSLoader* ITSloader =  (AliITSLoader*) runLoader->GetLoader("ITSLoader");
  68   ITSloader->LoadRecPoints("read");
  69
  70   AliMagF * magf = gAlice->Field();
  71   AliTracker::SetFieldMap(magf,kTRUE);
  72   if(optdebug) printf("MagneticField=%f\n",AliTracker::GetBz());
  73
  74   Int_t totev=runLoader->GetNumberOfEvents();
  75   if(optdebug)  printf("Number of events= %d\n",totev);
  76
  77
  78   Double_t xnom=0.,ynom=0.;
  79
  80   AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
  81   vertz->Init("default");
  82   AliITSVertexer3D *vert3d = new AliITSVertexer3D();
  83   vert3d->Init("default");
  84   vert3d->SetWideFiducialRegion(40.,1.);
  85   vert3d->PrintStatus();
  86
  87   /* uncomment these lines to use diamond constrain */
  88 //   Double_t posdiam[3]={0.03,0.1,0.};
  89 //   Double_t sigdiam[3]={0.01,0.01,10.0};
  90 //   AliESDVertex* diam=new AliESDVertex(posdiam,sigdiam);
  91 //   vertz->SetVtxStart(diam);
  92 //   vert3d->SetVtxStart(diam);
  93   /* end lines to be uncommented to use diamond constrain */
  94
  95   Int_t goodz=0,good3d=0;
  96
  97   for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
  98     TArrayF mcVertex(3);
  99     runLoader->GetEvent(iEvent);
 100     runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 101     AliGenPythiaEventHeader *evh=(AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
 102     Int_t ptype = evh->ProcessType();
 103     if(optdebug){
 104       printf("==============================================================\n");
 105       printf("\nEvent: %d ---- Process Type = %d \n",iEvent,ptype);
 106     }
 107
 108     AliStack* stack = runLoader->Stack();
 109     TTree *treek=(TTree*)runLoader->TreeK();
 110     Int_t npart = (Int_t)treek->GetEntries();
 111     if(optdebug) printf("particles  %d\n",npart);
 112
 113     Double_t dNchdy = 0.;
 114
 115    // loop on particles to get generated dN/dy
 116     for(Int_t iPart=0; iPart<npart; iPart++) {
 117       if(!stack->IsPhysicalPrimary(iPart)) continue;
 118       TParticle* part = (TParticle*)stack->Particle(iPart);
 119       if(part->GetPDG()->Charge() == 0) continue;
 120       Double_t eta=part->Eta();
 121
 122       if(TMath::Abs(eta)<1.5) dNchdy+=1.;
 123     }
 124     if(optdebug) printf(" dNch/dy = %f\n",dNchdy);
 125
 126     TTree* cltree = ITSloader->TreeR();
 127
 128     AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
 129     AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
 130     AliMultiplicity *alimult = vert3d->GetMultiplicity();
 131     Int_t  ntrklets=0;
 132     Int_t nrecp1=0;
 133     if(alimult){
 134       ntrklets=alimult->GetNumberOfTracklets() ;
 135       nrecp1=ntrklets+alimult->GetNumberOfSingleClusters();
 136     }
 137     TString tit=vtx3d->GetTitle();
 138     Bool_t is3d=kFALSE;
 139     if(tit.Contains("3D")) is3d=kTRUE;
 140
 141
 142
 143     TDirectory *current = gDirectory;
 144     fint->cd();
 145     Float_t xnt[22];
 146
 147     Double_t zz = 0.;
 148     Double_t zresz = 0.;
 149     Double_t zdiffz = 0.;
 150
 151     Int_t ntrkz = 0;
 152     if(vtxz){
 153
 154       ntrkz = vtxz->GetNContributors();
 155       if(ntrkz>0)goodz++;
 156       zz=vtxz->GetZv();
 157       zresz =vtxz->GetZRes(); // microns
 158       zdiffz = 10000.*(zz-mcVertex[2]); // microns
 159     }
 160
 161     Double_t x3d = 0.;
 162     Double_t xerr3d = 0.;
 163     Double_t xdiff3d = 0.;
 164
 165     Double_t y3d = 0.;
 166     Double_t yerr3d = 0.;
 167     Double_t ydiff3d = 0.;
 168
 169     Double_t z3d = 0.;
 170     Double_t zerr3d = 0.;
 171     Double_t zdiff3d = 0.;
 172
 173     Int_t ntrk3d = 0;
 174     if(vtx3d){
 175
 176       ntrk3d = vtx3d->GetNContributors();
 177       if(is3d && ntrk3d>0)good3d++;
 178       x3d = vtx3d->GetXv();
 179       xerr3d=vtx3d->GetXRes();
 180       xdiff3d = 10000.*(x3d-mcVertex[0]);  // microns
 181
 182       y3d = vtx3d->GetYv();
 183       yerr3d=vtx3d->GetYRes();
 184       ydiff3d = 10000.*(y3d-mcVertex[1]);  // microns
 185
 186       z3d = vtx3d->GetZv();
 187       zerr3d=vtx3d->GetZRes();
 188       zdiff3d = 10000.*(z3d-mcVertex[2]);  // microns
 189
 190     }
 191
 192     xnt[0]=mcVertex[0];//x
 193     xnt[1]=mcVertex[1];//y
 194     xnt[2]=mcVertex[2];//z
 195
 196     xnt[3]=zz;
 197     xnt[4]=zdiffz;
 198     xnt[5]=zresz;
 199     xnt[6]=ntrkz;
 200
 201     xnt[7]=x3d;
 202     xnt[8]=xdiff3d;
 203     xnt[9]=xerr3d;
 204     xnt[10]=y3d;
 205     xnt[11]=ydiff3d;
 206     xnt[12]=yerr3d;
 207     xnt[13]=z3d;
 208     xnt[14]=zdiff3d;
 209     xnt[15]=zerr3d;
 210     xnt[16]=ntrk3d;
 211
 212     xnt[17]=dNchdy;
 213     xnt[18]=float(ntrklets);
 214     xnt[19]=float(nrecp1);
 215     xnt[20]=float(ptype);
 216     xnt[21]=float(is3d);
 217     nt->Fill(xnt);
 218     current->cd();
 219
 220     if(optdebug){
 221       printf("\nVertexerZ:  \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",zz,ntrkz,mcVertex[2],zdiffz);
 222       printf("Vertexer3D: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",z3d,ntrk3d,mcVertex[2],zdiff3d);
 223     }
 224
 225   }
 226   fint->cd();
 227   nt->Write();
 228   fint->Close();
 229   delete fint;
 230   if(optdebug){
 231    printf("***********************************************\n");
 232    printf("Number of Z vertices found= %d\n",goodz);
 233    printf("efficiency=%f\n",float(goodz)/float(totev));
 234    printf("Number of 3D vertices found= %d\n",good3d);
 235    printf("efficiency=%f\n",float(good3d)/float(totev));
 236   }
 237   return kTRUE;
 238 }