X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FDoVerticesSPD.C;h=187a63b1d7409923732794092774e020c49ae991;hb=e63dc6442e960bfb55149a971c344c7152e515e6;hp=0e92f384d77a3904762cd2726ea33c473fc98a13;hpb=5951029fa093300f404f85e5dbbc6b6583326b99;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/DoVerticesSPD.C b/ITS/DoVerticesSPD.C index 0e92f384d77..187a63b1d74 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, 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, 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, Int_t optdebug=0){ if(!gGeoManager){ AliGeomManager::LoadGeometry("geometry.root"); } + AliGeomManager::ApplyAlignObjsFromCDB("ITS"); AliRunLoader* runLoader = AliRunLoader::Open("galice.root"); if (!runLoader) { @@ -58,18 +62,21 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, 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); TFile* esdFile = TFile::Open("AliESDs.root"); @@ -85,31 +92,26 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, Int_t optdebug=0){ return kFALSE; } esd->ReadFromTree(tree); - tree->GetEvent(0); - Int_t fieldkG=(Int_t)(esd->GetMagneticField()+0.001); - - if(fieldkG==5){ - TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 2, 1., 1., 10., AliMagF::k5kG)); - }else if(fieldkG==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)); - } - AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); - printf("Magnetic field set to %f\n",-fld->SolenoidField()); + 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}; @@ -127,32 +129,33 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, 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(); @@ -168,7 +171,6 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, Int_t optdebug=0){ ntrklets=alimult->GetNumberOfTracklets() ; nrecp1=ntrklets+alimult->GetNumberOfSingleClusters(); } - TDirectory *current = gDirectory; @@ -242,7 +244,7 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, 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); @@ -252,8 +254,10 @@ Bool_t DoVerticesSPD(Int_t pileupalgo=1, 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); } }