From d4945a3be156626446bb3b80e18256f45caccca7 Mon Sep 17 00:00:00 2001 From: amastros Date: Thu, 15 Sep 2011 10:19:59 +0000 Subject: [PATCH] Fast Tool updates to compute efficiency with at least 3 clusters (SetAtLeastCorr(3)). new Config.C as a reference settings for upgrade studies FastVsSlowSim.C : macro to compare the two simulation methods (-> cdr plots) --- ITS/UPGRADE/Config.C | 442 ++++++++- ITS/UPGRADE/DetectorK.cxx | 402 ++++++-- ITS/UPGRADE/DetectorK.h | 329 ++++--- ITS/UPGRADE/macros/FastVsSlowSim.C | 1367 ++++++++++++++++++++++++++++ 4 files changed, 2309 insertions(+), 231 deletions(-) create mode 100644 ITS/UPGRADE/macros/FastVsSlowSim.C diff --git a/ITS/UPGRADE/Config.C b/ITS/UPGRADE/Config.C index a7533a25fee..7fba6efcad0 100644 --- a/ITS/UPGRADE/Config.C +++ b/ITS/UPGRADE/Config.C @@ -46,6 +46,12 @@ #include #endif + +Int_t generatorFlag = 2; +Int_t ITSflag = 3; + + + /* $Id$ */ enum PprTrigConf_t { @@ -57,6 +63,7 @@ const char * pprTrigConfName[] = { }; Float_t EtaToTheta(Float_t arg); +AliGenerator *Hijing(); static PprTrigConf_t strig = kDefaultPPTrig;// default PP trigger configuration @@ -76,6 +83,9 @@ void Config() gSystem->Load("libpythia6"); gSystem->Load("libAliPythia6"); gSystem->Load("libgeant321"); + gSystem->Load("libhijing"); + gSystem->Load("libTHijing"); + #endif new TGeant3TGeo("C++ Interface to Geant3"); @@ -151,25 +161,109 @@ void Config() // The cocktail itself - AliGenCocktail *gener = new AliGenCocktail(); - gener->SetPhiRange(0, 360); - // Set pseudorapidity range from -8 to 8. - Float_t thmin = EtaToTheta(8); // theta min. <---> eta max - Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min - gener->SetThetaRange(thmin,thmax); - gener->SetOrigin(0., 0., 0); //vertex position - gener->SetSigma(0., 0., 0); //Sigma in (X,Y,Z) (cm) on IP position - - AliGenBox *gbox2 = new AliGenBox(1); - gbox2->SetPtRange(0.99,1.0001); // - gbox2->SetPhiRange(0,360.); - gbox2->SetThetaRange(40.,140.); - gbox2->SetPart(kPiMinus); - gener->AddGenerator(gbox2,"GENBOX PIONS for ITS",1); - + if (generatorFlag==0) { + + AliGenCocktail *gener = new AliGenCocktail(); + gener->SetPhiRange(0, 360); + // Set pseudorapidity range from -8 to 8. + Float_t thmin = EtaToTheta(8); // theta min. <---> eta max + Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min + gener->SetThetaRange(thmin,thmax); + gener->SetOrigin(0., 0., 0); //vertex position + gener->SetSigma(0., 0., 0); //Sigma in (X,Y,Z) (cm) on IP position + /* + AliGenBox *gbox2 = new AliGenBox(2300*1.8); + gbox2->SetPtRange(0.1,2.0); // + gbox2->SetPhiRange(0,360); + gbox2->SetThetaRange(44.,135.); // +/- 0.9 eta + //gbox2->SetPhiRange(250,250.5); + // gbox2->SetThetaRange(50.,50.1); + gbox2->SetPart(kPiMinus); + gener->AddGenerator(gbox2,"GENBOX PIONS for ITS",1); + */ + + + + AliGenHIJINGpara *gbox2=new AliGenHIJINGpara(2300*1.8); // 2300 dNdy in 2*0.9eta + gbox2->SetPtRange(0.06,9999999); + gbox2->SetPhiRange(0,360.); + gbox2->SetThetaRange(44.,135.); // +/- 0.9 eta + gener->AddGenerator(gbox2,"GenHIJINGpara PIONS for ITS",1); - gener->Init(); + gener->Init(); + + + } else if (generatorFlag==1) { + + + AliGenHijing *generHijing = new AliGenHijing(-1); + // centre of mass energy + generHijing->SetEnergyCMS(5500.); + generHijing->SetImpactParameterRange(0,1); + // reference frame + generHijing->SetReferenceFrame("CMS"); + // projectile + generHijing->SetProjectile("A", 208, 82); + generHijing->SetTarget ("A", 208, 82); + // tell hijing to keep the full parent child chain + generHijing->KeepFullEvent(); + // enable jet quenching + generHijing->SetJetQuenching(1); + // enable shadowing + generHijing->SetShadowing(1); + // neutral pion and heavy particle decays switched off + // generHijing->SetDecaysOff(3); // 3 requested by Ana Marin + // Don't track spectators + generHijing->SetSpectators(0); + // kinematic selection + generHijing->SetSelectAll(0); + + AliGenerator* gener = generHijing; + gener->SetSigma(0, 0, 0); // Sigma in (X,Y,Z) (cm) on IP position + gener->SetVertexSmear(kPerEvent); + gener->Init(); + + + } else if (generatorFlag==2) { + + + // Francesco ... + int nParticles = 14022; + AliGenHIJINGpara *gener = new AliGenHIJINGpara(nParticles); + gener->SetMomentumRange(0.1, 10.); + gener->SetPhiRange(0., 360.); + Float_t thmin = EtaToTheta(2.5); // theta min. <---> eta max + Float_t thmax = EtaToTheta(-2.5); // theta max. <---> eta min + gener->SetThetaRange(thmin,thmax); + gener->SetOrigin(0, 0, 0); //vertex position + gener->SetSigma(0, 0, 0); //Sigma in (X,Y,Z) (cm) on IP position + gener->Init(); + + + } else if (generatorFlag==3) { + + // Annalisa ... + AliGenHijing *generHijing = new AliGenHijing(-1); + generHijing->SetEnergyCMS(5500.); + generHijing->SetImpactParameterRange(0,5); + generHijing->SetReferenceFrame("CMS"); + generHijing->SetProjectile("A", 208, 82); + generHijing->SetTarget ("A", 208, 82); + generHijing->KeepFullEvent(); + generHijing->SetJetQuenching(1); + generHijing->SetShadowing(1); + generHijing->SetSpectators(0); + generHijing->SetSelectAll(0); + generHijing->SetPtHardMin(4.5); + + AliGenerator* gener = generHijing; + gener->SetSigma(0, 0, 0); // Sigma in (X,Y,Z) (cm) on IP position + gener->SetVertexSmear(kPerEvent); + gener->Init(); + + + } // @@ -189,14 +283,15 @@ void Config() Int_t iFRAME = 0; Int_t iHALL = 0; Int_t iITS = 1; + Int_t isUpgrade = ITSflag; Int_t iMAG = 0; Int_t iMUON = 0; Int_t iPHOS = 0; - Int_t iPIPE = 1; + Int_t iPIPE = 0; Int_t iPMD = 0; - Int_t iHMPID = 0; + Int_t iHMPID = 0; Int_t iSHIL = 0; - Int_t iT0 = 0; + Int_t iT0 = 0; Int_t iTOF = 0; Int_t iTPC = 0; Int_t iTRD = 0; @@ -264,14 +359,302 @@ void Config() if (iITS) { //=================== ITS parameters ============================ - Bool_t isUpgrade = kTRUE; - AliITS *ITS = 0x0; - if(isUpgrade) ITS = new AliITSupgrade("ITS","ITS Upgrade"); - else ITS = new AliITSv11Hybrid("ITS","ITS v11Hybrid"); - - } + + if (isUpgrade==1) { // current ITS geometry ... + + + AliITSupgrade *ITSu = new AliITSupgrade("ITS","ITS upgrade",kTRUE); + + + } else if(isUpgrade==2) { // "New SPDs" Configuration + + // BeamPipe + Bool_t bp=kTRUE; + Double_t widthBP = 0.08; + Double_t radiusBP = 2.0-widthBP; + Double_t halflengthBP = 200.; + + + // BeamPipe + Int_t nlayers =7; + TArrayD xsize(nlayers); TArrayD zsize(nlayers); + TArrayD width(nlayers); TArrayD radii(nlayers); + + TArrayD widthCu(nlayers); TArrayD radiiCu(nlayers); + TArrayS copper(nlayers); + TArrayD halfLengths(nlayers); + + + for(Int_t i=0; iSetFullSegmentation(xsize,zsize); + + } else if (isUpgrade==3) { // "All New" Configuration + + // BeamPipe + Bool_t bp=kTRUE; + Double_t widthBP = 0.08; + Double_t radiusBP = 2.0-widthBP; + Double_t halflengthBP = 200.; + + + // BeamPipe + Int_t nlayers =7; + TArrayD xsize(nlayers); TArrayD zsize(nlayers); + TArrayD width(nlayers); TArrayD radii(nlayers); + + TArrayD widthCu(nlayers); TArrayD radiiCu(nlayers); + TArrayS copper(nlayers); + TArrayD halfLengths(nlayers); + + + for(Int_t i=0; iSetFullSegmentation(xsize,zsize); + + } else if (isUpgrade==4) { // "All New" Configuration with 8 layers + + // BeamPipe + Bool_t bp=kTRUE; + Double_t widthBP = 0.08; + Double_t radiusBP = 2.0-widthBP; + Double_t halflengthBP = 200.; + + + // BeamPipe + Int_t nlayers =8; + TArrayD xsize(nlayers); TArrayD zsize(nlayers); + TArrayD width(nlayers); TArrayD radii(nlayers); + + TArrayD widthCu(nlayers); TArrayD radiiCu(nlayers); + TArrayS copper(nlayers); + TArrayD halfLengths(nlayers); + + + for(Int_t i=0; iSetFullSegmentation(xsize,zsize); + + } else if (isUpgrade==5) { // "All New" Configuration with 8 layers + + // BeamPipe + Bool_t bp=kTRUE; + Double_t widthBP = 0.08; + Double_t radiusBP = 2.0-widthBP; + Double_t halflengthBP = 200.; + + + // BeamPipe + Int_t nlayers =9; + TArrayD xsize(nlayers); TArrayD zsize(nlayers); + TArrayD width(nlayers); TArrayD radii(nlayers); + + TArrayD widthCu(nlayers); TArrayD radiiCu(nlayers); + TArrayS copper(nlayers); + TArrayD halfLengths(nlayers); + + + for(Int_t i=0; iSetFullSegmentation(xsize,zsize); + + + } else { + + } else { + AliITS *ITS = new AliITSv11Hybrid("ITS","ITS v11Hybrid"); + } + } + + if (iTPC) { //============================ TPC parameters =================== @@ -376,3 +759,10 @@ void Config() Float_t EtaToTheta(Float_t arg){ return (180./TMath::Pi())*2.*atan(exp(-arg)); } + + +AliGenerator* Hijing() +{ + return gener; +} + diff --git a/ITS/UPGRADE/DetectorK.cxx b/ITS/UPGRADE/DetectorK.cxx index 75bf57722d1..0c671887f58 100644 --- a/ITS/UPGRADE/DetectorK.cxx +++ b/ITS/UPGRADE/DetectorK.cxx @@ -28,8 +28,8 @@ Changes by S. Rossegger -> see header file #define Luminosity 1.e27 // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 ) #define SigmaD 6.0 // Size of the interaction diamond (cm) (LHC = 6.0 cm) -#define dNdEtaMinB 950//660//950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700) -#define dNdEtaCent 2300 //1600//2300 // Multiplicity per unit Eta (LHC at 5.5 TeV not known) +#define dNdEtaMinB 1//950//660//950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700) +// #define dNdEtaCent 2300//15000 //1600//2300 // Multiplicity per unit Eta (LHC at 5.5 TeV not known) #define CrossSectionMinB 8 // minB Cross section for event under study (PbPb MinBias ~ 8 Barns) #define AcceptanceOfTpcAndSi 1 //1//0.60 //0.35 // Assumed geometric acceptance (efficiency) of the TPC and Si detectors @@ -42,6 +42,7 @@ Changes by S. Rossegger -> see header file #define KaonMass 0.498 // Mass of the Kaon #define D0Mass 1.865 // Mass of the D0 +//TMatrixD *probKomb; // table for efficiency kombinatorics class CylLayerK : public TNamed { @@ -92,13 +93,17 @@ DetectorK::DetectorK() : TNamed("test_detector","detector"), fNumberOfLayers(0), fNumberOfActiveLayers(0), + fNumberOfActiveITSLayers(0), fBField(0.5), fLhcUPCscale(1.0), fIntegrationTime(0.02), // in ms fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle fParticleMass(0.140), // Standard: pion mass - fMaxRadiusSlowDet(10.) + fMaxRadiusSlowDet(10.), + fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers + fAtLeastFake(1), // if at least x fakes, track is considered fake ... + fdNdEtaCent(2300) { // // default constructor @@ -111,13 +116,17 @@ DetectorK::DetectorK(char *name, char *title) : TNamed(name,title), fNumberOfLayers(0), fNumberOfActiveLayers(0), + fNumberOfActiveITSLayers(0), fBField(0.5), fLhcUPCscale(1.0), fIntegrationTime(0.02), // in ms fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle fParticleMass(0.140), // Standard: pion mass - fMaxRadiusSlowDet(10.) + fMaxRadiusSlowDet(10.), + fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers + fAtLeastFake(1), // if at least x fakes, track is considered fake ... + fdNdEtaCent(2300) { // // default constructor, that set the name and title @@ -168,7 +177,11 @@ void DetectorK::AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRe } fNumberOfLayers += 1; - if (!(newLayer->isDead)) fNumberOfActiveLayers += 1; + if (!(newLayer->isDead)) { + fNumberOfActiveLayers += 1; + TString lname(newLayer->GetName()); + if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1; + } } else { @@ -192,6 +205,8 @@ void DetectorK::KillLayer(char *name) { if (!(tmp->isDead)) { tmp->isDead = kTRUE; fNumberOfActiveLayers -= 1; + TString lname(tmp->GetName()); + if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1; } } } @@ -273,13 +288,20 @@ void DetectorK::SetResolution(char *name, Float_t phiRes, Float_t zRes) { tmp->phiRes = phiRes; tmp->zRes = zRes; - + TString lname(tmp->GetName()); + if (zRes==RIDICULOUS && phiRes==RIDICULOUS) { tmp->isDead = kTRUE; - if (!wasDead) fNumberOfActiveLayers -= 1; + if (!wasDead) { + fNumberOfActiveLayers -= 1; + if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1; + } } else { tmp->isDead = kFALSE; - if (wasDead) fNumberOfActiveLayers += 1; + if (wasDead) { + fNumberOfActiveLayers += 1; + if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1; + } } @@ -344,7 +366,12 @@ void DetectorK::RemoveLayer(char *name) { Bool_t wasDead = tmp->isDead; fLayers.Remove(tmp); fNumberOfLayers -= 1; - if (!wasDead) fNumberOfActiveLayers -= 1; + if (!wasDead) { + fNumberOfActiveLayers -= 1; + TString lname(tmp->GetName()); + if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1; + + } } } @@ -526,13 +553,22 @@ Double_t DetectorK::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, D // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account // Following is correct for 2 DOF - Double_t gamma = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level + Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2; - Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*gamma)); - + Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*c)); return ( goodHit ) ; } +Double_t DetectorK::ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) +{ + // Based on work by Ruben Shahoyen + // This is the probability to not have any match to the track (see also :ProbGoodChiSqPlusConfHit:) + + Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level + Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2; + Double_t nullHit = (1-leff+fConfLevel*leff)*TMath::Exp(-c*(alpha-1./2)); + return ( nullHit ) ; +} Double_t DetectorK::HitDensity ( Double_t radius ) { @@ -548,13 +584,13 @@ Double_t DetectorK::HitDensity ( Double_t radius ) if ( radius > fMaxRadiusSlowDet ) { - arealDensity = OneEventHitDensity(dNdEtaCent,radius) ; // Fast detectors see central collision density (only) + arealDensity = OneEventHitDensity(fdNdEtaCent,radius) ; // Fast detectors see central collision density (only) arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius) ; // Increase density due to background } if (radius < fMaxRadiusSlowDet ) { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero. - arealDensity = OneEventHitDensity(dNdEtaCent,radius) + arealDensity = OneEventHitDensity(fdNdEtaCent,radius) + IntegratedHitDensity(dNdEtaMinB,radius) + UpcHitDensity(radius) ; arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ; @@ -660,12 +696,12 @@ Double_t DetectorK::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBField)) ; Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBField)) ; - if ( pionindex >= 400 ) pionindex = 399 ; + if ( pionindex >= kNptBins ) pionindex = 399 ; if ( pionindex >= 0 ) effp = corrEfficiency[0][pionindex] ; if ( pionindex < 0 ) effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd if ( effp < 0 ) effp = 0 ; - if ( kaonindex >= 400 ) kaonindex = 399 ; + if ( kaonindex >= kNptBins ) kaonindex = 399 ; if ( kaonindex >= 0 ) effk = corrEfficiency[1][kaonindex] ; if ( kaonindex < 0 ) effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd if ( effk < 0 ) effk = 0 ; @@ -689,9 +725,9 @@ void DetectorK::SolveDOFminusOneAverage() { // Note: Obviously, does not work for the Telescope equation // - Double_t fMomentumResM[400], fResolutionRPhiM[400], fResolutionZM[400]; - Double_t efficiencyM[3][400]; - for (Int_t i=0; i<400; i++) { + Double_t fMomentumResM[kNptBins], fResolutionRPhiM[kNptBins], fResolutionZM[kNptBins]; + Double_t efficiencyM[3][kNptBins]; + for (Int_t i=0; i fast mode + + + + // Prepare Probability Kombinations + Int_t nLayer = fNumberOfActiveITSLayers; + Int_t base = 3; // null, fake, correct + + Int_t komb = TMath::Power(base,nLayer); + + TMatrixD probLay(base,fNumberOfActiveITSLayers); + TMatrixD probKomb(komb,nLayer); + for (Int_t num=0; num=0; l--) { + Int_t pow = ((Int_t)TMath::Power(base,l+1)); + probKomb(num,nLayer-1-l)=(num%pow)/((Int_t)TMath::Power(base,l)); + } + } + for ( Int_t massloop = mStart ; massloop < 3 ; massloop++ ) { // PseudoRapidity OK, used as an angle @@ -800,10 +854,9 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t //if (bigRad<61) bigRad=61; // -> min pt around 100 MeV for Bz=0.5T (don't overdo it ... ;-) ) fTransMomenta[i] = ( 0.3*bigRad*TMath::Abs(fBField)*1e-2 ) - 0.08 + TMath::Power(10,2.3*i/nPt) / 10.0 ; if (!allPt) { // just 3 points around meanPt - fTransMomenta[i] = meanPt-0.1+i*0.1; + fTransMomenta[i] = meanPt-0.01+(Double_t)(i)*0.01; } - - + // New from here ................ // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis @@ -823,7 +876,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t trPars[kZ] = 0; // Z = 0 trPars[kSnp] = 0; // track along X axis at the vertex trPars[kTgl] = TMath::Tan(lambda); // dip - trPars[kPtI] = charge/pt; // q/pt + trPars[kPtI] = charge/pt; // q/pt // // put tiny errors to propagate to the outer radius trCov[kY2] = trCov[kZ2] = trCov[kSnp2] = trCov[kTgl2] = trCov[kPtI2] = 1e-9; @@ -841,7 +894,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t trCov[kY2] = trCov[kZ2] = kLargeErr2Coord; trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir; trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI]; - // + // // printf("%d - pt %lf r%lf | %lf %lf\n",massloop,fTransMomenta[i],(last->radius)/100,momentum, d); // Set Detector-Efficiency Storage area to unity @@ -931,13 +984,14 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { printf("Number of active layers: %d\n",fNumberOfActiveLayers) ; + if (fAtLeastCorr != -1) printf("Number of combinatorics for probabilities: %d\n",komb); printf("Mass of tracked particle: %f (at pt=%5.0lf MeV)\n",fParticleMass,fTransMomenta[i]*1000); printf("Name Radius Thickness PointResOn PointResOnZ DetRes DetResZ Density Efficiency\n") ; // printOnce =0; } // print out and efficiency calculation - // for (Int_t j=fLayers.GetEntries(); j--;) { // Layer loop + Int_t iLayActive=0; for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop layer = (CylLayerK*)fLayers.At(j); @@ -952,11 +1006,13 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t if ( (!isDead && radLength >0) ) { + Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] + phiRes * phiRes ) * 100. ; // work in cm Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] + zRes * zRes ) * 100. ; // work in cm + Double_t layerEfficiency = 0; if ( EfficiencySearchFlag == 0 ) layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ; @@ -964,8 +1020,14 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ; else if ( EfficiencySearchFlag == 2 ) layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; - + TString name(layer->GetName()); + if (!name.Contains("tpc")) { + probLay(2,iLayActive)= layerEfficiency ; // Pcorr + probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; // Pnull + probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive); // Pfake + iLayActive++; + } if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue; if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { @@ -982,9 +1044,18 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t printf(" - \n"); } - if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency; - + if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency; + + + } + + if (fAtLeastCorr != -1) { + // Calculate probabilities from Kombinatorics tree ... + Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay); + fEfficiency[massloop][i] = probs[0]; // efficiency + fFake[massloop][i] = probs[1]; // fake } + /* // vertex print if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1 && radius==0) { @@ -1014,7 +1085,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t // use a weighted estimate Cw = (Cf^-1 + Cb^-1)^-1, pw = Cw (pf Cf^-1 + pb Cb^-1). // Surely, for the most extreme point, where one error matrices is infinite, this does not change anything. - Bool_t doLikeAliRoot = kFALSE; // don't do the "combined info" but do like in Aliroot + Bool_t doLikeAliRoot = 0; // don't do the "combined info" but do like in Aliroot if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { printf("- Numbers of active layer is low (%d):\n -> \"outward\" fitting done as well to get reliable eff.estimates\n", @@ -1024,7 +1095,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t // RESET Covariance Matrix ( to 10 x the estimate -> as it is done in AliExternalTrackParam) // mIstar.UnitMatrix(); // start with unity if (doLikeAliRoot) { - probTr.ResetCovariance(10.); + probTr.ResetCovariance(10); } else { // cannot do complete reset, set to very large errors trCov[kY2] = trCov[kZ2] = kLargeErr2Coord; @@ -1033,7 +1104,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t // cout<=0; j--) { // Layer loop layer = (CylLayerK*)fLayers.At(j); @@ -1160,8 +1232,14 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ; else if ( EfficiencySearchFlag == 2 ) layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; - + TString name(layer->GetName()); + if (!name.Contains("tpc")) { + probLay(2,iLayActive)= layerEfficiency ; // Pcorr + probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; // Pnull + probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive); // Pfake + iLayActive++; + } if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue; if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { @@ -1181,6 +1259,12 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency; } + if (fAtLeastCorr != -1) { + // Calculate probabilities from Kombinatorics tree ... + Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay); + fEfficiency[massloop][i] = probs[0]; // efficiency + fFake[massloop][i] = probs[1]; // fake + } } if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { printOnce = 0 ; @@ -1201,7 +1285,7 @@ TGraph * DetectorK::GetGraphMomentumResolution(Int_t color, Int_t linewidth) { // returns the momentum resolution // - TGraph *graph = new TGraph(400, fTransMomenta, fMomentumRes); + TGraph *graph = new TGraph(kNptBins, fTransMomenta, fMomentumRes); graph->SetTitle("Momentum Resolution .vs. Pt" ) ; // graph->GetXaxis()->SetRangeUser(0.,5.0) ; graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ; @@ -1231,11 +1315,11 @@ TGraph * DetectorK::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t li TGraph * graph = 0; if (axis==0) { - graph = new TGraph ( 400, fTransMomenta, fResolutionRPhi ) ; + graph = new TGraph ( kNptBins, fTransMomenta, fResolutionRPhi ) ; graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ; graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ; } else { - graph = new TGraph ( 400, fTransMomenta, fResolutionZ ) ; + graph = new TGraph ( kNptBins, fTransMomenta, fResolutionZ ) ; graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ; graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ; } @@ -1264,7 +1348,7 @@ TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, In // axis =1 ... in z // - Double_t resolution[400]; + Double_t resolution[kNptBins]; Double_t layerResolution[2]; Double_t layerRadius[2]; @@ -1292,7 +1376,7 @@ TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, In Double_t pt, momentum, thickness,aMCS ; Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); - for ( Int_t i = 0 ; i < 400 ; i++ ) { + for ( Int_t i = 0 ; i < kNptBins ; i++ ) { // Reference data as if first two layers were acting all alone pt = fTransMomenta[i] ; momentum = pt / TMath::Cos(lambda) ; // Total momentum @@ -1307,7 +1391,7 @@ TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, In - TGraph* graph = new TGraph ( 400, fTransMomenta, resolution ) ; + TGraph* graph = new TGraph ( kNptBins, fTransMomenta, resolution ) ; if (axis==0) { graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ; @@ -1342,8 +1426,8 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin // Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); - Double_t particleEfficiency[400]; // with chosen particle mass - Double_t kaonEfficiency[400], pionEfficiency[400], d0efficiency[400]; + Double_t particleEfficiency[kNptBins]; // with chosen particle mass + Double_t kaonEfficiency[kNptBins], pionEfficiency[kNptBins], d0efficiency[kNptBins]; Double_t partEfficiency[2][400]; if (particle != 0) { @@ -1351,7 +1435,7 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin Double_t doNotDecayFactor; for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon - for ( Int_t j = 0 ; j < 400 ; j++ ) { + for ( Int_t j = 0 ; j < kNptBins ; j++ ) { // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm. Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity if ( massloop == 1 ) { // KAON @@ -1367,17 +1451,17 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin } // resulting estimate of the D0 efficiency - for ( Int_t j = 0 ; j < 400 ; j++ ) { + for ( Int_t j = 0 ; j < kNptBins ; j++ ) { d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency); } } else { - for ( Int_t j = 0 ; j < 400 ; j++ ) { + for ( Int_t j = 0 ; j < kNptBins ; j++ ) { particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi; // NOTE: Decay factor (see kaon) should be included to be realiable } } - for ( Int_t j = 0 ; j < 400 ; j++ ) { + for ( Int_t j = 0 ; j < kNptBins ; j++ ) { pionEfficiency[j] *= 100; kaonEfficiency[j] *= 100; d0efficiency[j] *= 100; @@ -1386,16 +1470,16 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin TGraph * graph = 0; if (particle==0) { - graph = new TGraph ( 400, fTransMomenta, particleEfficiency ) ; // choosen mass + graph = new TGraph ( kNptBins, fTransMomenta, particleEfficiency ) ; // choosen mass graph->SetLineWidth(1); } else if (particle==1) { - graph = new TGraph ( 400, fTransMomenta, pionEfficiency ) ; + graph = new TGraph ( kNptBins, fTransMomenta, pionEfficiency ) ; graph->SetLineWidth(1); } else if (particle ==2) { - graph = new TGraph ( 400, fTransMomenta, kaonEfficiency ) ; + graph = new TGraph ( kNptBins, fTransMomenta, kaonEfficiency ) ; graph->SetLineWidth(1); } else if (particle ==3) { - graph = new TGraph ( 400, fTransMomenta, d0efficiency ) ; + graph = new TGraph ( kNptBins, fTransMomenta, d0efficiency ) ; graph->SetLineStyle(kDashed); } else return 0; @@ -1417,6 +1501,138 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin return graph; } +TGraph * DetectorK::GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth) { + // + // particle = 0 ... choosen particle (setted particleMass) + // particle = 1 ... Pion + // particle = 2 ... Kaon + // + + Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); + + Double_t particleFake[kNptBins]; // with chosen particle mass + Double_t kaonFake[kNptBins], pionFake[kNptBins]; + Double_t partFake[2][kNptBins]; + + if (particle != 0) { + // resulting Pion and Kaon efficiency scaled with overall efficiency + Double_t doNotDecayFactor; + for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon + + for ( Int_t j = 0 ; j < kNptBins ; j++ ) { + // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm. + Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity + if ( massloop == 1 ) { // KAON + doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm. + kaonFake[j] = fFake[1][j] /( doNotDecayFactor) ; + } else { // PION + pionFake[j] = fFake[0][j] ; + } + partFake[0][j] = pionFake[j]; + partFake[1][j] = kaonFake[j]; + } + } + + } else { + for ( Int_t j = 0 ; j < kNptBins ; j++ ) { + particleFake[j] = fFake[2][j]; + // NOTE: Decay factor (see kaon) should be included to be realiable + } + } + + for ( Int_t j = 0 ; j < kNptBins ; j++ ) { + pionFake[j] *= 100; + kaonFake[j] *= 100; + particleFake[j] *= 100; + } + + TGraph * graph = 0; + if (particle==0) { + graph = new TGraph ( kNptBins, fTransMomenta, particleFake ) ; // choosen mass + graph->SetLineWidth(1); + } else if (particle==1) { + graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ; + graph->SetLineWidth(1); + } else if (particle ==2) { + graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ; + graph->SetLineWidth(1); + } + + graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ; + graph->GetXaxis()->CenterTitle(); + graph->GetXaxis()->SetNoExponent(1) ; + graph->GetXaxis()->SetMoreLogLabels(1) ; + graph->GetYaxis()->SetTitle("Fake (%)") ; + graph->GetYaxis()->CenterTitle(); + + graph->SetMinimum(0.01) ; + graph->SetMaximum(100) ; + + graph->SetLineColor(color); + graph->SetMarkerColor(color); + graph->SetLineWidth(linewidth); + + return graph; +} +TGraph * DetectorK::GetGraphRecoPurity(Int_t particle,Int_t color, Int_t linewidth) { + // + // particle = 0 ... choosen particle (setted particleMass) + // particle = 1 ... Pion + // particle = 2 ... Kaon + // + + // Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); + + Double_t particleFake[kNptBins]; // with chosen particle mass + Double_t kaonFake[kNptBins], pionFake[kNptBins]; + // Double_t partFake[2][kNptBins]; + + if (particle != 0) { + cout <<" not implemented"<SetLineWidth(1); + } else if (particle==1) { + graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ; + graph->SetLineWidth(1); + } else if (particle ==2) { + graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ; + graph->SetLineWidth(1); + } + + graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ; + graph->GetXaxis()->CenterTitle(); + graph->GetXaxis()->SetNoExponent(1) ; + graph->GetXaxis()->SetMoreLogLabels(1) ; + graph->GetYaxis()->SetTitle("Purity (%)") ; + graph->GetYaxis()->CenterTitle(); + + graph->SetMinimum(0.01) ; + graph->SetMaximum(100) ; + + graph->SetLineColor(color); + graph->SetMarkerColor(color); + graph->SetLineWidth(linewidth); + + return graph; +} + TGraph* DetectorK::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) { // @@ -1492,6 +1708,12 @@ TGraph* DetectorK::GetGraph(Int_t number, Int_t color, Int_t linewidth) { return GetGraphRecoEfficiency(2, color, linewidth); // eff. kaon case 13: return GetGraphRecoEfficiency(3, color, linewidth); // eff. D0 + case 15: + return GetGraphRecoFakes(0, color, linewidth); // Fake tracked particle + case 16: + return GetGraphRecoFakes(1, color, linewidth); // Fake pion + case 17: + return GetGraphRecoFakes(2, color, linewidth); // Fake kaon default: printf(" Error: chosen graph number not valid\n"); } @@ -1499,7 +1721,7 @@ TGraph* DetectorK::GetGraph(Int_t number, Int_t color, Int_t linewidth) { } -void DetectorK::MakeAliceAllNew(Bool_t flagTPC) { +void DetectorK::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon) { // All New configuration with X0 = 0.3 and resolution = 4 microns @@ -1507,9 +1729,15 @@ void DetectorK::MakeAliceAllNew(Bool_t flagTPC) { AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation // new ideal Pixel properties? - Double_t x0 = 0.0030; - Double_t resRPhi = 0.0004; - Double_t resZ = 0.0004; + Double_t x0 = 0.0050; + Double_t resRPhi = 0.0006; + Double_t resZ = 0.0006; + + if (flagMon) { + x0 = 0.0030; + resRPhi = 0.0004; + resZ = 0.0004; + } AddLayer((char*)"ddd1", 2.2 , x0, resRPhi, resZ); AddLayer((char*)"ddd2", 3.8 , x0, resRPhi, resZ); @@ -1791,3 +2019,67 @@ Bool_t DetectorK::GetXatLabR(AliExternalTrackParam* tr,Double_t r,Double_t &x, D // return kTRUE; } + + + +Double_t* DetectorK::PrepareEffFakeKombinations(TMatrixD *probKomb, TMatrixD *probLay) { + + if (!probLay) { + printf("Error: Layer tracking efficiencies not set \n"); + return 0; + } + + TMatrixD &tProbKomb = *probKomb; + TMatrixD &tProbLay = *probLay; + + + // Int_t base = tProbLay.GetNcols(); // 3? null, fake, correct + Int_t nLayer = tProbKomb.GetNcols(); // nlayer? - number of ITS layers + Int_t komb = tProbKomb.GetNrows(); // 3^nlayer? - number of kombinations + + // Fill probabilities + + Double_t probEff =0; + Double_t probFake =0; + for (Int_t num=0; num=fkAtLeastCorr && flFake==0) { // at least correct but zero fake + Double_t probEffLayer = 1; + for (Int_t l=0; l=fAtLeastFake) { + Double_t probFakeLayer = 1; + for (Int_t l=0; l -#include -#include -#include - -/*********************************************************** - -Fast Simulation tool for Inner Tracker Systems - -original code of using the billoir technique was developed -for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov -http://rnc.lbl.gov/~jhthomas - -Changes by S. Rossegger - -April 2011 - Now uses the Kalman method (aliroot implementation) instead of the Billoir - technique ... (changes by Ruben Shahoyan) - -March 2011 - Changes to comply with the Alice Offline coding conventions - -Feb. 2011 - Improvement in "lowest pt allowed" -> now uses helix param. for calc. of a,b - - - Adding a more sophisticaed efficiency calculation which includes - the possibility to make chi2 cuts via Confidence Levels (method of Ruben Shahoyan) - plus adding 'basic' detection efficiencies per layer ... - - - Adding "ITS Stand alone" tracking capabilities via - forward+backward tracking -> Kalman smoothing is then - used for the parameter estimates (important for efficiencies) - -Jan. 2011 - Inclusion of ImpactParameter Plots (with vtx estimates) - to allow comparison with ITS performances in pp data - -Dec. 2010 - Translation into C++ class format - - Adding various Setters and Getters to build the geometry - (based on cylinders) plus handling of the layer properties - - -***********************************************************/ - -class AliExternalTrackParam; - -class DetectorK : public TNamed { - public: - DetectorK(); - DetectorK(char *name,char *title); - virtual ~DetectorK(); - - void AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes=999999, Float_t zRes=999999, Float_t eff=1.); - - void KillLayer(char *name); - void SetRadius(char *name, Float_t radius); - void SetRadiationLength(char *name, Float_t radL); - void SetResolution(char *name, Float_t phiRes=999999, Float_t zRes=999999); - void SetLayerEfficiency(char *name, Float_t eff=1.0); - void RemoveLayer(char *name); - - Float_t GetRadius(char *name); - Float_t GetRadiationLength(char *name); - Float_t GetResolution(char *name, Int_t axis=0); - Float_t GetLayerEfficiency(char *name); - - void PrintLayout(); - void PlotLayout(Int_t plotDead = kTRUE); - - void MakeAliceAllNew(Bool_t flagTPC =1); - void MakeAliceCurrent(Int_t AlignResiduals = 0, Bool_t flagTPC =1); - void AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip=1); - void RemoveTPC(); - - void SetBField(Float_t bfield) {fBField = bfield; } - Float_t GetBField() const {return fBField; } - void SetLhcUPCscale(Float_t lhcUPCscale) {fLhcUPCscale = lhcUPCscale; } - Float_t GetLhcUPCscale() const { return fLhcUPCscale; } - void SetParticleMass(Float_t particleMass) {fParticleMass = particleMass; } - Float_t GetParticleMass() const { return fParticleMass; } - void SetIntegrationTime(Float_t integrationTime) {fIntegrationTime = integrationTime; } - Float_t GetIntegrationTime() const { return fIntegrationTime; } - void SetMaxRadiusOfSlowDetectors(Float_t maxRadiusSlowDet) {fMaxRadiusSlowDet = maxRadiusSlowDet; } - Float_t GetMaxRadiusOfSlowDetectors() const { return fMaxRadiusSlowDet; } - void SetAvgRapidity(Float_t avgRapidity) {fAvgRapidity = avgRapidity; } - Float_t GetAvgRapidity() const { return fAvgRapidity; } - void SetConfidenceLevel(Float_t confLevel) {fConfLevel = confLevel; } - Float_t GetConfidenceLevel() const { return fConfLevel; } - - Float_t GetNumberOfActiveLayers() const {return fNumberOfActiveLayers; } - - void SolveViaBilloir(Int_t flagD0=1,Int_t print=1, Bool_t allPt=1, Double_t meanPt =0.750); - void SolveDOFminusOneAverage(); - - // Helper functions - Double_t ThetaMCS ( Double_t mass, Double_t RadLength, Double_t momentum ) const; - Double_t ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) ; - Double_t ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) ; - Double_t ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) ; - - // Howard W. hit distribution and convolution integral - Double_t Dist ( Double_t Z, Double_t radius ) ; - Double_t HitDensity ( Double_t radius ) ; - Double_t UpcHitDensity ( Double_t radius ) ; - Double_t IntegratedHitDensity ( Double_t multiplicity, Double_t radius ) ; - Double_t OneEventHitDensity ( Double_t multiplicity, Double_t radius ) const ; - Double_t D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) const ; - - TGraph* GetGraphMomentumResolution(Int_t color, Int_t linewidth=1); - TGraph* GetGraphPointingResolution(Int_t axis,Int_t color, Int_t linewidth=1); - TGraph* GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth=1); - - TGraph* GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth=1); - - TGraph* GetGraphRecoEfficiency(Int_t particle, Int_t color, Int_t linewidth=1); - - TGraph* GetGraph(Int_t number, Int_t color, Int_t linewidth=1); - - void MakeStandardPlots(Bool_t add =0, Int_t color=1, Int_t linewidth=1,Bool_t onlyPionEff=0); - - // method to extend AliExternalTrackParam functionality - Bool_t GetXatLabR(AliExternalTrackParam* tr,Double_t r,Double_t &x, Double_t bz, Int_t dir=0) const; - - protected: - - Int_t fNumberOfLayers; // total number of layers in the model - Int_t fNumberOfActiveLayers; // number of active layers in the model - TList fLayers; // List of layer pointers - Float_t fBField; // Magnetic Field in Tesla - Float_t fLhcUPCscale; // UltraPeripheralElectrons: scale from RHIC to LHC - Float_t fIntegrationTime; // electronics integration time - Float_t fConfLevel; // Confidence Level for the tracking - Float_t fAvgRapidity; // rapidity of the track (= mean) - Float_t fParticleMass; // Particle used for tracking. Standard: mass of pion - Double_t fMaxRadiusSlowDet; // Maximum radius for slow detectors. Fast detectors - // and only fast detectors reside outside this radius. - - enum {kMaxNumberOfDetectors = 200}; - - Double_t fTransMomenta[400]; // array of transverse momenta - Double_t fMomentumRes[400]; // array of momentum resolution - Double_t fResolutionRPhi[400]; // array of rphi resolution - Double_t fResolutionZ[400]; // array of z resolution - Double_t fDetPointRes[kMaxNumberOfDetectors][400]; // array of rphi resolution per layer - Double_t fDetPointZRes[kMaxNumberOfDetectors][400]; // array of z resolution per layer - Double_t fEfficiency[3][400]; // efficiency for different particles - - ClassDef(DetectorK,1); -}; - -#endif +#ifndef DETECTORK_H +#define DETECTORK_H + +#include +#include +#include +#include + +/*********************************************************** + +Fast Simulation tool for Inner Tracker Systems + +original code of using the billoir technique was developed +for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov +http://rnc.lbl.gov/~jhthomas + +Changes by S. Rossegger + +July 2011 - Adding the possibility of "fake calculation" and "efficiency calculation" + with a number of "at least correct clusters on the track" + Done using the complete combinatorics table with 3^nLayer track outcomes. + +April 2011 - Now uses the Kalman method (aliroot implementation) instead of the Billoir + technique ... (changes by Ruben Shahoyan) + +March 2011 - Changes to comply with the Alice Offline coding conventions + +Feb. 2011 - Improvement in "lowest pt allowed" -> now uses helix param. for calc. of a,b + + - Adding a more sophisticaed efficiency calculation which includes + the possibility to make chi2 cuts via Confidence Levels (method of Ruben Shahoyan) + plus adding 'basic' detection efficiencies per layer ... + + - Adding "ITS Stand alone" tracking capabilities via + forward+backward tracking -> Kalman smoothing is then + used for the parameter estimates (important for efficiencies) + +Jan. 2011 - Inclusion of ImpactParameter Plots (with vtx estimates) + to allow comparison with ITS performances in pp data + +Dec. 2010 - Translation into C++ class format + - Adding various Setters and Getters to build the geometry + (based on cylinders) plus handling of the layer properties + + +***********************************************************/ + +class AliExternalTrackParam; +#include + +class DetectorK : public TNamed { + + public: + + DetectorK(); + DetectorK(char *name,char *title); + virtual ~DetectorK(); + + enum {kNptBins = 100}; // less then 400 !! + + void AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes=999999, Float_t zRes=999999, Float_t eff=0.99); + + void KillLayer(char *name); + void SetRadius(char *name, Float_t radius); + void SetRadiationLength(char *name, Float_t radL); + void SetResolution(char *name, Float_t phiRes=999999, Float_t zRes=999999); + void SetLayerEfficiency(char *name, Float_t eff=1.0); + void RemoveLayer(char *name); + + Float_t GetRadius(char *name); + Float_t GetRadiationLength(char *name); + Float_t GetResolution(char *name, Int_t axis=0); + Float_t GetLayerEfficiency(char *name); + + void PrintLayout(); + void PlotLayout(Int_t plotDead = kTRUE); + + void MakeAliceAllNew(Bool_t flagTPC =1,Bool_t flagMon=1); + void MakeAliceCurrent(Int_t AlignResiduals = 0, Bool_t flagTPC =1); + void AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip=1); + void RemoveTPC(); + + void SetBField(Float_t bfield) {fBField = bfield; } + Float_t GetBField() const {return fBField; } + void SetLhcUPCscale(Float_t lhcUPCscale) {fLhcUPCscale = lhcUPCscale; } + Float_t GetLhcUPCscale() const { return fLhcUPCscale; } + void SetParticleMass(Float_t particleMass) {fParticleMass = particleMass; } + Float_t GetParticleMass() const { return fParticleMass; } + void SetIntegrationTime(Float_t integrationTime) {fIntegrationTime = integrationTime; } + Float_t GetIntegrationTime() const { return fIntegrationTime; } + void SetMaxRadiusOfSlowDetectors(Float_t maxRadiusSlowDet) {fMaxRadiusSlowDet = maxRadiusSlowDet; } + Float_t GetMaxRadiusOfSlowDetectors() const { return fMaxRadiusSlowDet; } + void SetAvgRapidity(Float_t avgRapidity) {fAvgRapidity = avgRapidity; } + Float_t GetAvgRapidity() const { return fAvgRapidity; } + void SetConfidenceLevel(Float_t confLevel) {fConfLevel = confLevel; } + Float_t GetConfidenceLevel() const { return fConfLevel; } + void SetAtLeastCorr(Int_t atLeastCorr ) {fAtLeastCorr = atLeastCorr; } + Int_t GetAtLeastCorr() const { return fAtLeastCorr; } + void SetAtLeastFake(Int_t atLeastFake ) {fAtLeastFake = atLeastFake; } + Int_t GetAtLeastFake() const { return fAtLeastFake; } + + void SetdNdEtaCent(Int_t dNdEtaCent ) {fdNdEtaCent = dNdEtaCent; } + Float_t GetdNdEtaCent() const { return fdNdEtaCent; } + + + Float_t GetNumberOfActiveLayers() const {return fNumberOfActiveLayers; } + Float_t GetNumberOfActiveITSLayers() const {return fNumberOfActiveITSLayers; } + + void SolveViaBilloir(Int_t flagD0=1,Int_t print=1, Bool_t allPt=1, Double_t meanPt =0.750); + void SolveDOFminusOneAverage(); + + // Helper functions + Double_t ThetaMCS ( Double_t mass, Double_t RadLength, Double_t momentum ) const; + Double_t ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) ; + Double_t ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) ; + Double_t ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) ; + Double_t ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) ; + + // Howard W. hit distribution and convolution integral + Double_t Dist ( Double_t Z, Double_t radius ) ; + Double_t HitDensity ( Double_t radius ) ; + Double_t UpcHitDensity ( Double_t radius ) ; + Double_t IntegratedHitDensity ( Double_t multiplicity, Double_t radius ) ; + Double_t OneEventHitDensity ( Double_t multiplicity, Double_t radius ) const ; + Double_t D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) const ; + + TGraph* GetGraphMomentumResolution(Int_t color, Int_t linewidth=1); + TGraph* GetGraphPointingResolution(Int_t axis,Int_t color, Int_t linewidth=1); + TGraph* GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth=1); + + TGraph* GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth=1); + + TGraph* GetGraphRecoEfficiency(Int_t particle, Int_t color, Int_t linewidth=1); + TGraph* GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth); + TGraph* GetGraphRecoPurity(Int_t particle,Int_t color, Int_t linewidth); + + TGraph* GetGraph(Int_t number, Int_t color, Int_t linewidth=1); + + void MakeStandardPlots(Bool_t add =0, Int_t color=1, Int_t linewidth=1,Bool_t onlyPionEff=0); + + // method to extend AliExternalTrackParam functionality + Bool_t GetXatLabR(AliExternalTrackParam* tr,Double_t r,Double_t &x, Double_t bz, Int_t dir=0) const; + + Double_t* PrepareEffFakeKombinations(TMatrixD *probKomb, TMatrixD *probLay); + + protected: + + Int_t fNumberOfLayers; // total number of layers in the model + Int_t fNumberOfActiveLayers; // number of active layers in the model + Int_t fNumberOfActiveITSLayers; // number of active ITS layers in the model + TList fLayers; // List of layer pointers + Float_t fBField; // Magnetic Field in Tesla + Float_t fLhcUPCscale; // UltraPeripheralElectrons: scale from RHIC to LHC + Float_t fIntegrationTime; // electronics integration time + Float_t fConfLevel; // Confidence Level for the tracking + Float_t fAvgRapidity; // rapidity of the track (= mean) + Float_t fParticleMass; // Particle used for tracking. Standard: mass of pion + Double_t fMaxRadiusSlowDet; // Maximum radius for slow detectors. Fast detectors + // and only fast detectors reside outside this radius. + Int_t fAtLeastCorr; // min. number of correct hits for the track to be "good" + Int_t fAtLeastFake; // min. number of fake hits for the track to be "fake" + + Int_t fdNdEtaCent; // Multiplicity + + enum {kMaxNumberOfDetectors = 200}; + + Double_t fTransMomenta[kNptBins]; // array of transverse momenta + Double_t fMomentumRes[kNptBins]; // array of momentum resolution + Double_t fResolutionRPhi[kNptBins]; // array of rphi resolution + Double_t fResolutionZ[kNptBins]; // array of z resolution + Double_t fDetPointRes[kMaxNumberOfDetectors][kNptBins]; // array of rphi resolution per layer + Double_t fDetPointZRes[kMaxNumberOfDetectors][kNptBins]; // array of z resolution per layer + Double_t fEfficiency[3][kNptBins]; // efficiency for different particles + Double_t fFake[3][kNptBins]; // fake prob for different particles + + ClassDef(DetectorK,1); +}; + +#endif diff --git a/ITS/UPGRADE/macros/FastVsSlowSim.C b/ITS/UPGRADE/macros/FastVsSlowSim.C new file mode 100644 index 00000000000..fc0cad1c35a --- /dev/null +++ b/ITS/UPGRADE/macros/FastVsSlowSim.C @@ -0,0 +1,1367 @@ +#ifndef __CINT__ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliCDBManager.h" +#include "AliCDBEntry.h" +#include "AliRawReader.h" +#include +#include "AliFMDRawReader.h" +#include "AliFMDParameters.h" +#include "AliFMDDigit.h" +#include "AliTPCAltroMapping.h" +#include "AliTriggerConfiguration.h" +#include "TAlienCollection.h" +#include "AliTriggerClass.h" +#include "TGrid.h" +#include "TPRegexp.h" +#include "AliVZERORawStream.h" + +#include "AliITSsegmentationUpgrade.h" +#include "AliRunLoader.h" +#include "AliStack.h" +#include "AliITSRecPointU.h" +#include "AliITSDigitUpgrade.h" +#include "TParticle.h" +#include "AliESDtrack.h" +#include "AliESDEvent.h" +#include "AliTrackReference.h" +#include +#include +#include "TGraphErrors.h" +#include +#include + +#endif +#include +#include +#include + +/* + .L ~/ITSupgrade/CDRpics/FastVsSlowSim.C + ExtractOutputHistos(0,1); + + + .L ~/ITSupgrade/BuildDetector/DetectorK.cxx++ + .L FastVsSlowSim.C + + FastVsSlowSimRes(); + FastVsSlowSimPtRes(); + FastVsSlowSimEff(0,1); + FastVsSlowSimEff(1,1); + FastVsSlowSimEff(2,1); + FastVsSlowSimEff(3,1); + + FastVsSlowSimEff(0); + +*/ + + +Int_t atLeastcorr = 3; + + +void FastVsSlowSim(Bool_t extract=1) { + + if (extract) + ExtractOutputHistos(0,1); + else + plotMerged(); +} + +TGraph *gr[5]; +Int_t colors[5]={1,2,3,4,6}; +Int_t width =2; + + +// new ideal Pixel properties? + +Double_t etaCut = 0.9; +Double_t X0 = 0.003; +Double_t resRPhi = 0.0004; +Double_t resZ = 0.0004; + +void FastVsSlowSimRes() { + + Int_t plusTPC =0; + + gROOT->LoadMacro("~/fig_template.C"); // figure style + myOptions(0); + gROOT->ForceStyle(); + + TCanvas *myCan = new TCanvas("myCan"); + myCan->Draw(); + myCan->cd(); + + TPad *myPad = new TPad("myPad", "The pad",0,0,1,1); + myPadSetUp(myPad,0.15,0.04,0.04,0.15); + myPad->Draw(); myPad->cd(); + myPad->SetGridx(); myPad->SetGridy(); myPad->SetLogx(); + + // TLegend *leg = new TLegend(0.7,160,20,290,"","brCDN"); + TLegend *leg = new TLegend(0.44,160,1.7,290,"","brCDN"); + + leg->SetFillColor(0); + + + // Current ITS +++++++++++++++++++++++++++++++++++++++++ + + DetectorK its("ALICE","ITS"); + its.MakeAliceCurrent(0,plusTPC); + its.SetMaxRadiusOfSlowDetectors(0.1); + its.SolveViaBilloir(0); + Int_t color=1; Int_t linewidth=2; + + TGraph *c[6]; + TGraph *d[6]; + + Int_t pi =0; + d[pi] = its.GetGraphPointingResolution(1,color,linewidth); + d[pi]->SetLineStyle(2); + // d[pi]->GetYaxis()->SetTitle("Pointing resolution #sigma [#mum]"); + // d[pi]->SetTitle("Pointing resolution .vs. Pt"); + // d[pi]->Draw("AC"); + + c[pi] = its.GetGraphPointingResolution(0,color,linewidth); + c[pi]->SetMinimum(-1); + c[pi]->Draw("AC"); + + leg->AddEntry(c[pi],"FastTool: Current ITS","l"); + // leg->AddEntry(d[pi],"in z - Current ITS","l"); + + + + + // Current ITS +++++++++++++++++++++++++++++++++++++++++ + + Int_t color=3; Int_t linewidth=2; + Int_t pi =2; + + DetectorK its("ALICE","ITS"); + its.MakeAliceCurrent(0,plusTPC); + + its.SetRadius("bpipe",2.0); + its.AddLayer("spd0", 2.2,1,1,1); + + its.SetRadius("spd0",2.2); its.SetRadiationLength("spd0",X0); its.SetResolution("spd0",resRPhi,resZ); + its.SetRadius("spd1",4.8); its.SetRadiationLength("spd1",X0); its.SetResolution("spd1",resRPhi,resZ); + its.SetRadius("spd2",9.1); its.SetRadiationLength("spd2",X0); its.SetResolution("spd2",resRPhi,resZ); + + its.SetMaxRadiusOfSlowDetectors(0.1); + its.SolveViaBilloir(0); + + d[pi] = its.GetGraphPointingResolution(1,color,linewidth); + d[pi]->SetLineStyle(2); + // d[pi]->Draw("C"); + + c[pi] = its.GetGraphPointingResolution(0,color,linewidth); + c[pi]->Draw("C"); + + leg->AddEntry(c[pi],"FastTool: \"New SPDs\"","l"); + // leg->AddEntry(d[pi],"in z - \"New SPDs\"","l"); + + + + // ALL NEW +++++++++++++++++++++++++++++++++++++++++++ + + color=2; Int_t linewidth=2; + Int_t pi =1; + + + // for a 0.8,0.2 weight configuration + + DetectorK *itsU = new DetectorK((char*)"ALICE",(char*)"ITS"); + + itsU->AddLayer((char*)"bpipe", 2.0,0.0022); // beam pipe + itsU->AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation + + itsU->AddLayer("ddd1", 2.2 , X0, resRPhi, resZ); + itsU->AddLayer("ddd2", 3.8 , X0, resRPhi, resZ); + itsU->AddLayer("ddd3", 6.8 , X0, resRPhi, resZ); + itsU->AddLayer("ddd4", 12.4 , X0, resRPhi, resZ); + itsU->AddLayer("ddd5", 23.5 , X0, resRPhi, resZ); + itsU->AddLayer("ddd6", 39.6 , X0, resRPhi, resZ); + itsU->AddLayer("ddd7", 43.0 , X0, resRPhi, resZ); + + if(plusTPC) itsU->AddTPC(0.1,0.1); + itsU->SetMaxRadiusOfSlowDetectors(0.1); + itsU->SolveViaBilloir(0); + itsU->PrintLayout(); + + + d[pi] = itsU->GetGraphPointingResolution(1,color,linewidth); + d[pi]->SetLineStyle(2); + // d[pi]->Draw("C"); + + c[pi] = itsU->GetGraphPointingResolution(0,color,linewidth); + c[pi]->SetMaximum(150); + c[pi]->Draw("C"); + + leg->AddEntry(c[pi],"FastTool: \"All New\" ","l"); + // leg->AddEntry(d[pi],"in z - \"All New\" ","l"); + + + // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + + TFile f1("root/FastVsSlow_CurrentITS-PbPb-fran.root"); + TFile f2("root/FastVsSlow_NewSPDs-PbPb-fran.root"); + TFile f3("root/FastVsSlow_AllNew-PbPb-fran.root"); + TGraphErrors *dca1 = (TGraphErrors*)f1.Get("dca"); + TGraphErrors *dca2 = (TGraphErrors*)f2.Get("dca"); + TGraphErrors *dca3 = (TGraphErrors*)f3.Get("dca"); + + dca1->SetMarkerStyle(21); dca1->SetMarkerColor(1); + dca2->SetMarkerStyle(21); dca2->SetMarkerColor(3); + dca3->SetMarkerStyle(21); dca3->SetMarkerColor(2); + + leg->AddEntry(dca1,"FullMC: Current ITS","PE"); + leg->AddEntry(dca2,"FullMC: \"New SPDs\"","PE"); + leg->AddEntry(dca3,"FullMC: \"All New\" ","PE"); + + dca1->Draw("APE"); dca1->SetMinimum(-1); dca1->SetMaximum(300); + dca2->Draw("PE"); + dca3->Draw("PE"); + c[0]->Draw("C"); + c[1]->Draw("C"); + c[2]->Draw("C"); + + leg->Draw(); + + myCan->SaveAs(Form("FastVsSlowSim-Res-%d.pdf",plusTPC)); + myCan->SaveAs(Form("FastVsSlowSim-Res-%d.eps",plusTPC)); + + +} + +// ============================================================================================== +// ============================================================================================== + +void FastVsSlowSimPtRes() { + + Int_t plusTPC =0; + + gROOT->LoadMacro("~/fig_template.C"); // figure style + myOptions(0); + gROOT->ForceStyle(); + + TCanvas *myCan = new TCanvas("myCan"); + myCan->Draw(); + myCan->cd(); + + TPad *myPad = new TPad("myPad", "The pad",0,0,1,1); + myPadSetUp(myPad,0.15,0.04,0.04,0.15); + myPad->Draw(); myPad->cd(); + myPad->SetGridx(); myPad->SetGridy(); myPad->SetLogx(); myPad->SetLogy(); + + + TLegend *leg = new TLegend(0.44,0.13,0.1.7,0.9,"","brCDN"); + leg->SetFillColor(0); + + TGraph *c[6]; + + // Current ITS +++++++++++++++++++++++++++++++++++++++++ + Int_t color=1; Int_t linewidth=2; + + Int_t pi =0; + + DetectorK its("ALICE","ITS"); + its.MakeAliceCurrent(0,plusTPC); + its.SetMaxRadiusOfSlowDetectors(0.1); + its.SolveViaBilloir(0); + Int_t color=1; Int_t linewidth=2; + + c[pi] = its.GetGraphMomentumResolution(color,linewidth); + c[pi]->Draw("AC"); + + leg->AddEntry(c[pi],"FastTool: Current ITS","l"); + + + // Current ITS +++++++++++++++++++++++++++++++++++++++++ + + Int_t color=3; Int_t linewidth=2; + Int_t pi =2; + + DetectorK its("ALICE","ITS"); + its.MakeAliceCurrent(0,plusTPC); + + its.SetRadius("bpipe",2.0); + its.AddLayer("spd0", 2.2,1,1,1); + + its.SetRadius("spd0",2.2); its.SetRadiationLength("spd0",X0); its.SetResolution("spd0",resRPhi,resZ); + its.SetRadius("spd1",4.8); its.SetRadiationLength("spd1",X0); its.SetResolution("spd1",resRPhi,resZ); + its.SetRadius("spd2",9.1); its.SetRadiationLength("spd2",X0); its.SetResolution("spd2",resRPhi,resZ); + + its.SetMaxRadiusOfSlowDetectors(0.1); + its.SolveViaBilloir(0); + + c[pi] = its.GetGraphMomentumResolution(color,linewidth); + c[pi]->Draw("C"); + + leg->AddEntry(c[pi],"FastTool: \"New SPDs\"","l"); + + + // ALL NEW +++++++++++++++++++++++++++++++++++++++++++ + + color=2; Int_t linewidth=2; + Int_t pi =1; + + + // for a 0.8,0.2 weight configuration + + DetectorK *itsU = new DetectorK((char*)"ALICE",(char*)"ITS"); + + itsU->AddLayer((char*)"bpipe", 2.0,0.0022); // beam pipe + itsU->AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation + + itsU->AddLayer("ddd1", 2.2 , X0, resRPhi, resZ); + itsU->AddLayer("ddd2", 3.8 , X0, resRPhi, resZ); + itsU->AddLayer("ddd3", 6.8 , X0, resRPhi, resZ); + itsU->AddLayer("ddd4", 12.4 , X0, resRPhi, resZ); + itsU->AddLayer("ddd5", 23.5 , X0, resRPhi, resZ); + itsU->AddLayer("ddd6", 39.6 , X0, resRPhi, resZ); + itsU->AddLayer("ddd7", 43.0 , X0, resRPhi, resZ); + + if(plusTPC) itsU->AddTPC(0.1,0.1); + itsU->SetMaxRadiusOfSlowDetectors(0.1); + itsU->SolveViaBilloir(0); + itsU->PrintLayout(); + + + c[pi] = itsU->GetGraphMomentumResolution(color,linewidth); + c[pi]->Draw("C"); + + leg->AddEntry(c[pi],"FastTool: \"All New\" ","l"); + + // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + + TFile f1("root/FastVsSlow_CurrentITS-PbPb-fran.root"); + TFile f2("root/FastVsSlow_NewSPDs-PbPb-fran.root"); + TFile f3("root/FastVsSlow_AllNew-PbPb-fran.root"); + + TGraphErrors *gpt1 = (TGraphErrors*)f1.Get("dPt"); + TGraphErrors *gpt2 = (TGraphErrors*)f2.Get("dPt"); + TGraphErrors *gpt3 = (TGraphErrors*)f3.Get("dPt"); + + gpt1->SetMarkerStyle(21); gpt1->SetMarkerColor(1); + gpt2->SetMarkerStyle(21); gpt2->SetMarkerColor(3); + gpt3->SetMarkerStyle(21); gpt3->SetMarkerColor(2); + + leg->AddEntry(gpt1,"FullMC: Current ITS","PE"); + leg->AddEntry(gpt2,"FullMC: \"New SPDs\"","PE"); + leg->AddEntry(gpt3,"FullMC: \"All New\" ","PE"); + + gpt1->Draw("APE"); gpt1->SetMinimum(0.1); gpt1->SetMaximum(20); + gpt2->Draw("PE"); + gpt3->Draw("PE"); + c[0]->Draw("C"); + c[1]->Draw("C"); + c[2]->Draw("C"); + + leg->Draw(); + + myCan->SaveAs(Form("FastVsSlowSim-PtRes-%d.pdf",plusTPC)); + myCan->SaveAs(Form("FastVsSlowSim-PtRes-%d.eps",plusTPC)); + + + +} + + + +// ============================================================================================== +// ============================================================================================== + +void FastVsSlowSimEff(Int_t id=0,Int_t PbPb=0) { + + Int_t mult = 2400; // 2800 // deducted from "Frackable" + if (PbPb) mult=2800; + + + Int_t plusTPC =0; + + gROOT->LoadMacro("~/fig_template.C"); // figure style + myOptions(0); + gROOT->ForceStyle(); + + TCanvas *myCan = new TCanvas("myCan"); + myCan->Draw(); + myCan->cd(); + + TPad *myPad = new TPad("myPad", "The pad",0,0,1,1); + myPadSetUp(myPad,0.15,0.04,0.04,0.15); + myPad->Draw(); myPad->cd(); + myPad->SetGridx(); myPad->SetGridy();// myPad->SetLogx(); + + + TLegend *leg = new TLegend(0.9,30,1.7,70,"","brCDN"); + leg->SetFillColor(0); + + TGraph *c[6]; + if (id!=2) { + + + // Current ITS +++++++++++++++++++++++++++++++++++++++++ + Int_t color=1; Int_t linewidth=2; + + Int_t pi =0; + + DetectorK its("ALICE","ITS"); + its.MakeAliceCurrent(0,plusTPC); + its.SetMaxRadiusOfSlowDetectors(0.01); + its.SetAtLeastCorr(atLeastcorr); + if (PbPb) its.SetdNdEtaCent(mult); + its.SolveViaBilloir(0); + Int_t color=1; Int_t linewidth=2; + + if (id==0) + c[pi] = its.GetGraphRecoEfficiency(0,color,linewidth); + else if (id==1) + c[pi] = its.GetGraphRecoPurity(0,color,linewidth); + else + c[pi] = its.GetGraphRecoFakes(0,color,linewidth); + + c[pi]->Draw("AC"); + + leg->AddEntry(c[pi],"FastTool: Current ITS","l"); + + + // NEW SPD +++++++++++++++++++++++++++++++++++++++++ + + Int_t color=3; Int_t linewidth=2; + Int_t pi =2; + + DetectorK its("ALICE","ITS"); + its.MakeAliceCurrent(0,plusTPC); + its.SetAtLeastCorr(atLeastcorr); + if (PbPb) its.SetdNdEtaCent(mult); + its.SetRadius("bpipe",2.0); + its.AddLayer("spd0", 2.2,1,1,1); + + its.SetRadius("spd0",2.2); its.SetRadiationLength("spd0",X0); its.SetResolution("spd0",resRPhi,resZ); + its.SetRadius("spd1",4.8); its.SetRadiationLength("spd1",X0); its.SetResolution("spd1",resRPhi,resZ); + its.SetRadius("spd2",9.1); its.SetRadiationLength("spd2",X0); its.SetResolution("spd2",resRPhi,resZ); + + its.SetMaxRadiusOfSlowDetectors(0.1); + its.SolveViaBilloir(0); + + if (id==0) + c[pi] = its.GetGraphRecoEfficiency(0,color,linewidth); + else if (id==1) + c[pi] = its.GetGraphRecoPurity(0,color,linewidth); + else + c[pi] = its.GetGraphRecoFakes(0,color,linewidth); + + c[pi]->Draw("C"); + + leg->AddEntry(c[pi],"FastTool: \"New SPDs\"","l"); + + + // ALL NEW +++++++++++++++++++++++++++++++++++++++++++ + + color=4; Int_t linewidth=2; + Int_t pi =1; + + + // for a 0.8,0.2 weight configuration + + DetectorK *itsU = new DetectorK((char*)"ALICE",(char*)"ITS"); + itsU->SetAtLeastCorr(atLeastcorr); + itsU->AddLayer((char*)"bpipe", 2.0,0.0022); // beam pipe + itsU->AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation + + itsU->AddLayer("ddd1", 2.2 , X0, resRPhi, resZ); + itsU->AddLayer("ddd2", 3.8 , X0, resRPhi, resZ); + itsU->AddLayer("ddd3", 6.8 , X0, resRPhi, resZ); + itsU->AddLayer("ddd4", 12.4 , X0, resRPhi, resZ); + itsU->AddLayer("ddd5", 23.5 , X0, resRPhi, resZ); + itsU->AddLayer("ddd6", 39.6 , X0, resRPhi, resZ); + // itsU->AddLayer("ddd6", 42.6 , X0, resRPhi, resZ); + itsU->AddLayer("ddd7", 43.0 , X0, resRPhi, resZ); + // itsU->AddLayer("ddd8", 43.4 , X0, resRPhi, resZ); + + + if (PbPb) itsU->SetdNdEtaCent(mult); + if(plusTPC) itsU->AddTPC(0.1,0.1); + itsU->SetMaxRadiusOfSlowDetectors(0.1); + itsU->SolveViaBilloir(0); + itsU->PrintLayout(); + + if (id==0) + c[pi] = itsU->GetGraphRecoEfficiency(0,color,linewidth); + else if (id==1) + c[pi] = itsU->GetGraphRecoPurity(0,color,linewidth); + else + c[pi] = itsU->GetGraphRecoFakes(0,color,linewidth); + c[pi]->Draw("C"); + + leg->AddEntry(c[pi],"FastTool: \"All New\" ","l"); + + // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + // ALL NEW - double outer layer +++++++++++++++++++++++++++++++++++++ + + color=2; Int_t linewidth=2; + Int_t pi =3; + + + // for a 0.8,0.2 weight configuration + + DetectorK *itsU = new DetectorK((char*)"ALICE",(char*)"ITS"); + itsU->SetAtLeastCorr(atLeastcorr); + itsU->AddLayer((char*)"bpipe", 2.0,0.0022); // beam pipe + itsU->AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation + + itsU->AddLayer("ddd1", 2.2 , X0, resRPhi, resZ); + itsU->AddLayer("ddd2", 3.8 , X0, resRPhi, resZ); + itsU->AddLayer("ddd3", 6.8 , X0, resRPhi, resZ); + itsU->AddLayer("ddd4", 12.4 , X0, resRPhi, resZ); + itsU->AddLayer("ddd5", 23.5 , X0, resRPhi, resZ); + itsU->AddLayer("ddd6", 39.6 , X0, resRPhi, resZ); + itsU->AddLayer("ddd8", 40.0 , X0, resRPhi, resZ); + itsU->AddLayer("ddd7", 43.0 , X0, resRPhi, resZ); + itsU->AddLayer("ddd9", 43.4 , X0, resRPhi, resZ); + + + if (PbPb) itsU->SetdNdEtaCent(mult); + if(plusTPC) itsU->AddTPC(0.1,0.1); + itsU->SetMaxRadiusOfSlowDetectors(0.1); + itsU->SolveViaBilloir(0); + itsU->PrintLayout(); + + if (id==0) + c[pi] = itsU->GetGraphRecoEfficiency(0,color,linewidth); + else if (id==1) + c[pi] = itsU->GetGraphRecoPurity(0,color,linewidth); + else + c[pi] = itsU->GetGraphRecoFakes(0,color,linewidth); + c[pi]->Draw("C"); + + leg->AddEntry(c[pi],"FastTool: \"All New\" (2x double layer)","l"); + + // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + } + + char h[100]; + if (PbPb==0) + sprintf(h,"-fran"); + else + sprintf(h,"-Anna"); + + TFile f1(Form("root/FastVsSlow_CurrentITS-PbPb%s.root",h)); + TFile f2(Form("root/FastVsSlow_NewSPDs-PbPb%s.root",h)); + TFile f3(Form("root/FastVsSlow_AllNew-PbPb%s.root",h)); + TFile f4(Form("root/FastVsSlow_AllNew-9-PbPb%s.root",h)); + + // TFile f1(Form("root/FastVsSlow_CurrentITS%s-fran.root",h)); + // TFile f2(Form("root/FastVsSlow_NewSPDs%s-fran.root",h)); + // TFile f3(Form("root/FastVsSlow_AllNew%s-fran.root",h)); + + TH1F *eff1 = 0; + TH1F *eff2 = 0; + TH1F *eff3 = 0; + TH1F *eff4 = 0; + if (id==0) { + eff1 = (TH1F*)f1.Get("efficiency"); + eff2 = (TH1F*)f2.Get("efficiency"); + eff3 = (TH1F*)f3.Get("efficiency"); + eff4 = (TH1F*)f4.Get("efficiency"); + eff1->GetYaxis()->SetTitle("efficiency (%)"); + } else if (id==1) { + eff1 = (TH1F*)f1.Get("purity"); + eff2 = (TH1F*)f2.Get("purity"); + eff3 = (TH1F*)f3.Get("purity"); + eff4 = (TH1F*)f4.Get("purity"); + eff1->GetYaxis()->SetTitle("purity (%)"); + } else if (id==2) { + eff1 = (TH1F*)f1.Get("annaEff"); + eff2 = (TH1F*)f2.Get("annaEff"); + eff3 = (TH1F*)f3.Get("annaEff"); + eff4 = (TH1F*)f4.Get("annaEff"); + eff1->GetYaxis()->SetTitle("Overall efficiency (%)"); + } else if (id==3) { + eff1 = (TH1F*)f1.Get("fake"); + eff2 = (TH1F*)f2.Get("fake"); + eff3 = (TH1F*)f3.Get("fake"); + eff4 = (TH1F*)f4.Get("fake"); + eff1->GetYaxis()->SetTitle("Fake ratio (%)"); + } + + eff1->SetMarkerStyle(21); eff1->SetMarkerColor(1); + eff2->SetMarkerStyle(21); eff2->SetMarkerColor(3); + eff3->SetMarkerStyle(21); eff3->SetMarkerColor(4); + eff4->SetMarkerStyle(21); eff4->SetMarkerColor(2); + + leg->AddEntry(eff1,"FullMC: Current ITS","PE"); + leg->AddEntry(eff2,"FullMC: \"New SPDs\"","PE"); + leg->AddEntry(eff3,"FullMC: \"All New\" ","PE"); + leg->AddEntry(eff4,"FullMC: \"All New\" (2x double layer)","PE"); + + eff1->SetMinimum(0.4); eff1->SetMaximum(100); + eff1->DrawCopy("E"); + eff2->DrawCopy("sameE"); + eff4->DrawCopy("sameE"); + eff3->DrawCopy("sameE"); + if (id!=2) { + c[0]->Draw("C"); + c[1]->Draw("C"); + c[2]->Draw("C"); + c[3]->Draw("C"); + } + eff2->DrawCopy("sameE"); + eff4->DrawCopy("sameE"); + eff3->DrawCopy("sameE"); + + + leg->Draw(); + + + + + TPaveText *pt = 0; + if (id!=3) + pt = new TPaveText(0.4,0.1,1.76,30); + else + pt = new TPaveText(0.4,70,1.76,100); + + pt->SetBorderSize(1); // no shadow + pt->SetTextFont(12); + TText *t1 = pt->AddText("FastTool settings: "); t1->SetTextFont(32); // bold + + pt->AddText(Form(" Tracked particles: Pions; Average rapidity: 0.45; dN_{ch}/d#eta = %d ",mult)); + + // pt->AddText("\"New SPDs\": layer radii: r = {2.2,4.8,9.1} cm"); + // pt->AddText("\"All New\: layer radii: r = {2.2,3.8,6.8,...} cm"); + // pt->AddText(Form(" New layer prop.: X/X_{0}=%1.1lf%%; #sigma_{r#phi,z}=%1.0lf#mum",X0*100,resZ*1e4)); + + TText *t2 = pt->AddText("FullMC settings: "); t2->SetTextFont(32); // bold + if (PbPb==0) { + pt->AddText(" Generator: AliGenHIJINGpara (parametrized PbPb event)"); + pt->AddText(" dN_{ch.pr.}/d#eta = 2070"); + pt->AddText(" Track selection: Pions, |#eta|<0.9"); + } else { + pt->AddText(" Generator: AliGenHijing (modified); #sqrt{s_{NN}} = 5.5 TeV"); + pt->AddText(" dN_{ch.pr.}/d#eta = 2410; Impactparameter range: b#in(0,5) #rightarrow central PbPb "); + pt->AddText(" Track selection: Pions, |#eta|<0.9"); + } + // pt->SetLabel("Settings"); + pt->SetTextAlign(12); + pt->SetFillColor(0); + pt->Draw(); + + + + + + if (PbPb==0) { + myCan->SaveAs(Form("FastVsSlowSim-Eff-%d.pdf",id)); + myCan->SaveAs(Form("FastVsSlowSim-Eff-%d.eps",id)); + }else{ + myCan->SaveAs(Form("FastVsSlowSim-Eff-PbPb-%d.pdf",id)); + myCan->SaveAs(Form("FastVsSlowSim-Eff-PbPb-%d.eps",id)); + } + + + + +} + + + +// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// Extraction + + +void GetDetectorRadii(TArrayD *rmin,TArrayD *rmax) { + // Loads geometry of the ITS upgrade + AliITSsegmentationUpgrade *seg=new AliITSsegmentationUpgrade(); + Int_t nlayers = seg->GetNLayers(); + rmin->Set(nlayers); rmax->Set(nlayers); + for (Int_t i=0; iAddAt(seg->GetRadius(i),i); + rmax->AddAt(seg->GetRadius(i)+seg->GetThickness(i),i); + } +} + +void PrintDetectorGeometry() { + gSystem->Load("libITSUpgradeBase"); + gSystem->Load("libITSUpgradeSim"); + + TArrayD rmin(0); + TArrayD rmax(0); + GetDetectorRadii(&rmin,&rmax); + for (Int_t i=0; iSetStyle("Plain"); + gStyle->SetPalette(1); + + const Int_t nbins=20; + Double_t ptmin=0.06;//04; + Double_t ptmax=2.0;//GeV + Double_t logxmin = TMath::Log10(ptmin); + Double_t logxmax = TMath::Log10(ptmax); + Double_t binwidth = (logxmax-logxmin)/(nbins+1); + enum {nb=nbins+1}; + Double_t xbins[nb]; + xbins[0] = ptmin; + for (Int_t i=1;i<=nbins;i++) { + xbins[i] = ptmin + TMath::Power(10,logxmin+(i)*binwidth); + // cout<GetXaxis()->SetTitle("eta"); + hMultCount->GetYaxis()->SetTitle("N/d#eta"); + + TH1F *hAllMC = new TH1F("allMC","All Tracks MC primaries",nbins,xbins); + TH1F *hAllFound = new TH1F("allFound","All Tracks found",nbins,xbins); + TH1F *hImperfect = new TH1F("imperfect","Imperfect tracks",nbins,xbins); + TH1F *hPerfect = new TH1F("perfect","Perfect tracks",nbins,xbins); + TH1F *hEff = new TH1F("efficiency","Efficiency (Perfect tracks in \"ALL MC\")",nbins,xbins); + TH1F *hFake = new TH1F("fake","Fake tracks (Inperfect tracks in \"ALL MC\")",nbins,xbins); + TH1F *hPurity = new TH1F("purity","Purity (Perfect tracks in \"All Found\")",nbins,xbins); + TH1F *hAnna = new TH1F("annaEff","AnnalisaEff ",nbins,xbins); + TH1F *hNoMCTrack = new TH1F("noMCtrack","noMCtrack ",nbins,xbins); + + TH1F *hEta = new TH1F("","",50,-2,2); + // TH1F *hEtaMC = new TH1F("","",50,-2,2); + + TH2D *h2Ddca = new TH2D("dca2D","DCAvsPt2D",nbins,xbins,50,-0.05,0.05); + TH2D *h2Dpt = new TH2D("dPt2D","dPtdvsPt2D",nbins,xbins,50,-25,25); + + // open run loader and load gAlice, kinematics and header + AliRunLoader* runLoader = AliRunLoader::Open("galice.root"); + if (!runLoader) { + Error("Check kine", "getting run loader from file %s failed", + "galice.root"); + return; + } + runLoader->LoadgAlice(); + gAlice = runLoader->GetAliRun(); + if (!gAlice) { + Error("Check kine", "no galice object found"); + return; + } + runLoader->LoadHeader(); + runLoader->LoadKinematics(); + + TFile* esdFile = TFile::Open("AliESDs.root"); + if (!esdFile || !esdFile->IsOpen()) { + Error("CheckESD", "opening ESD file %s failed", "AliESDs.root"); + return; + } + AliESDEvent *esd = new AliESDEvent(); + TTree* tree = (TTree*) esdFile->Get("esdTree"); + if (!tree) { + Error("CheckESD", "no ESD tree found"); + return; + } + esd->ReadFromTree(tree); + + Int_t nTrackTotalMC = 0; + Int_t nTrackFound = 0; + Int_t nTrackImperfect = 0; + Int_t nTrackPerfect = 0; + Int_t nNoMCTrack = 0; + + + for(Int_t iEv =0; iEvGetEntries(); iEv++){ + tree->GetEvent(iEv); + runLoader->GetEvent(iEv); + + printf("+++ event %i (of %lld) +++++++++++++++++++++++ # ESDtracks: %d \n",iEv,tree->GetEntries()-1,esd->GetNumberOfTracks()); + Int_t nESDtracks = esd->GetNumberOfTracks(); + for (Int_t iTrack = 0; iTrack < nESDtracks; iTrack++) { + AliESDtrack* track = esd->GetTrack(iTrack); + if (!(iTrack%1000)) printf("event %i: ESD track count %d (of %d)\n",iEv,iTrack,nESDtracks); + + Int_t label = track->GetLabel(); + + Int_t idx[12]; + // Int_t ncl = track->GetITSclusters(idx); + + if(label<0) { + // cout<< " ESD track label " << label; + // cout<<" ---> imperfect track (label "< track Pt: "<< track->Pt() << endl; + } + + AliStack* stack = runLoader->Stack(); + // nTrackTotalMC += stack->GetNprimary(); + + + TParticle* particle = stack->Particle(TMath::Abs(label)); + Double_t pt = track->Pt(); + + if(particle) { + + if (TMath::Abs(particle->Eta())>etaCut) continue; + + Double_t ptMC = particle->Pt(); + + // Efficiencies + if (onlyPion && TMath::Abs(particle->GetPdgCode())!=211) continue; + + if ( (!onlyPrims) || stack->IsPhysicalPrimary(TMath::Abs(label))) { + // cout<<" # clusters "<Fill(ptMC); + hEta->Fill(track->Eta()); + + if (label<0) { + nTrackImperfect++; + hImperfect->Fill(ptMC); + } else { + nTrackPerfect++; + hPerfect->Fill(ptMC); + } + + } + + + // following only for "true tracks, pions + + if(particle->Pt() < 0.001)continue; + if (TMath::Abs(particle->GetPdgCode())!=211) continue; + if (label>0) { + + // Impact parameters for Pions only + Double_t dca = track->GetD(0,0,0.5); + h2Ddca->Fill(ptMC,dca); + + // Pt resolution for Pions only + Double_t dPt = (pt-ptMC)/ptMC*100; + h2Dpt->Fill(ptMC,dPt); + } + + } else { + nNoMCTrackFound++; + hNoMCTrack->Fill(pt); + cout<<" according MC particle not found"<UnloadHeader(); + runLoader->UnloadKinematics(); + delete runLoader; + + + // Count trackable MC tracks + CountTrackableMCs(hAllMC, onlyPrims, onlyPion); + + + // Count trackable MC tracks + CountPrimaries(hMultCount); + + + + + // Get Errors right + hMultCount->Sumw2(); + hAllMC->Sumw2(); + hAllFound->Sumw2(); + hPerfect->Sumw2(); + hImperfect->Sumw2(); + h2Dpt->Sumw2(); + h2Ddca->Sumw2(); + + // -- Global efficienies + + nTrackTotalMC = hAllMC->GetEntries(); + Double_t eff = ((Double_t)nTrackPerfect)/nTrackTotalMC; + printf("-> Total number of events: %lld -> MCtracks %d -> nPerfect %d -> Eff: %3.2lf \n", + tree->GetEntries(),nTrackTotalMC,nTrackPerfect,eff); + + Double_t purity = ((Double_t)nTrackPerfect)/nTrackFound; + printf("-> Total number of events: %lld -> FoundTracks %d -> nPerfect %d -> Purity: %3.2lf \n", + tree->GetEntries(),nTrackFound,nTrackPerfect,purity); + + // Efficiencies - and normalize to 100% + + TF1 f1("f1","100+x*0",0.,1.e3); + + hPurity->Divide(hPerfect,hAllFound,1,1,"b"); + hPurity->Multiply(&f1); + hPurity->SetMarkerColor(kGreen); + hPurity->SetMarkerStyle(21); + hPurity->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + hPurity->SetStats(0); + + hPurity->GetYaxis()->SetRangeUser(0,100); + hPurity->SetTitle("Efficiency & Purity"); + + hEff->Divide(hPerfect,hAllMC,1,1,"b"); + hEff->Multiply(&f1); + hEff->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + hEff->SetMarkerColor(kBlue); + hEff->SetMarkerStyle(21); + hEff->SetStats(0); + + hFake->Divide(hImperfect,hAllMC,1,1,"b"); + hFake->Multiply(&f1); + hFake->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + hFake->SetMarkerColor(kRed); + hFake->SetMarkerStyle(21); + hFake->SetStats(0); + + + hAnna->Divide(hAllFound,hAllMC,1,1,"b"); + hAnna->Multiply(&f1); + hAnna->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + hAnna->SetMarkerColor(kBlack); + hAnna->SetMarkerStyle(21); + hAnna->SetStats(0); + + TCanvas *c1 = new TCanvas("c1","NoMCTrackFound");//,200,10,900,900); + TVirtualPad *pad = c1->cd(); + pad->SetGridx(); pad->SetGridy(); + hNoMCTrack->Draw(); + + TCanvas *c2 = new TCanvas("c2","Eff&Purity");//,200,10,900,900); + TVirtualPad *pad = c2->cd(); + pad->SetGridx(); pad->SetGridy(); + // pad->SetLogx(); + + hPurity->Draw("E"); + hEff->Draw("Same E"); + hFake->Draw("Same E"); + hAnna->Draw("Same E"); + + TLegend *leg = new TLegend(0.1,0.8,0.6,0.9);leg->SetFillColor(0); + leg->AddEntry(hPurity,"Purity (\"Perfect tracks\" within \"Found Tracks\")","PE"); + leg->AddEntry(hEff,"Efficiency (\"Perfect tracks\" within \"MC findable Tracks\")","PE"); + leg->AddEntry(hFake,"Fake (\"Inperfect tracks\" within \"MC findable Tracks\")","PE"); + leg->AddEntry(hAnna,"AnnaLisa - Efficiency (\"Found tracks\" within \"MC findable Tracks\")","PE"); + leg->Draw(); + + + if (plotFlag==1){ + hAllMC->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + hAllMC->Draw(); // MC pt distribution + hAllFound->SetLineColor(2); + hAllFound->Draw("same"); // MC pt distribution + } + + + /* + + .L ~/ITSupgrade/BuildDetector/DetectorK.cxx+ + + // All NEW + DetectorK its("ALICE","ITS"); + its.MakeAliceAllNew(0); + its.SetMaxRadiusOfSlowDetectors(0.01); + its.SolveViaBilloir(0); + TGraph *c = its.GetGraphRecoEfficiency(0,3,2); + c->Draw("C"); + + + // Current + DetectorK its("ALICE","ITS"); + its.MakeAliceCurrent(0,0); + its.SetMaxRadiusOfSlowDetectors(0.01); + its.SolveViaBilloir(0); + TGraph *c = its.GetGraphRecoEfficiency(0,4,2); + c->Draw("C"); + + */ + + TCanvas *c3 = new TCanvas("c3","impact");//,200,10,900,900); + c3->Divide(2,1); c3->cd(1); + // Impact parameter + + // Impact parameter resolution --------------- + h2Ddca->Draw("colz"); + h2Ddca->FitSlicesY() ; + TH2D *dcaM = (TH2D*)gDirectory->Get("dca2D_1"); dcaM->Draw("same"); + TH2D *dcaRMS = (TH2D*)gDirectory->Get("dca2D_2"); //dcaRMS->Draw(); + TGraphErrors *d0 = new TGraphErrors(); + for (Int_t ibin =1; ibin<=dcaRMS->GetXaxis()->GetNbins(); ibin++) { + d0->SetPoint( ibin-1,dcaRMS->GetBinCenter(ibin),dcaRMS->GetBinContent(ibin)*1e4); // microns + d0->SetPointError(ibin-1,0,dcaRMS->GetBinError(ibin)*1e4); // microns + } + d0->SetMarkerStyle(21); + d0->SetMaximum(200); d0->SetMinimum(0); + d0->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + d0->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)"); + d0->SetName("dca"); d0->SetTitle("DCAvsPt"); + + c3->cd(1); h2Ddca->Draw("surf2"); + c3->cd(2); d0->Draw("APE"); + + // PT RESOLUTION ------------ + TCanvas *c4 = new TCanvas("c4","pt resolution");//,200,10,900,900); + c4->Divide(2,1); c4->cd(1); + // Impact parameter + h2Dpt->Draw("colz"); + h2Dpt->FitSlicesY() ; + TH2D *dPtM = (TH2D*)gDirectory->Get("dPt2D_1"); dPtM->Draw("same"); + TH2D *dPtRMS = (TH2D*)gDirectory->Get("dPt2D_2"); // dPtRMS->Draw(""); + TGraphErrors *gPt = new TGraphErrors(); + for (Int_t ibin =1; ibin<=dPtRMS->GetXaxis()->GetNbins(); ibin++) { + gPt->SetPoint( ibin-1,dPtRMS->GetBinCenter(ibin),dPtRMS->GetBinContent(ibin)); + gPt->SetPointError(ibin-1,0,dPtRMS->GetBinError(ibin)); + } + gPt->SetMarkerStyle(21); + gPt->SetMaximum(20); gPt->SetMinimum(0); + gPt->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + gPt->GetYaxis()->SetTitle("relative momentum resolution (%)"); + gPt->SetName("dPt"); gPt->SetTitle("DPTvsPt"); + + c4->cd(1); h2Dpt->Draw("surf2"); + c4->cd(2); gPt->Draw("APE"); + + + // EXPORT -------- + + TFile f("histos.root","RECREATE"); + + hMultCount->Write(); + hAllMC->Write(); + hAllFound->Write(); + hImperfect->Write(); + hPerfect->Write(); + hNoMCTrack->Write(); + + hPurity->Write(); + hEff->Write(); + hFake->Write(); + hAnna->Write(); + + h2Ddca->Write(); + d0->Write(); + + h2Dpt->Write(); + gPt->Write(); + + f.Close(); + + return; + +} + +void CountTrackableMCs(TH1F *hAllMC, Bool_t onlyPrims,Bool_t onlyPion) { + + gSystem->Load("libITSUpgradeBase"); + gSystem->Load("libITSUpgradeSim"); + + // open run loader and load gAlice, kinematics and header + AliRunLoader* runLoader = AliRunLoader::Open("galice.root"); + if (!runLoader) { + Error("Check kine", "getting run loader from file %s failed", + "galice.root"); + return; + } + runLoader->LoadHeader(); + runLoader->LoadKinematics(); + runLoader->LoadTrackRefs(); + + AliLoader *dl = runLoader->GetDetectorLoader("ITS"); + + //Trackf + TTree *trackRefTree = 0x0; + TClonesArray *trackRef = new TClonesArray("AliTrackReference",1000); + + // TH1F *hRef = new TH1F("","",100,0,100); + TH1F *hR = new TH1F("","",100,0,100); + if (hAllMC==0) hAllMC = new TH1F("","",100,0.1,2); + Float_t ptmin = hAllMC->GetBinCenter(1)-hAllMC->GetBinWidth(1)/2; + Float_t ptmax = hAllMC->GetBinCenter(hAllMC->GetNbinsX())+hAllMC->GetBinWidth(hAllMC->GetNbinsX())/2; + // Int_t nAllMC = 0; + + // Detector geometry + TArrayD rmin(0); TArrayD rmax(0); + GetDetectorRadii(&rmin,&rmax); + TArrayI nLaySigs(rmin.GetSize()); + + printf("Counting trackable MC tracks ...\n"); + + for(Int_t iEv =0; iEvGetNumberOfEvents(); iEv++){ + Int_t nTrackableTracks = 0; + runLoader->GetEvent(iEv); + AliStack* stack = runLoader->Stack(); + printf("+++ event %i (of %d) +++++++++++++++++++++++ # total MCtracks: %d \n",iEv,runLoader->GetNumberOfEvents()-1,stack->GetNtrack()); + + trackRefTree=runLoader->TreeTR(); + TBranch *br = trackRefTree->GetBranch("TrackReferences"); + if(!br) { + printf("no TR branch available , exiting \n"); + return; + } + br->SetAddress(&trackRef); + + // init the trackRef tree + trackRefTree=runLoader->TreeTR(); + trackRefTree->SetBranchAddress("TrackReferences",&trackRef); + + // Count trackable MC tracks + for (Int_t iMC=0; iMCGetNtrack(); iMC++) { + + TParticle* particle = stack->Particle(iMC); + if (TMath::Abs(particle->Eta())>etaCut) continue; + if (onlyPrims && !stack->IsPhysicalPrimary(iMC)) continue; + if (onlyPion && TMath::Abs(particle->GetPdgCode())!=211) continue; + + + Bool_t isTrackable = 0; + nLaySigs.Reset(0); + + trackRefTree->GetEntry(stack->TreeKEntry(iMC)); + Int_t nref=trackRef->GetEntriesFast(); + for(Int_t iref =0; irefAt(iref); + if(!trR) continue; + if(trR->DetectorId()!=AliTrackReference::kITS) continue; + Float_t radPos = trR->R(); + hR->Fill(radPos); + for (Int_t il=0; il=rmin.At(il)-0.1 && radPos<=rmax.At(il)+0.1) { + // cout<<" in Layer "<=3) { + isTrackable =1; + // cout<Pt(); + // Double_t etaMC = particle->Eta(); + // if (ptMC>ptmin&&ptMCFill(ptMC);} + if (ptMC>ptmin) {nTrackableTracks++;hAllMC->Fill(ptMC);} + + } + + + } // entries tracks MC + printf(" -> trackable MC tracks: %d (%d)\n",nTrackableTracks,hAllMC->GetEntries()); + }//entries Events + + + hR->DrawCopy(); + hAllMC->DrawCopy(); + runLoader->UnloadHeader(); + runLoader->UnloadKinematics(); + delete runLoader; + +} + +void CountPrimaries(TH1F *hMultCount) { + + if (hMultCount==0) hMultCount = new TH1F("mult","averaged multiplicity (charg. prim)",80,-4.,4.); + + AliRunLoader *rl = AliRunLoader::Open("galice.root"); + rl->SetKineFileName("Kinematics.root"); + rl->LoadHeader(); + rl->LoadKinematics(); + Int_t nEvents = rl->GetNumberOfEvents(); + cout<< "N events "<GetEvent(iEv); + AliStack *s = rl->Stack(); + for(Int_t iP=0; iPGetNtrack(); iP++ ){ + TParticle *p = s->Particle(iP); + if (!(s->IsPhysicalPrimary(iP))) continue; + Float_t eta = p->Eta(); + if (p->Pt()>0.06) { + hMultCount->Fill(eta); + } + } + } + + hMultCount->DrawCopy(); + rl->UnloadHeader(); + rl->UnloadKinematics(); + delete rl; + + + +} + + +void plotMerged(Bool_t onlyPlot=0) { + + gStyle->SetPalette(1); + + TFile f("histoSum.root","UPDATE"); + + TH1F* hAllMC = f.Get("allMC"); + TH1F* hAllFound= f.Get("allFound"); + TH1F* hImperfect= f.Get("imperfect"); + TH1F* hPerfect= f.Get("perfect"); + TH1F* hNoMCTrack= f.Get("noMCtrack"); + + + // have to be recalculated + TH1F* hPurity = f.Get("purity"); + TH1F* hEff= f.Get("efficiency"); + TH1F* hFake= f.Get("fake"); + TH1F* hAnna= f.Get("annaEff"); + + TH2D* h2Ddca= f.Get("dca2D"); + TGraphErrors *d0= f.Get("dca"); + + TH2D* h2Dpt= f.Get("dPt2D"); + TGraphErrors *gPt= f.Get("dPt"); + + + if (!onlyPlot) { + /* // Get Errors right + hAllMC->Sumw2(); + hAllFound->Sumw2(); + hPerfect->Sumw2(); + hImperfect->Sumw2(); + h2Dpt->Sumw2(); + h2Ddca->Sumw2(); + */ + + // Efficiencies - and normalize to 100% + + TF1 f1("f1","100+x*0",0.,1.e3); + + hPurity->Divide(hPerfect,hAllFound,1,1,"b"); + hPurity->Multiply(&f1); + hPurity->SetMarkerColor(kGreen); + hPurity->SetMarkerStyle(21); + hPurity->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + hPurity->SetStats(0); + + hPurity->GetYaxis()->SetRangeUser(0,100); + hPurity->SetTitle("Efficiency & Purity"); + + hEff->Divide(hPerfect,hAllMC,1,1,"b"); + hEff->Multiply(&f1); + hEff->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + hEff->SetMarkerColor(kBlue); + hEff->SetMarkerStyle(21); + hEff->SetStats(0); + + hFake->Divide(hImperfect,hAllMC,1,1,"b"); + hFake->Multiply(&f1); + hFake->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + hFake->SetMarkerColor(kRed); + hFake->SetMarkerStyle(21); + hFake->SetStats(0); + + hAnna->Divide(hAllFound,hAllMC,1,1,"b"); + hAnna->Multiply(&f1); + hAnna->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + hAnna->SetMarkerColor(kBlack); + hAnna->SetMarkerStyle(21); + hAnna->SetStats(0); + + + // Impact parameter resolution --------------- + TCanvas *c3 = new TCanvas("c3","impact");//,200,10,900,900); + c3->Divide(2,1); c3->cd(1); + h2Ddca->DrawCopy("colz"); + h2Ddca->FitSlicesY() ; + TH2D *dcaM = (TH2D*)gDirectory->Get("dca2D_1"); dcaM->Draw("same"); + TH2D *dcaRMS = (TH2D*)gDirectory->Get("dca2D_2"); //dcaRMS->Draw(); + TGraphErrors *d0 = new TGraphErrors(); + for (Int_t ibin =1; ibin<=dcaRMS->GetXaxis()->GetNbins(); ibin++) { + d0->SetPoint( ibin-1,dcaRMS->GetBinCenter(ibin),dcaRMS->GetBinContent(ibin)*1e4); // microns + d0->SetPointError(ibin-1,0,dcaRMS->GetBinError(ibin)*1e4); // microns + } + d0->SetMarkerStyle(21); + d0->SetMaximum(200); d0->SetMinimum(0); + d0->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + d0->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)"); + d0->SetName("dca"); d0->SetTitle("DCAvsPt"); + // c3->cd(1); h2Ddca->Draw("surf2"); + c3->cd(2); d0->Draw("APE"); + + // PT RESOLUTION ------------ + TCanvas *c4 = new TCanvas("c4","pt resolution");//,200,10,900,900); + c4->Divide(2,1); c4->cd(1); + h2Dpt->DrawCopy("colz"); + h2Dpt->FitSlicesY() ; + TH2D *dPtM = (TH2D*)gDirectory->Get("dPt2D_1"); dPtM->Draw("same"); + TH2D *dPtRMS = (TH2D*)gDirectory->Get("dPt2D_2"); // dPtRMS->Draw(""); + TGraphErrors *gPt = new TGraphErrors(); + for (Int_t ibin =1; ibin<=dPtRMS->GetXaxis()->GetNbins(); ibin++) { + gPt->SetPoint( ibin-1,dPtRMS->GetBinCenter(ibin),dPtRMS->GetBinContent(ibin)); + gPt->SetPointError(ibin-1,0,dPtRMS->GetBinError(ibin)); + } + gPt->SetMarkerStyle(21); + gPt->SetMaximum(20); gPt->SetMinimum(0); + gPt->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)"); + gPt->GetYaxis()->SetTitle("relative momentum resolution (%)"); + gPt->SetName("dPt"); gPt->SetTitle("DPTvsPt"); + // c4->cd(1); h2Dpt->Draw("surf2"); + c4->cd(2); gPt->Draw("APE"); + + + // overwrite with normalized graphs + hPurity->Write(); + hEff->Write(); + hFake->Write(); + hAnna->Write(); + h2Ddca->Write(); + d0->Write(); + h2Dpt->Write(); + gPt->Write(); + + } + + // Plots + + TCanvas *c2 = new TCanvas("c2","Eff&Purity");//,200,10,900,900); + TVirtualPad *pad = c2->cd(); + pad->SetGridx(); pad->SetGridy(); + // pad->SetLogx(); + + TLegend *leg = new TLegend(0.1,0.8,0.6,0.9);leg->SetFillColor(0); + leg->AddEntry(hPurity,"Purity (\"Perfect tracks\" within \"Found Tracks\")","PE"); + leg->AddEntry(hEff,"Efficiency (\"Perfect tracks\" within \"MC findable Tracks\")","PE"); + leg->AddEntry(hFake,"Fake (\"Inperfect tracks\" within \"MC findable Tracks\")","PE"); + leg->AddEntry(hAnna,"AnnaLisa - Efficiency (\"Found tracks\" within \"MC findable Tracks\")","PE"); + + + hPurity->DrawCopy("E"); + hEff->DrawCopy("Same E"); + hFake->DrawCopy("Same E"); + hAnna->DrawCopy("Same E"); + leg->Draw(); + + c2->SaveAs("EffPlot.png"); + + f.Close(); + + + +} + -- 2.39.3