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 "AliITSDetTypeRec.h"
  19 #include "AliITSVertexer3D.h"
  20 #include "AliITSVertexerZ.h"
  21 #include "AliESDVertex.h"
  22 #include "AliESDEvent.h"
  23 #include "AliRun.h"
  24 #include "AliRunLoader.h"
  25 #include "AliGenEventHeader.h"
  26 #include "AliStack.h"
  27 #endif
  28
  29 /*  $Id$    */
  30
  31 Bool_t DoVerticesSPD(Int_t pileupalgo=1, Int_t optdebug=0){
  32
  33   TFile *fint = new TFile("VertexSPD.root","recreate");
  34   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:isTriggered");
  35
  36   if (gClassTable->GetID("AliRun") < 0) {
  37     gInterpreter->ExecuteMacro("loadlibs.C");
  38   }
  39   // Set OCDB if needed
  40   AliCDBManager* man = AliCDBManager::Instance();
  41   if (!man->IsDefaultStorageSet()) {
  42     printf("Setting a local default storage and run number 0\n");
  43     man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  44     man->SetRun(0);
  45   }
  46   else {
  47     printf("Using deafult storage \n");
  48   }
  49
  50   // retrives geometry
  51   if(!gGeoManager){
  52     AliGeomManager::LoadGeometry("geometry.root");
  53   }
  54   AliGeomManager::ApplyAlignObjsFromCDB("ITS");
  55
  56   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  57   if (!runLoader) {
  58     Error("DoVertices", "getting run loader from file %s failed",
  59           "galice.root");
  60     return kFALSE;
  61   }
  62   runLoader->LoadgAlice();
  63   gAlice = runLoader->GetAliRun();
  64   if (!gAlice) {
  65     Error("DoVertices", "no galice object found");
  66     return kFALSE;
  67   }
  68   runLoader->LoadKinematics();
  69   runLoader->LoadHeader();
  70   AliITSLoader* ITSloader =  (AliITSLoader*) runLoader->GetLoader("ITSLoader");
  71   ITSloader->LoadRecPoints("read");
  72
  73   Int_t totev=runLoader->GetNumberOfEvents();
  74   if(optdebug)  printf("Number of events= %d\n",totev);
  75
  76   TFile* esdFile = TFile::Open("AliESDs.root");
  77   if (!esdFile || !esdFile->IsOpen()) {
  78     printf("Error in opening ESD file");
  79     return kFALSE;
  80   }
  81
  82   AliESDEvent * esd = new AliESDEvent;
  83   TTree* tree = (TTree*) esdFile->Get("esdTree");
  84   if (!tree) {
  85     printf("Error: no ESD tree found");
  86     return kFALSE;
  87   }
  88   esd->ReadFromTree(tree);
  89   tree->GetEvent(0);
  90   Int_t fieldkG=(Int_t)(TMath::Abs(esd->GetMagneticField())+0.001);
  91
  92   if(fieldkG==5){
  93     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 2, 1., 1., 10., AliMagF::k5kG));
  94   }else if(fieldkG==2){
  95     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 2, 1., 1., 10., AliMagF::k2kG));
  96   }else{
  97     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 2, 0., 1., 10., AliMagF::k5kG));
  98   }
  99   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 100   printf("Magnetic field set to %f\n",-fld->SolenoidField());
 101
 102   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
 103   Double_t xnom=0.,ynom=0.;
 104   AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
 105   vertz->Init("default");
 106   AliITSVertexer3D *vert3d = new AliITSVertexer3D();
 107   vert3d->Init("default");
 108   vert3d->SetWideFiducialRegion(40.,1.);
 109   vert3d->SetPileupAlgo(pileupalgo);
 110   vert3d->PrintStatus();
 111   vertz->SetDetTypeRec(detTypeRec);
 112   vert3d->SetDetTypeRec(detTypeRec);
 113
 114   /* uncomment these lines to use diamond constrain */
 115 //   Double_t posdiam[3]={0.03,0.1,0.};
 116 //   Double_t sigdiam[3]={0.01,0.01,10.0};
 117 //   AliESDVertex* diam=new AliESDVertex(posdiam,sigdiam);
 118 //   vertz->SetVtxStart(diam);
 119 //   vert3d->SetVtxStart(diam);
 120   /* end lines to be uncommented to use diamond constrain */
 121
 122   Int_t goodz=0,good3d=0;
 123   // Trigger mask
 124   ULong64_t spdFO = (1 << 14);
 125   ULong64_t v0left = (1 << 11);
 126   ULong64_t v0right = (1 << 12);
 127
 128   for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
 129     TArrayF mcVertex(3);
 130     runLoader->GetEvent(iEvent);
 131     runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 132     AliGenPythiaEventHeader *evh=(AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
 133     Int_t ptype = evh->ProcessType();
 134     if(optdebug){
 135       printf("==============================================================\n");
 136       printf("\nEvent: %d ---- Process Type = %d \n",iEvent,ptype);
 137     }
 138
 139     AliStack* stack = runLoader->Stack();
 140     TTree *treek=(TTree*)runLoader->TreeK();
 141     Int_t npart = (Int_t)treek->GetEntries();
 142     if(optdebug) printf("particles  %d\n",npart);
 143
 144     Double_t dNchdy = 0.;
 145
 146    // loop on particles to get generated dN/dy
 147     for(Int_t iPart=0; iPart<npart; iPart++) {
 148       if(!stack->IsPhysicalPrimary(iPart)) continue;
 149       TParticle* part = (TParticle*)stack->Particle(iPart);
 150       if(part->GetPDG()->Charge() == 0) continue;
 151       Double_t eta=part->Eta();
 152
 153       if(TMath::Abs(eta)<1.5) dNchdy+=1.;
 154     }
 155     if(optdebug) printf(" dNch/dy = %f\n",dNchdy);
 156
 157     tree->GetEvent(iEvent);
 158
 159     TTree* cltree = ITSloader->TreeR();
 160     ULong64_t triggerMask=esd->GetTriggerMask();
 161     // MB1: SPDFO || V0L || V0R
 162     Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)));
 163     AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
 164     AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
 165     AliMultiplicity *alimult = vert3d->GetMultiplicity();
 166     Int_t  ntrklets=0;
 167     Int_t nrecp1=0;
 168     if(alimult){
 169       ntrklets=alimult->GetNumberOfTracklets() ;
 170       nrecp1=ntrklets+alimult->GetNumberOfSingleClusters();
 171     }
 172
 173
 174
 175     TDirectory *current = gDirectory;
 176     fint->cd();
 177     Float_t xnt[23];
 178
 179     Double_t zz = 0.;
 180     Double_t zresz = 0.;
 181     Double_t zdiffz = 0.;
 182
 183     Int_t ntrkz = 0;
 184     if(vtxz){
 185
 186       ntrkz = vtxz->GetNContributors();
 187       if(ntrkz>0)goodz++;
 188       zz=vtxz->GetZv();
 189       zresz =vtxz->GetZRes(); // microns
 190       zdiffz = 10000.*(zz-mcVertex[2]); // microns
 191     }
 192
 193     Double_t x3d = 0.;
 194     Double_t xerr3d = 0.;
 195     Double_t xdiff3d = 0.;
 196
 197     Double_t y3d = 0.;
 198     Double_t yerr3d = 0.;
 199     Double_t ydiff3d = 0.;
 200
 201     Double_t z3d = 0.;
 202     Double_t zerr3d = 0.;
 203     Double_t zdiff3d = 0.;
 204
 205     Int_t ntrk3d = 0;
 206     Bool_t is3d=kFALSE;
 207     if(vtx3d){
 208
 209       if(vtx3d->IsFromVertexer3D()) is3d=kTRUE;
 210       ntrk3d = vtx3d->GetNContributors();
 211       if(is3d && ntrk3d>0)good3d++;
 212       x3d = vtx3d->GetXv();
 213       xerr3d=vtx3d->GetXRes();
 214       xdiff3d = 10000.*(x3d-mcVertex[0]);  // microns
 215
 216       y3d = vtx3d->GetYv();
 217       yerr3d=vtx3d->GetYRes();
 218       ydiff3d = 10000.*(y3d-mcVertex[1]);  // microns
 219
 220       z3d = vtx3d->GetZv();
 221       zerr3d=vtx3d->GetZRes();
 222       zdiff3d = 10000.*(z3d-mcVertex[2]);  // microns
 223
 224     }
 225
 226     xnt[0]=mcVertex[0];//x
 227     xnt[1]=mcVertex[1];//y
 228     xnt[2]=mcVertex[2];//z
 229
 230     xnt[3]=zz;
 231     xnt[4]=zdiffz;
 232     xnt[5]=zresz;
 233     xnt[6]=ntrkz;
 234
 235     xnt[7]=x3d;
 236     xnt[8]=xdiff3d;
 237     xnt[9]=xerr3d;
 238     xnt[10]=y3d;
 239     xnt[11]=ydiff3d;
 240     xnt[12]=yerr3d;
 241     xnt[13]=z3d;
 242     xnt[14]=zdiff3d;
 243     xnt[15]=zerr3d;
 244     xnt[16]=ntrk3d;
 245
 246     xnt[17]=dNchdy;
 247     xnt[18]=float(ntrklets);
 248     xnt[19]=float(nrecp1);
 249     xnt[20]=float(ptype);
 250     xnt[21]=float(is3d);
 251     xnt[22]=float(eventTriggered);
 252     nt->Fill(xnt);
 253     current->cd();
 254
 255     if(optdebug){
 256       printf("\nVertexerZ:  \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",zz,ntrkz,mcVertex[2],zdiffz);
 257       printf("Vertexer3D: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",z3d,ntrk3d,mcVertex[2],zdiff3d);
 258     }
 259
 260   }
 261   fint->cd();
 262   nt->Write();
 263   fint->Close();
 264   delete fint;
 265   if(optdebug){
 266    printf("***********************************************\n");
 267    printf("Number of Z vertices found= %d\n",goodz);
 268    printf("efficiency=%f\n",float(goodz)/float(totev));
 269    printf("Number of 3D vertices found= %d\n",good3d);
 270    printf("efficiency=%f\n",float(good3d)/float(totev));
 271   }
 272   return kTRUE;
 273 }