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   TFile *fint = new TFile("VertexSPD.root","recreate");
  31   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");
  32
  33   if (gClassTable->GetID("AliRun") < 0) {
  34     gInterpreter->ExecuteMacro("loadlibs.C");
  35   }
  36   // Set OCDB if needed
  37   AliCDBManager* man = AliCDBManager::Instance();
  38   if (!man->IsDefaultStorageSet()) {
  39     printf("Setting a local default storage and run number 0\n");
  40     man->SetDefaultStorage("local://$ALICE_ROOT");
  41     man->SetRun(0);
  42   }
  43   else {
  44     printf("Using deafult storage \n");
  45   }
  46
  47   // retrives geometry
  48   if(!gGeoManager){
  49     AliGeomManager::LoadGeometry("geometry.root");
  50   }
  51
  52   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  53   if (!runLoader) {
  54     Error("DoVertices", "getting run loader from file %s failed",
  55           "galice.root");
  56     return kFALSE;
  57   }
  58   runLoader->LoadgAlice();
  59   gAlice = runLoader->GetAliRun();
  60   if (!gAlice) {
  61     Error("DoVertices", "no galice object found");
  62     return kFALSE;
  63   }
  64   runLoader->LoadKinematics();
  65   runLoader->LoadHeader();
  66   AliITSLoader* ITSloader =  (AliITSLoader*) runLoader->GetLoader("ITSLoader");
  67   ITSloader->LoadRecPoints("read");
  68
  69   AliMagF * magf = gAlice->Field();
  70   AliTracker::SetFieldMap(magf,kTRUE);
  71   if(optdebug) printf("MagneticField=%f\n",AliTracker::GetBz());
  72
  73   Int_t totev=runLoader->GetNumberOfEvents();
  74   if(optdebug)  printf("Number of events= %d\n",totev);
  75
  76
  77 //   TFile* esdFile = TFile::Open("AliESDs.root");
  78 //   if (!esdFile || !esdFile->IsOpen()) {
  79 //     Error("DoVertices", "opening ESD file %s failed", "AliESDs.root");
  80 //     return kFALSE;
  81 //   }
  82 //   AliESD* esd = new AliESD;
  83 //   TTree* tree = (TTree*) esdFile->Get("esdTree");
  84 //   if (!tree) {
  85 //     Error("DoVertices", "no ESD tree found");
  86 //     return kFALSE;
  87 //   }
  88 //   tree->SetBranchAddress("ESD", &esd);
  89
  90   Double_t xnom=0.,ynom=0.;
  91   AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
  92   vertz->Init("default");
  93   AliITSVertexer3D *vert3d = new AliITSVertexer3D();
  94   vert3d->Init("default");
  95   //  vert3d->SetDebug(10);
  96   //  vertz->ConfigIterations(5);
  97
  98   Int_t goodz=0,good3d=0;
  99
 100   for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
 101     TArrayF mcVertex(3);
 102     runLoader->GetEvent(iEvent);
 103     runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 104     AliGenPythiaEventHeader *evh=(AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
 105     Int_t ptype = evh->ProcessType();
 106     if(optdebug){
 107       printf("==============================================================\n");
 108       printf("\nEvent: %d ---- Process Type = %d \n",iEvent,ptype);
 109     }
 110
 111     AliStack* stack = runLoader->Stack();
 112     TTree *treek=(TTree*)runLoader->TreeK();
 113     Int_t npart = (Int_t)treek->GetEntries();
 114     if(optdebug) printf("particles  %d\n",npart);
 115
 116     Double_t dNchdy = 0.;
 117
 118    // loop on particles
 119     for(Int_t pa=0; pa<npart; pa++) {
 120       TParticle* part = (TParticle*)stack->Particle(pa);
 121       Int_t pdg = part->GetPdgCode();
 122       Int_t apdg = TMath::Abs(pdg);
 123       Double_t energy  = part->Energy();
 124       if(energy>6900.) continue; // reject incoming protons
 125       Double_t pz = part->Pz();
 126       Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
 127
 128
 129       if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;      // reject secondaries
 130       if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
 131       if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
 132     }
 133     if(optdebug) printf(" dNch/dy = %f\n",dNchdy);
 134
 135     TTree* cltree = ITSloader->TreeR();
 136
 137     AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
 138     AliMultiplicity *alimult = vertz->GetMultiplicity();
 139     Int_t ntrklets=0,nrecp1=0;
 140     if(alimult) {
 141       nrecp1=alimult->GetNumberOfTracklets() ;
 142       ntrklets = alimult->GetNumberOfTracklets();
 143       for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
 144         if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
 145       }
 146     }
 147
 148     AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
 149
 150     TDirectory *current = gDirectory;
 151     fint->cd();
 152     Float_t xnt[21];
 153
 154     Double_t zz = 0.;
 155     Double_t zresz = 0.;
 156     Double_t zdiffz = 0.;
 157
 158     Int_t ntrkz = 0;
 159     if(vtxz){
 160
 161       ntrkz = vtxz->GetNContributors();
 162       if(ntrkz>0)goodz++;
 163       zz=vtxz->GetZv();
 164       zresz =vtxz->GetZRes(); // microns
 165       zdiffz = 10000.*(zz-mcVertex[2]); // microns
 166     }
 167
 168     Double_t x3d = 0.;
 169     Double_t xerr3d = 0.;
 170     Double_t xdiff3d = 0.;
 171
 172     Double_t y3d = 0.;
 173     Double_t yerr3d = 0.;
 174     Double_t ydiff3d = 0.;
 175
 176     Double_t z3d = 0.;
 177     Double_t zerr3d = 0.;
 178     Double_t zdiff3d = 0.;
 179
 180     Int_t ntrk3d = 0;
 181     if(vtx3d){
 182
 183       ntrk3d = vtx3d->GetNContributors();
 184       if(ntrk3d>0)good3d++;
 185       x3d = vtx3d->GetXv();
 186       xerr3d=vtx3d->GetXRes();
 187       xdiff3d = 10000.*(x3d-mcVertex[0]);  // microns
 188
 189       y3d = vtx3d->GetYv();
 190       yerr3d=vtx3d->GetYRes();
 191       ydiff3d = 10000.*(y3d-mcVertex[1]);  // microns
 192
 193       z3d = vtx3d->GetZv();
 194       zerr3d=vtx3d->GetZRes();
 195       zdiff3d = 10000.*(z3d-mcVertex[2]);  // microns
 196
 197     }
 198
 199     xnt[0]=mcVertex[0];//x
 200     xnt[1]=mcVertex[1];//y
 201     xnt[2]=mcVertex[2];//z
 202
 203     xnt[3]=zz;
 204     xnt[4]=zdiffz;
 205     xnt[5]=zresz;
 206     xnt[6]=ntrkz;
 207
 208     xnt[7]=x3d;
 209     xnt[8]=xdiff3d;
 210     xnt[9]=xerr3d;
 211     xnt[11]=y3d;
 212     xnt[11]=ydiff3d;
 213     xnt[12]=yerr3d;
 214     xnt[13]=z3d;
 215     xnt[14]=zdiff3d;
 216     xnt[15]=zerr3d;
 217     xnt[16]=ntrk3d;
 218
 219     xnt[17]=dNchdy;
 220     xnt[18]=float(ntrklets);
 221     xnt[19]=float(nrecp1);
 222     xnt[20]=float(ptype);
 223     nt->Fill(xnt);
 224     current->cd();
 225
 226     if(optdebug){
 227       printf("\nVertexerZ:  \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",zz,ntrkz,mcVertex[2],zdiffz);
 228       printf("Vertexer3D: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",z3d,ntrk3d,mcVertex[2],zdiff3d);
 229     }
 230
 231   }
 232   fint->cd();
 233   nt->Write();
 234   fint->Close();
 235   delete fint;
 236   if(optdebug){
 237    printf("***********************************************\n");
 238    printf("Number of Z vertices found= %d\n",goodz);
 239    printf("efficiency=%f\n",float(goodz)/float(totev));
 240    printf("Number of 3D vertices found= %d\n",good3d);
 241    printf("efficiency=%f\n",float(good3d)/float(totev));
 242   }
 243   return kTRUE;
 244 }