X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FDoVerticesSPD.C;h=ffbe6d60de505059c13bdf62945116c4a964efaa;hb=dd028c53090c8bd4eb79119d6a7aae3928ed4410;hp=a4616468894ec6e53899f200134e9b8797fdbe47;hpb=7d766abb6347bab43de96c8aa697d034a859c8e8;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/DoVerticesSPD.C b/ITS/DoVerticesSPD.C index a4616468894..ffbe6d60de5 100644 --- a/ITS/DoVerticesSPD.C +++ b/ITS/DoVerticesSPD.C @@ -12,9 +12,11 @@ #include "AliTracker.h" #include "AliHeader.h" #include "AliITSLoader.h" +#include "AliITSsegmentationSPD.h" #include "AliVertexerTracks.h" #include "AliCDBManager.h" #include "AliGeomManager.h" +#include "AliGRPManager.h" #include "AliITSDetTypeRec.h" #include "AliITSVertexer3D.h" #include "AliITSVertexerZ.h" @@ -28,10 +30,10 @@ /* $Id$ */ -Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0){ +Bool_t DoVerticesSPD(Bool_t isMC=kFALSE, Int_t pileupalgo=1, Int_t optdebug=0){ TFile *fint = new TFile("VertexSPD.root","recreate"); - 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"); + TNtuple *nt = new TNtuple("ntvert","vertices","xtrue:ytrue:ztrue:zZ:zdiffZ:zerrZ:ntrksZ:x3D:xdiff3D:xerr3D:y3D:ydiff3D:yerr3D:z3D:zdiff3D:zerr3D:ntrks3D:nchgen:ntrklets:nrp1:ptyp:is3D:isTriggered"); if (gClassTable->GetID("AliRun") < 0) { gInterpreter->ExecuteMacro("loadlibs.C"); @@ -41,6 +43,7 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) if (!man->IsDefaultStorageSet()) { printf("Setting a local default storage and run number 0\n"); man->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); + // man->SetSpecificStorage("GRP/GRP/Data",Form("local://%s",gSystem->pwd())); man->SetRun(0); } else { @@ -51,6 +54,7 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) if(!gGeoManager){ AliGeomManager::LoadGeometry("geometry.root"); } + AliGeomManager::ApplyAlignObjsFromCDB("ITS"); AliRunLoader* runLoader = AliRunLoader::Open("galice.root"); if (!runLoader) { @@ -58,41 +62,56 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) "galice.root"); return kFALSE; } - runLoader->LoadgAlice(); - gAlice = runLoader->GetAliRun(); - if (!gAlice) { - Error("DoVertices", "no galice object found"); - return kFALSE; + + if(isMC){ + runLoader->LoadgAlice(); + gAlice = runLoader->GetAliRun(); + if (!gAlice) { + Error("DoVertices", "no galice object found"); + return kFALSE; + } + runLoader->LoadKinematics(); + runLoader->LoadHeader(); } - runLoader->LoadKinematics(); - runLoader->LoadHeader(); AliITSLoader* ITSloader = (AliITSLoader*) runLoader->GetLoader("ITSLoader"); ITSloader->LoadRecPoints("read"); - Int_t totev=runLoader->GetNumberOfEvents(); + Int_t totev=runLoader->GetNumberOfEvents(); if(optdebug) printf("Number of events= %d\n",totev); - AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec(); - if(bfield.Contains("5")){ - TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 2, 1., 1., 10., AliMagF::k5kG)); - }else if(bfield.Contains("2")){ - TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 2, 1., 1., 10., AliMagF::k2kG)); - }else{ - TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 2, 0., 1., 10., AliMagF::k5kG)); + TFile* esdFile = TFile::Open("AliESDs.root"); + if (!esdFile || !esdFile->IsOpen()) { + printf("Error in opening ESD file"); + return kFALSE; } - AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); - printf("Magnetic field set to %f\n",-fld->SolenoidField()); + + AliESDEvent * esd = new AliESDEvent; + TTree* tree = (TTree*) esdFile->Get("esdTree"); + if (!tree) { + printf("Error: no ESD tree found"); + return kFALSE; + } + esd->ReadFromTree(tree); + AliGRPManager * grpMan=new AliGRPManager(); + grpMan->ReadGRPEntry(); + grpMan->SetMagField(); + printf("Magnetic field set to %f T\n",AliTracker::GetBz()/10.); + + AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec(); + AliITSsegmentation* seg = new AliITSsegmentationSPD(); + detTypeRec->SetSegmentationModel(0,seg); + Double_t xnom=0.,ynom=0.; AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom); vertz->Init("default"); - AliITSVertexer3D *vert3d = new AliITSVertexer3D(); + AliITSVertexer3D *vert3d = new AliITSVertexer3D(); vert3d->Init("default"); vert3d->SetWideFiducialRegion(40.,1.); vert3d->SetPileupAlgo(pileupalgo); vert3d->PrintStatus(); vertz->SetDetTypeRec(detTypeRec); vert3d->SetDetTypeRec(detTypeRec); - + vert3d->SetComputeMultiplicity(kTRUE); /* uncomment these lines to use diamond constrain */ // Double_t posdiam[3]={0.03,0.1,0.}; // Double_t sigdiam[3]={0.01,0.01,10.0}; @@ -102,20 +121,6 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) /* end lines to be uncommented to use diamond constrain */ Int_t goodz=0,good3d=0; - TFile* esdFile = TFile::Open("AliESDs.root"); - if (!esdFile || !esdFile->IsOpen()) { - printf("Error in opening ESD file"); - return kFALSE; - } - - AliESDEvent * esd = new AliESDEvent; - TTree* tree = (TTree*) esdFile->Get("esdTree"); - if (!tree) { - printf("Error: no ESD tree found"); - return kFALSE; - } - esd->ReadFromTree(tree); - // Trigger mask ULong64_t spdFO = (1 << 14); ULong64_t v0left = (1 << 11); @@ -124,32 +129,33 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) for (Int_t iEvent = 0; iEvent < totev; iEvent++) { TArrayF mcVertex(3); runLoader->GetEvent(iEvent); - runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex); - AliGenPythiaEventHeader *evh=(AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader(); - Int_t ptype = evh->ProcessType(); + Double_t nChGen = 0.; + Int_t ptype = 0; if(optdebug){ printf("==============================================================\n"); - printf("\nEvent: %d ---- Process Type = %d \n",iEvent,ptype); + printf("Event: %d \n",iEvent); } - - AliStack* stack = runLoader->Stack(); - TTree *treek=(TTree*)runLoader->TreeK(); - Int_t npart = (Int_t)treek->GetEntries(); - if(optdebug) printf("particles %d\n",npart); - - Double_t dNchdy = 0.; - - // loop on particles to get generated dN/dy - for(Int_t iPart=0; iPartIsPhysicalPrimary(iPart)) continue; - TParticle* part = (TParticle*)stack->Particle(iPart); - if(part->GetPDG()->Charge() == 0) continue; - Double_t eta=part->Eta(); - - if(TMath::Abs(eta)<1.5) dNchdy+=1.; + if(isMC){ + runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex); + AliGenPythiaEventHeader *evh=(AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader(); + ptype = evh->ProcessType(); + + AliStack* stack = runLoader->Stack(); + TTree *treek=(TTree*)runLoader->TreeK(); + Int_t npart = (Int_t)treek->GetEntries(); + if(optdebug) printf("Process Type = %d --- Particles %d\n",ptype,npart); + + // loop on particles to get generated dN/dy + for(Int_t iPart=0; iPartIsPhysicalPrimary(iPart)) continue; + TParticle* part = (TParticle*)stack->Particle(iPart); + if(part->GetPDG()->Charge() == 0) continue; + Double_t eta=part->Eta(); + + if(TMath::Abs(eta)<1.5) nChGen+=1.; + } + if(optdebug) printf(" N. of generated charged primaries in |eta|<1.5 = %f\n",nChGen); } - if(optdebug) printf(" dNch/dy = %f\n",dNchdy); - tree->GetEvent(iEvent); TTree* cltree = ITSloader->TreeR(); @@ -165,7 +171,6 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) ntrklets=alimult->GetNumberOfTracklets() ; nrecp1=ntrklets+alimult->GetNumberOfSingleClusters(); } - TDirectory *current = gDirectory; @@ -181,7 +186,7 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) ntrkz = vtxz->GetNContributors(); if(ntrkz>0)goodz++; - zz=vtxz->GetZv(); + zz=vtxz->GetZ(); zresz =vtxz->GetZRes(); // microns zdiffz = 10000.*(zz-mcVertex[2]); // microns } @@ -205,15 +210,15 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) if(vtx3d->IsFromVertexer3D()) is3d=kTRUE; ntrk3d = vtx3d->GetNContributors(); if(is3d && ntrk3d>0)good3d++; - x3d = vtx3d->GetXv(); + x3d = vtx3d->GetX(); xerr3d=vtx3d->GetXRes(); xdiff3d = 10000.*(x3d-mcVertex[0]); // microns - y3d = vtx3d->GetYv(); + y3d = vtx3d->GetY(); yerr3d=vtx3d->GetYRes(); ydiff3d = 10000.*(y3d-mcVertex[1]); // microns - z3d = vtx3d->GetZv(); + z3d = vtx3d->GetZ(); zerr3d=vtx3d->GetZRes(); zdiff3d = 10000.*(z3d-mcVertex[2]); // microns @@ -239,7 +244,7 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) xnt[15]=zerr3d; xnt[16]=ntrk3d; - xnt[17]=dNchdy; + xnt[17]=nChGen; xnt[18]=float(ntrklets); xnt[19]=float(nrecp1); xnt[20]=float(ptype); @@ -249,8 +254,10 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, TString bfield="5kG", Int_t optdebug=0) current->cd(); if(optdebug){ - printf("\nVertexerZ: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",zz,ntrkz,mcVertex[2],zdiffz); - printf("Vertexer3D: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",z3d,ntrk3d,mcVertex[2],zdiff3d); + printf("Vertexer3D: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tContributors=%d \n",x3d,y3d,z3d,ntrk3d); + printf("VertexerZ: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tContributors=%d \n",0.,0.,zz,ntrkz); + if(isMC) printf("True Pos.: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tdN/dy=%.1f\n",mcVertex[0],mcVertex[1],mcVertex[2],nChGen); + printf("Multiplicity: Tracklets=%d ClustersLay1=%d\n",ntrklets,nrecp1); } }