Fast Tool updates to compute efficiency with at least 3 clusters (SetAtLeastCorr...
authoramastros <amastros@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 15 Sep 2011 10:19:59 +0000 (10:19 +0000)
committeramastros <amastros@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 15 Sep 2011 10:19:59 +0000 (10:19 +0000)
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
ITS/UPGRADE/DetectorK.cxx
ITS/UPGRADE/DetectorK.h
ITS/UPGRADE/macros/FastVsSlowSim.C [new file with mode: 0644]

index a7533a2..7fba6ef 100644 (file)
 #include <TVirtualMagField.h>
 #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; i<nlayers; i++){
+
+         Double_t xsz[7]={14*1e-04, 14*1e-04, 14*1e-04, 121.2e-04,121.4*1e-04,69.3*1e-04,69.3*1e-04};
+         Double_t zsz[7]={14*1e-04, 14*1e-04, 14*1e-04,  97*1e-04,97*1e-04,2875*1e-04,2875*1e-04};
+
+         Double_t halfL[7]={21./2, 25./2., 32./2., 22.2,29.7,43.1,48.9};
+
+
+         Double_t r[7]={2.2,3.8,6.8,14.9,23.8,39.1,43.6};
+         radii.AddAt(r[i],i);
+
+         Double_t thick[7]={50.*1e-04,50.*1e-04,50.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04};
+         width.AddAt(thick[i],i);
+
+         // Copper radius & thickness
+         radiiCu.AddAt(r[i]+thick[i],i);
+         if (i<3) 
+           widthCu.AddAt(0.0036,i);
+         else
+           widthCu.AddAt(0.0150,i);//micron
+         
+         // no idea ????????
+         Int_t c[7]={1,1,1,1,1,1,1};
+         copper.AddAt(c[i],i);
+
+
+         // Silicon pixel sizes and length in Z
+
+         Int_t npixHalf[7];
+         npixHalf[i]=(Int_t)(halfL[i]/zsz[i]); 
+         Double_t HalfL[7];
+         HalfL[i]=npixHalf[i]*zsz[i];
+       
+         Int_t npixR[7];
+         npixR[i] = (Int_t)(2*TMath::Pi()*r[i]/xsz[i]);
+         Double_t xszInt[7];
+         xszInt[i]= 2*TMath::Pi()*r[i]/npixR[i];
+         xsize.AddAt(xszInt[i],i);
+         zsize.AddAt(zsz[i],i);
+
+         halfLengths.AddAt(HalfL[i],i);
+
+
+       }
+       AliITSupgrade *ITSu  = 
+         new AliITSupgrade("ITS","ITS upgrade",
+                           width,radii,halfLengths,
+                           radiiCu,widthCu,copper,
+                           bp,radiusBP, widthBP, halflengthBP );
+     
+
+       ITSu->SetFullSegmentation(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; i<nlayers; i++){
+
+         Double_t xsz[7]={14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04};
+         Double_t zsz[7]={14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04};
+
+         Double_t halfL[7]={21./2, 25./2., 32./2., 45./2, 67./2, 107./2, 114./2};
+
+
+         Double_t r[7]={2.2,3.8,6.8,12.4,23.5,39.6,43.0};
+         radii.AddAt(r[i],i);
+
+         Double_t thick[7]={50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04};
+         width.AddAt(thick[i],i);
+
+         // Copper radius & thickness
+         radiiCu.AddAt(r[i]+thick[i],i);
+         widthCu.AddAt(0.0036,i);
+         
+         // no idea ????????
+         Int_t c[7]={1,1,1,1,1,1,1};
+         copper.AddAt(c[i],i);
+
+
+         // Silicon pixel sizes and length in Z
+
+         Int_t npixHalf[7];
+         npixHalf[i]=(Int_t)(halfL[i]/zsz[i]); 
+         Double_t HalfL[7];
+         HalfL[i]=npixHalf[i]*zsz[i];
+       
+         Int_t npixR[7];
+         npixR[i] = (Int_t)(2*TMath::Pi()*r[i]/xsz[i]);
+         Double_t xszInt[7];
+         xszInt[i]= 2*TMath::Pi()*r[i]/npixR[i];
+         xsize.AddAt(xszInt[i],i);
+         zsize.AddAt(zsz[i],i);
+
+         halfLengths.AddAt(HalfL[i],i);
 
+
+       }
+       AliITSupgrade *ITSu  = 
+         new AliITSupgrade("ITS","ITS upgrade",
+                           width,radii,halfLengths,
+                           radiiCu,widthCu,copper,
+                           bp,radiusBP, widthBP, halflengthBP );
+     
+
+       ITSu->SetFullSegmentation(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; i<nlayers; i++){
+
+         Double_t xsz[8]={14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04};
+         Double_t zsz[8]={14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04};
+
+         Double_t halfL[8]={21./2, 25./2., 32./2., 45./2, 67./2, 114./2, 114./2, 114./2};
+
+
+         Double_t r[8]={2.2,3.8,6.8,12.4,23.5, 42.6,43.0,43.4};
+         radii.AddAt(r[i],i);
+
+         Double_t thick[8]={50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04};
+         width.AddAt(thick[i],i);
+
+         // Copper radius & thickness
+         radiiCu.AddAt(r[i]+thick[i],i);
+         widthCu.AddAt(0.0036,i);
+         
+         // no idea ????????
+         Int_t c[8]={1,1,1,1,1,1,1};
+         copper.AddAt(c[i],i);
+
+
+         // Silicon pixel sizes and length in Z
+
+         Int_t npixHalf[8];
+         npixHalf[i]=(Int_t)(halfL[i]/zsz[i]); 
+         Double_t HalfL[8];
+         HalfL[i]=npixHalf[i]*zsz[i];
+       
+         Int_t npixR[8];
+         npixR[i] = (Int_t)(2*TMath::Pi()*r[i]/xsz[i]);
+         Double_t xszInt[8];
+         xszInt[i]= 2*TMath::Pi()*r[i]/npixR[i];
+         xsize.AddAt(xszInt[i],i);
+         zsize.AddAt(zsz[i],i);
+
+         halfLengths.AddAt(HalfL[i],i);
+
+
+       }
+       AliITSupgrade *ITSu  = 
+         new AliITSupgrade("ITS","ITS upgrade",
+                           width,radii,halfLengths,
+                           radiiCu,widthCu,copper,
+                           bp,radiusBP, widthBP, halflengthBP );
+     
+
+       ITSu->SetFullSegmentation(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; i<nlayers; i++){
+
+         Double_t xsz[9]={14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04};
+         Double_t zsz[9]={14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04, 14*1e-04};
+
+         Double_t halfL[9]={21./2, 25./2., 32./2., 45./2, 67./2, 107./2, 107./2, 114./2, 114./2};
+
+
+         Double_t r[9]={2.2,3.8,6.8,12.4,23.5, 39.6,40.0,43.0,43.4};
+         radii.AddAt(r[i],i);
+
+         Double_t thick[9]={50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04,50.*1e-04};
+         width.AddAt(thick[i],i);
+
+         // Copper radius & thickness
+         radiiCu.AddAt(r[i]+thick[i],i);
+         widthCu.AddAt(0.0036,i);
+         
+         // no idea ????????
+         Int_t c[9]={1,1,1,1,1,1,1};
+         copper.AddAt(c[i],i);
+
+
+         // Silicon pixel sizes and length in Z
+
+         Int_t npixHalf[9];
+         npixHalf[i]=(Int_t)(halfL[i]/zsz[i]); 
+         Double_t HalfL[9];
+         HalfL[i]=npixHalf[i]*zsz[i];
+       
+         Int_t npixR[9];
+         npixR[i] = (Int_t)(2*TMath::Pi()*r[i]/xsz[i]);
+         Double_t xszInt[9];
+         xszInt[i]= 2*TMath::Pi()*r[i]/npixR[i];
+         xsize.AddAt(xszInt[i],i);
+         zsize.AddAt(zsz[i],i);
+
+         halfLengths.AddAt(HalfL[i],i);
+
+
+       }
+       AliITSupgrade *ITSu  = 
+         new AliITSupgrade("ITS","ITS upgrade",
+                           width,radii,halfLengths,
+                           radiiCu,widthCu,copper,
+                           bp,radiusBP, widthBP, halflengthBP );
+     
+
+       ITSu->SetFullSegmentation(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;
+}
+
index 75bf577..0c67188 100644 (file)
@@ -28,8 +28,8 @@ Changes by S. Rossegger -> see header file
 \r
 #define Luminosity    1.e27       // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 )\r
 #define SigmaD        6.0         // Size of the interaction diamond (cm) (LHC = 6.0 cm)\r
-#define dNdEtaMinB    950//660//950           // Multiplicity per unit Eta  (AuAu MinBias = 170, Central = 700)\r
-#define dNdEtaCent    2300 //1600//2300        // Multiplicity per unit Eta  (LHC at 5.5 TeV not known)\r
+#define dNdEtaMinB    1//950//660//950           // Multiplicity per unit Eta  (AuAu MinBias = 170, Central = 700)\r
+// #define dNdEtaCent    2300//15000 //1600//2300        // Multiplicity per unit Eta  (LHC at 5.5 TeV not known)\r
 \r
 #define CrossSectionMinB         8    // minB Cross section for event under study (PbPb MinBias ~ 8 Barns)\r
 #define AcceptanceOfTpcAndSi     1 //1//0.60 //0.35  // Assumed geometric acceptance (efficiency) of the TPC and Si detectors\r
@@ -42,6 +42,7 @@ Changes by S. Rossegger -> see header file
 #define KaonMass                 0.498  // Mass of the Kaon\r
 #define D0Mass                   1.865  // Mass of the D0\r
 \r
+//TMatrixD *probKomb; // table for efficiency kombinatorics\r
 \r
 \r
 class CylLayerK : public TNamed {\r
@@ -92,13 +93,17 @@ DetectorK::DetectorK()
   : TNamed("test_detector","detector"),\r
     fNumberOfLayers(0),\r
     fNumberOfActiveLayers(0),\r
+    fNumberOfActiveITSLayers(0),\r
     fBField(0.5),\r
     fLhcUPCscale(1.0),    \r
     fIntegrationTime(0.02), // in ms\r
     fConfLevel(0.0027),      // 0.27 % -> 3 sigma confidence\r
     fAvgRapidity(0.45),      // Avg rapidity, MCS calc is a function of crossing angle\r
     fParticleMass(0.140),    // Standard: pion mass \r
-    fMaxRadiusSlowDet(10.)\r
+    fMaxRadiusSlowDet(10.),\r
+    fAtLeastCorr(-1),     // if -1, then correct hit on all ITS layers\r
+    fAtLeastFake(1),       // if at least x fakes, track is considered fake ...\r
+    fdNdEtaCent(2300)\r
 {\r
   //\r
   // default constructor\r
@@ -111,13 +116,17 @@ DetectorK::DetectorK(char *name, char *title)
   : TNamed(name,title),\r
     fNumberOfLayers(0),\r
     fNumberOfActiveLayers(0),\r
+    fNumberOfActiveITSLayers(0),\r
     fBField(0.5),\r
     fLhcUPCscale(1.0),\r
     fIntegrationTime(0.02),  // in ms\r
     fConfLevel(0.0027),      // 0.27 % -> 3 sigma confidence\r
     fAvgRapidity(0.45),      // Avg rapidity, MCS calc is a function of crossing angle\r
     fParticleMass(0.140),     // Standard: pion mass\r
-    fMaxRadiusSlowDet(10.)\r
+    fMaxRadiusSlowDet(10.),\r
+    fAtLeastCorr(-1),     // if -1, then correct hit on all ITS layers\r
+    fAtLeastFake(1),       // if at least x fakes, track is considered fake ...\r
+    fdNdEtaCent(2300)\r
 {\r
   //\r
   // default constructor, that set the name and title\r
@@ -168,7 +177,11 @@ void DetectorK::AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRe
       \r
     }\r
     fNumberOfLayers += 1;\r
-    if (!(newLayer->isDead)) fNumberOfActiveLayers += 1;\r
+    if (!(newLayer->isDead)) {\r
+      fNumberOfActiveLayers += 1;\r
+      TString lname(newLayer->GetName());\r
+      if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1;\r
+    }\r
 \r
 \r
   } else {\r
@@ -192,6 +205,8 @@ void DetectorK::KillLayer(char *name) {
      if (!(tmp->isDead)) {\r
        tmp->isDead = kTRUE;\r
        fNumberOfActiveLayers -= 1; \r
+       TString lname(tmp->GetName());\r
+       if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;\r
      }     \r
   }\r
 }\r
@@ -273,13 +288,20 @@ void DetectorK::SetResolution(char *name, Float_t phiRes, Float_t zRes) {
     \r
     tmp->phiRes = phiRes;\r
     tmp->zRes = zRes;\r
-    \r
+    TString lname(tmp->GetName());\r
+\r
     if (zRes==RIDICULOUS && phiRes==RIDICULOUS) {\r
       tmp->isDead = kTRUE;\r
-      if (!wasDead) fNumberOfActiveLayers -= 1;\r
+      if (!wasDead) {\r
+       fNumberOfActiveLayers -= 1;\r
+       if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;\r
+      }\r
     } else {\r
       tmp->isDead = kFALSE;\r
-      if (wasDead) fNumberOfActiveLayers += 1;\r
+      if (wasDead) {\r
+       fNumberOfActiveLayers += 1;\r
+       if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1;\r
+      }\r
     }\r
 \r
 \r
@@ -344,7 +366,12 @@ void DetectorK::RemoveLayer(char *name) {
     Bool_t wasDead = tmp->isDead;\r
     fLayers.Remove(tmp);\r
     fNumberOfLayers -= 1;\r
-    if (!wasDead) fNumberOfActiveLayers -= 1;\r
+    if (!wasDead) {\r
+      fNumberOfActiveLayers -= 1;\r
+      TString lname(tmp->GetName());\r
+      if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;\r
+      \r
+    }\r
   }\r
 }\r
 \r
@@ -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 \r
   // Following is correct for 2 DOF\r
 \r
-  Double_t gamma = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level\r
+  Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level\r
   Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2; \r
-  Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*gamma));\r
-\r
+  Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*c));\r
   return ( goodHit ) ;  \r
 }\r
 \r
+Double_t DetectorK::ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) \r
+{\r
+  // Based on work by Ruben Shahoyen \r
+  // This is the probability to not have any match to the track (see also :ProbGoodChiSqPlusConfHit:)\r
+\r
+  Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level\r
+  Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2; \r
+  Double_t nullHit = (1-leff+fConfLevel*leff)*TMath::Exp(-c*(alpha-1./2));\r
+  return ( nullHit ) ;  \r
+}\r
 \r
 Double_t DetectorK::HitDensity ( Double_t radius ) \r
 {\r
@@ -548,13 +584,13 @@ Double_t DetectorK::HitDensity ( Double_t radius )
 \r
   if ( radius > fMaxRadiusSlowDet ) \r
     {\r
-      arealDensity  = OneEventHitDensity(dNdEtaCent,radius)  ; // Fast detectors see central collision density (only)\r
+      arealDensity  = OneEventHitDensity(fdNdEtaCent,radius)  ; // Fast detectors see central collision density (only)\r
       arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius)  ;  // Increase density due to background \r
     }\r
 \r
   if (radius < fMaxRadiusSlowDet )\r
     { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero.\r
-      arealDensity  = OneEventHitDensity(dNdEtaCent,radius) \r
+      arealDensity  = OneEventHitDensity(fdNdEtaCent,radius) \r
                    + IntegratedHitDensity(dNdEtaMinB,radius) \r
                    + UpcHitDensity(radius) ;\r
       arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ;  \r
@@ -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)) ; \r
     Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBField)) ; \r
       \r
-    if ( pionindex >= 400 ) pionindex = 399 ;\r
+    if ( pionindex >= kNptBins ) pionindex = 399 ;\r
     if ( pionindex >= 0 )   effp = corrEfficiency[0][pionindex] ;\r
     if ( pionindex <  0 )   effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd\r
     if ( effp < 0 )         effp = 0 ;\r
 \r
-    if ( kaonindex >= 400 ) kaonindex = 399 ;\r
+    if ( kaonindex >= kNptBins ) kaonindex = 399 ;\r
     if ( kaonindex >= 0 )   effk = corrEfficiency[1][kaonindex] ;\r
     if ( kaonindex <  0 )   effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd\r
     if ( effk < 0 )         effk = 0 ;\r
@@ -689,9 +725,9 @@ void DetectorK::SolveDOFminusOneAverage() {
   // Note: Obviously, does not work for the Telescope equation \r
   //\r
   \r
-  Double_t fMomentumResM[400], fResolutionRPhiM[400], fResolutionZM[400]; \r
-  Double_t efficiencyM[3][400];\r
-  for (Int_t i=0; i<400; i++) {\r
+  Double_t fMomentumResM[kNptBins], fResolutionRPhiM[kNptBins], fResolutionZM[kNptBins]; \r
+  Double_t efficiencyM[3][kNptBins];\r
+  for (Int_t i=0; i<kNptBins; i++) {\r
     fMomentumResM[i] = 0;   // Momentum resolution\r
     fResolutionRPhiM[i] = 0;   // Resolution in R\r
     fResolutionZM[i] = 0; // Resolution in Z\r
@@ -717,7 +753,7 @@ void DetectorK::SolveDOFminusOneAverage() {
       SolveViaBilloir(1,0); \r
 \r
       // produce sum for the mean calculation\r
-      for (Int_t i=0; i<400; i++) {\r
+      for (Int_t i=0; i<kNptBins; i++) {\r
        fMomentumResM[i] += fMomentumRes[i];   // Momentum resolution\r
        fResolutionRPhiM[i] += fResolutionRPhi[i];   // Resolution in R\r
        fResolutionZM[i] += fResolutionZ[i]; // Resolution in Z\r
@@ -732,7 +768,7 @@ void DetectorK::SolveDOFminusOneAverage() {
   }\r
   \r
   // save means in "std. Arrays"\r
-  for (Int_t i=0; i<400; i++) {\r
+  for (Int_t i=0; i<kNptBins; i++) {\r
     fMomentumRes[i] = fMomentumResM[i]/nITSLayers;   // Momentum resolution\r
     fResolutionRPhi[i] = fResolutionRPhiM[i]/nITSLayers;   // Resolution in R\r
     fResolutionZ[i] = fResolutionZM[i]/nITSLayers; // Resolution in Z\r
@@ -754,7 +790,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
   probTr.SetUseLogTermMS(kTRUE);\r
 \r
 \r
-  Int_t nPt = 400;\r
+  Int_t nPt = kNptBins;\r
   // Clean up ......\r
   for (Int_t i=0; i<kMaxNumberOfDetectors; i++) {\r
     for (Int_t j=0; j<nPt; j++) {\r
@@ -785,6 +821,24 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
   Int_t mStart =0; \r
   if (!flagD0) mStart = 2; // pion and kaon is skipped -> fast mode\r
 \r
\r
+  \r
+\r
+  // Prepare Probability Kombinations\r
+  Int_t nLayer = fNumberOfActiveITSLayers;\r
+  Int_t base = 3; // null, fake, correct\r
+\r
+  Int_t komb = TMath::Power(base,nLayer);\r
+\r
+  TMatrixD probLay(base,fNumberOfActiveITSLayers);\r
+  TMatrixD probKomb(komb,nLayer);\r
+  for (Int_t num=0; num<komb; num++) {\r
+    for (Int_t l=nLayer-1; l>=0; l--) {\r
+      Int_t pow = ((Int_t)TMath::Power(base,l+1));\r
+      probKomb(num,nLayer-1-l)=(num%pow)/((Int_t)TMath::Power(base,l));\r
+    }\r
+  }\r
+\r
   for ( Int_t massloop = mStart ; massloop < 3 ; massloop++ )  { \r
     \r
     // PseudoRapidity OK, used as an angle\r
@@ -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 ... ;-) )\r
       fTransMomenta[i] =  ( 0.3*bigRad*TMath::Abs(fBField)*1e-2 ) - 0.08 + TMath::Power(10,2.3*i/nPt) / 10.0 ; \r
       if (!allPt) { // just 3 points around meanPt\r
-       fTransMomenta[i] = meanPt-0.1+i*0.1;\r
+       fTransMomenta[i] = meanPt-0.01+(Double_t)(i)*0.01;\r
       }\r
-   \r
-\r
+  \r
       // New from here ................\r
 \r
       // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis\r
@@ -823,7 +876,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
       trPars[kZ] = 0;                         //            Z = 0 \r
       trPars[kSnp] = 0;                       //            track along X axis at the vertex\r
       trPars[kTgl] = TMath::Tan(lambda);      //            dip\r
-      trPars[kPtI] = charge/pt;               //            q/pt\r
+      trPars[kPtI] = charge/pt;               //            q/pt      \r
       //\r
       // put tiny errors to propagate to the outer radius\r
       trCov[kY2] = trCov[kZ2] = trCov[kSnp2] = trCov[kTgl2] = trCov[kPtI2] = 1e-9;\r
@@ -841,7 +894,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
       trCov[kY2]   = trCov[kZ2]   = kLargeErr2Coord; \r
       trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;\r
       trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];\r
-      //\r
+       //\r
       //      printf("%d - pt %lf r%lf | %lf %lf\n",massloop,fTransMomenta[i],(last->radius)/100,momentum, d);\r
 \r
       // Set Detector-Efficiency Storage area to unity\r
@@ -931,13 +984,14 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
   \r
       if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {\r
        printf("Number of active layers: %d\n",fNumberOfActiveLayers) ;\r
+       if (fAtLeastCorr != -1) printf("Number of combinatorics for probabilities: %d\n",komb);\r
        printf("Mass of tracked particle: %f (at pt=%5.0lf MeV)\n",fParticleMass,fTransMomenta[i]*1000);\r
        printf("Name   Radius Thickness PointResOn PointResOnZ  DetRes  DetResZ  Density Efficiency\n") ;\r
        //      printOnce =0;\r
       }\r
       \r
       // print out and efficiency calculation\r
-      // for (Int_t j=fLayers.GetEntries(); j--;) {  // Layer loop\r
+      Int_t iLayActive=0;\r
       for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) {  // Layer loop\r
        \r
        layer = (CylLayerK*)fLayers.At(j);\r
@@ -952,11 +1006,13 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
        \r
 \r
        if ( (!isDead && radLength >0) )  { \r
+\r
            Double_t rphiError  =  TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] + \r
                                                phiRes * phiRes ) * 100.  ; // work in cm\r
            Double_t zError     =  TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +\r
                                                zRes * zRes ) * 100.  ; // work in cm\r
            \r
+           \r
            Double_t layerEfficiency = 0;\r
            if ( EfficiencySearchFlag == 0 )\r
              layerEfficiency =  ProbGoodHit( radius*100, rphiError , zError  ) ;\r
@@ -964,8 +1020,14 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
              layerEfficiency =  ProbGoodChiSqHit( radius*100, rphiError , zError  ) ;\r
            else if ( EfficiencySearchFlag == 2 )\r
              layerEfficiency =  ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError  ) ;\r
-\r
+             \r
            TString name(layer->GetName());\r
+           if (!name.Contains("tpc")) {\r
+             probLay(2,iLayActive)= layerEfficiency ; // Pcorr\r
+             probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError  ) ; // Pnull\r
+             probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive);                 // Pfake\r
+             iLayActive++;    \r
+           }\r
            if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;\r
 \r
            if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {\r
@@ -982,9 +1044,18 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
                printf("        -  \n");\r
            }\r
 \r
-           if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency;\r
-\r
+           if (!name.Contains("tpc"))   fEfficiency[massloop][i] *= layerEfficiency;\r
+           \r
+           \r
+       }\r
+       \r
+       if (fAtLeastCorr != -1) {\r
+         // Calculate probabilities from Kombinatorics tree ...\r
+         Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay);\r
+         fEfficiency[massloop][i] = probs[0]; // efficiency\r
+         fFake[massloop][i] = probs[1];       // fake\r
        }\r
+\r
        /*\r
        // vertex print\r
        if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1 && radius==0) {\r
@@ -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).\r
        // Surely, for the most extreme point, where one error matrices is infinite, this does not change anything.\r
 \r
-       Bool_t doLikeAliRoot = kFALSE; // don't do the "combined info" but do like in Aliroot\r
+       Bool_t doLikeAliRoot = 0; // don't do the "combined info" but do like in Aliroot\r
 \r
        if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {        \r
          printf("- Numbers of active layer is low (%d):\n    -> \"outward\" fitting done as well to get reliable eff.estimates\n",\r
@@ -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)\r
        //      mIstar.UnitMatrix(); // start with unity\r
        if (doLikeAliRoot) {\r
-         probTr.ResetCovariance(10.);\r
+         probTr.ResetCovariance(10);\r
        } else {\r
          // cannot do complete reset, set to very large errors\r
          trCov[kY2] = trCov[kZ2] = kLargeErr2Coord; \r
@@ -1033,7 +1104,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
          //      cout<<pt<<": "<<kLargeErr2Coord<<" "<<kLargeErr2Dir<<" "<<kLargeErr2PtI*trPars[kPtI]*trPars[kPtI]<<endl;\r
        }\r
        // Clean up and storing of "forward estimates"\r
-       Double_t detPointResForw[kMaxNumberOfDetectors][400], detPointZResForw[kMaxNumberOfDetectors][400] ; \r
+       Double_t detPointResForw[kMaxNumberOfDetectors][kNptBins], detPointZResForw[kMaxNumberOfDetectors][kNptBins] ; \r
        for (Int_t k=0; k<kMaxNumberOfDetectors; k++) {\r
          for (Int_t l=0; l<nPt; l++) {\r
            detPointResForw[k][l]  = fDetPointRes[k][l];\r
@@ -1134,6 +1205,7 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
        fEfficiency[massloop][i] = 1.0 ;\r
      \r
        // print out and efficiency calculation\r
+       iLayActive=0;\r
        for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) {  // Layer loop\r
          \r
          layer = (CylLayerK*)fLayers.At(j);\r
@@ -1160,8 +1232,14 @@ void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t
              layerEfficiency =  ProbGoodChiSqHit( radius*100, rphiError , zError  ) ;\r
            else if ( EfficiencySearchFlag == 2 )\r
              layerEfficiency =  ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError  ) ;\r
-\r
+             \r
            TString name(layer->GetName());\r
+           if (!name.Contains("tpc")) {\r
+             probLay(2,iLayActive)= layerEfficiency ; // Pcorr\r
+             probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError  ) ; // Pnull\r
+             probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive);                 // Pfake\r
+             iLayActive++;    \r
+           }\r
            if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;\r
 \r
            if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {\r
@@ -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;\r
            \r
          }\r
+         if (fAtLeastCorr != -1) {\r
+           // Calculate probabilities from Kombinatorics tree ...\r
+           Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay);\r
+           fEfficiency[massloop][i] = probs[0]; // efficiency\r
+           fFake[massloop][i] = probs[1];       // fake\r
+         }\r
        }\r
        if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {\r
          printOnce = 0 ;\r
@@ -1201,7 +1285,7 @@ TGraph * DetectorK::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {
   // returns the momentum resolution \r
   //\r
   \r
-  TGraph *graph = new TGraph(400, fTransMomenta, fMomentumRes);\r
+  TGraph *graph = new TGraph(kNptBins, fTransMomenta, fMomentumRes);\r
   graph->SetTitle("Momentum Resolution .vs. Pt" ) ;\r
   //  graph->GetXaxis()->SetRangeUser(0.,5.0) ;\r
   graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
@@ -1231,11 +1315,11 @@ TGraph * DetectorK::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t li
   TGraph * graph =  0;\r
 \r
   if (axis==0) {\r
-    graph = new TGraph ( 400, fTransMomenta, fResolutionRPhi ) ;\r
+    graph = new TGraph ( kNptBins, fTransMomenta, fResolutionRPhi ) ;\r
     graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;\r
     graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;\r
   } else {\r
-    graph =  new TGraph ( 400, fTransMomenta, fResolutionZ ) ;\r
+    graph =  new TGraph ( kNptBins, fTransMomenta, fResolutionZ ) ;\r
     graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;\r
     graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;\r
   }\r
@@ -1264,7 +1348,7 @@ TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, In
   // axis =1 ... in z\r
   //\r
   \r
-  Double_t resolution[400];\r
+  Double_t resolution[kNptBins];\r
 \r
   Double_t layerResolution[2];\r
   Double_t layerRadius[2];\r
@@ -1292,7 +1376,7 @@ TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, In
   Double_t pt, momentum, thickness,aMCS ;\r
   Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
 \r
-  for ( Int_t i = 0 ; i < 400 ; i++ ) { \r
+  for ( Int_t i = 0 ; i < kNptBins ; i++ ) { \r
     // Reference data as if first two layers were acting all alone \r
     pt  =  fTransMomenta[i]  ;\r
     momentum = pt / TMath::Cos(lambda)   ;  // Total momentum\r
@@ -1307,7 +1391,7 @@ TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, In
 \r
 \r
 \r
-  TGraph* graph = new TGraph ( 400, fTransMomenta, resolution ) ;\r
+  TGraph* graph = new TGraph ( kNptBins, fTransMomenta, resolution ) ;\r
    \r
   if (axis==0) {\r
     graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ;\r
@@ -1342,8 +1426,8 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin
   //\r
   Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
   \r
-  Double_t particleEfficiency[400]; // with chosen particle mass\r
-  Double_t kaonEfficiency[400], pionEfficiency[400], d0efficiency[400]; \r
+  Double_t particleEfficiency[kNptBins]; // with chosen particle mass\r
+  Double_t kaonEfficiency[kNptBins], pionEfficiency[kNptBins], d0efficiency[kNptBins]; \r
   Double_t partEfficiency[2][400];\r
   \r
   if (particle != 0) {\r
@@ -1351,7 +1435,7 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin
     Double_t doNotDecayFactor;\r
     for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon\r
       \r
-      for ( Int_t j = 0 ; j < 400 ; j++ ) { \r
+      for ( Int_t j = 0 ; j < kNptBins ; j++ ) { \r
        // JT Test Let the kaon decay.  If it decays inside the TPC ... then it is gone; for all decays < 130 cm.\r
        Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda)           ;  // Total momentum at average rapidity\r
        if ( massloop == 1 ) { // KAON\r
@@ -1367,17 +1451,17 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin
     }\r
     \r
     // resulting estimate of the D0 efficiency\r
-    for ( Int_t j = 0 ; j < 400 ; j++ ) {\r
+    for ( Int_t j = 0 ; j < kNptBins ; j++ ) {\r
       d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);\r
     }\r
   } else { \r
-    for ( Int_t j = 0 ; j < 400 ; j++ ) { \r
+    for ( Int_t j = 0 ; j < kNptBins ; j++ ) { \r
       particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;\r
       // NOTE: Decay factor (see kaon) should be included to be realiable\r
     }\r
   }\r
 \r
-  for ( Int_t j = 0 ; j < 400 ; j++ ) { \r
+  for ( Int_t j = 0 ; j < kNptBins ; j++ ) { \r
     pionEfficiency[j]     *= 100;\r
     kaonEfficiency[j]     *= 100;\r
     d0efficiency[j]       *= 100;\r
@@ -1386,16 +1470,16 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin
  \r
   TGraph * graph =  0;\r
   if (particle==0) {\r
-    graph = new TGraph ( 400, fTransMomenta, particleEfficiency ) ; // choosen mass\r
+    graph = new TGraph ( kNptBins, fTransMomenta, particleEfficiency ) ; // choosen mass\r
     graph->SetLineWidth(1);\r
   }  else if (particle==1) {\r
-    graph = new TGraph ( 400, fTransMomenta, pionEfficiency ) ;\r
+    graph = new TGraph ( kNptBins, fTransMomenta, pionEfficiency ) ;\r
     graph->SetLineWidth(1);\r
   }  else if (particle ==2) {\r
-    graph = new TGraph ( 400, fTransMomenta, kaonEfficiency ) ;\r
+    graph = new TGraph ( kNptBins, fTransMomenta, kaonEfficiency ) ;\r
     graph->SetLineWidth(1);\r
   }  else if (particle ==3) {\r
-    graph = new TGraph ( 400, fTransMomenta, d0efficiency ) ;\r
+    graph = new TGraph ( kNptBins, fTransMomenta, d0efficiency ) ;\r
     graph->SetLineStyle(kDashed);\r
   } else \r
     return 0;\r
@@ -1417,6 +1501,138 @@ TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t lin
   return graph;\r
 }\r
 \r
+TGraph * DetectorK::GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth) {\r
+  //\r
+  // particle = 0 ... choosen particle (setted particleMass)\r
+  // particle = 1 ... Pion\r
+  // particle = 2 ... Kaon\r
+  //\r
+\r
+  Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
+  \r
+  Double_t particleFake[kNptBins]; // with chosen particle mass\r
+  Double_t kaonFake[kNptBins], pionFake[kNptBins];\r
+  Double_t partFake[2][kNptBins];\r
+  \r
+  if (particle != 0) {\r
+    // resulting Pion and Kaon efficiency scaled with overall efficiency\r
+    Double_t doNotDecayFactor;\r
+    for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon\r
+      \r
+      for ( Int_t j = 0 ; j < kNptBins ; j++ ) { \r
+       // JT Test Let the kaon decay.  If it decays inside the TPC ... then it is gone; for all decays < 130 cm.\r
+       Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda)           ;  // Total momentum at average rapidity\r
+       if ( massloop == 1 ) { // KAON\r
+         doNotDecayFactor  = TMath::Exp(-130/(371*momentum/KaonMass)) ;  // Decay length for kaon is 371 cm.\r
+         kaonFake[j] = fFake[1][j] /( doNotDecayFactor) ;\r
+       } else { // PION\r
+         pionFake[j] = fFake[0][j] ;   \r
+       }\r
+       partFake[0][j] = pionFake[j];\r
+       partFake[1][j] = kaonFake[j];\r
+      }      \r
+    }\r
+    \r
+  } else { \r
+    for ( Int_t j = 0 ; j < kNptBins ; j++ ) { \r
+      particleFake[j] = fFake[2][j];\r
+      // NOTE: Decay factor (see kaon) should be included to be realiable\r
+    }\r
+  }\r
+\r
+  for ( Int_t j = 0 ; j < kNptBins ; j++ ) { \r
+    pionFake[j]     *= 100;\r
+    kaonFake[j]     *= 100;\r
+    particleFake[j] *= 100;\r
+  }\r
\r
+  TGraph * graph =  0;\r
+  if (particle==0) {\r
+    graph = new TGraph ( kNptBins, fTransMomenta, particleFake ) ; // choosen mass\r
+    graph->SetLineWidth(1);\r
+  }  else if (particle==1) {\r
+    graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ;\r
+    graph->SetLineWidth(1);\r
+  }  else if (particle ==2) {\r
+    graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ;\r
+    graph->SetLineWidth(1);\r
+  } \r
+  \r
+  graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
+  graph->GetXaxis()->CenterTitle();\r
+  graph->GetXaxis()->SetNoExponent(1) ;\r
+  graph->GetXaxis()->SetMoreLogLabels(1) ;\r
+  graph->GetYaxis()->SetTitle("Fake (%)") ;\r
+  graph->GetYaxis()->CenterTitle();\r
+         \r
+  graph->SetMinimum(0.01) ; \r
+  graph->SetMaximum(100)  ; \r
+\r
+  graph->SetLineColor(color);\r
+  graph->SetMarkerColor(color);\r
+  graph->SetLineWidth(linewidth);\r
+\r
+  return graph;\r
+}\r
+TGraph * DetectorK::GetGraphRecoPurity(Int_t particle,Int_t color, Int_t linewidth) {\r
+  //\r
+  // particle = 0 ... choosen particle (setted particleMass)\r
+  // particle = 1 ... Pion\r
+  // particle = 2 ... Kaon\r
+  //\r
+\r
+  //  Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
+  \r
+  Double_t particleFake[kNptBins]; // with chosen particle mass\r
+  Double_t kaonFake[kNptBins], pionFake[kNptBins];\r
+  //  Double_t partFake[2][kNptBins];\r
+  \r
+  if (particle != 0) {\r
+    cout <<" not implemented"<<endl;\r
+      \r
+  } else { \r
+    for ( Int_t j = 0 ; j < kNptBins ; j++ ) { \r
+      particleFake[j] = fFake[2][j];\r
+      // NOTE: Decay factor (see kaon) should be included to be realiable\r
+    }\r
+  }\r
+\r
+  // Get Purity\r
+  for ( Int_t j = 0 ; j < kNptBins ; j++ ) { \r
+    pionFake[j]     = (1-pionFake[j])*100;\r
+    kaonFake[j]     = (1-kaonFake[j])*100;\r
+    particleFake[j] = (1-particleFake[j])*100;\r
+  }\r
\r
+  TGraph * graph =  0;\r
+  if (particle==0) {\r
+    graph = new TGraph ( kNptBins, fTransMomenta, particleFake ) ; // choosen mass\r
+    graph->SetLineWidth(1);\r
+  }  else if (particle==1) {\r
+    graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ;\r
+    graph->SetLineWidth(1);\r
+  }  else if (particle ==2) {\r
+    graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ;\r
+    graph->SetLineWidth(1);\r
+  } \r
+  \r
+  graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
+  graph->GetXaxis()->CenterTitle();\r
+  graph->GetXaxis()->SetNoExponent(1) ;\r
+  graph->GetXaxis()->SetMoreLogLabels(1) ;\r
+  graph->GetYaxis()->SetTitle("Purity (%)") ;\r
+  graph->GetYaxis()->CenterTitle();\r
+         \r
+  graph->SetMinimum(0.01) ; \r
+  graph->SetMaximum(100)  ; \r
+\r
+  graph->SetLineColor(color);\r
+  graph->SetMarkerColor(color);\r
+  graph->SetLineWidth(linewidth);\r
+\r
+  return graph;\r
+}\r
+\r
 \r
 TGraph* DetectorK::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) {\r
   //\r
@@ -1492,6 +1708,12 @@ TGraph* DetectorK::GetGraph(Int_t number, Int_t color, Int_t linewidth) {
     return GetGraphRecoEfficiency(2, color, linewidth);  // eff. kaon\r
   case 13: \r
     return GetGraphRecoEfficiency(3, color, linewidth);  // eff. D0\r
+  case 15:\r
+    return GetGraphRecoFakes(0, color, linewidth);  // Fake tracked particle\r
+  case 16:\r
+    return GetGraphRecoFakes(1, color, linewidth);  // Fake pion\r
+  case 17:\r
+    return GetGraphRecoFakes(2, color, linewidth);  // Fake kaon\r
   default:\r
     printf(" Error: chosen graph number not valid\n");\r
   }\r
@@ -1499,7 +1721,7 @@ TGraph* DetectorK::GetGraph(Int_t number, Int_t color, Int_t linewidth) {
 \r
 }\r
 \r
-void DetectorK::MakeAliceAllNew(Bool_t flagTPC) {\r
+void DetectorK::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon) {\r
   \r
   // All New configuration with X0 = 0.3 and resolution = 4 microns\r
   \r
@@ -1507,9 +1729,15 @@ void DetectorK::MakeAliceAllNew(Bool_t flagTPC) {
   AddLayer((char*)"vertex",     0,     0); // dummy vertex for matrix calculation\r
 \r
   // new ideal Pixel properties?\r
-  Double_t x0     = 0.0030;\r
-  Double_t resRPhi = 0.0004;\r
-  Double_t resZ   = 0.0004;\r
+  Double_t x0     = 0.0050;\r
+  Double_t resRPhi = 0.0006;\r
+  Double_t resZ   = 0.0006;\r
\r
+  if (flagMon) {\r
+    x0     = 0.0030;\r
+    resRPhi = 0.0004;\r
+    resZ   = 0.0004;\r
+  }\r
   \r
   AddLayer((char*)"ddd1",  2.2 ,  x0, resRPhi, resZ); \r
   AddLayer((char*)"ddd2",  3.8 ,  x0, resRPhi, resZ); \r
@@ -1791,3 +2019,67 @@ Bool_t DetectorK::GetXatLabR(AliExternalTrackParam* tr,Double_t r,Double_t &x, D
   //\r
   return kTRUE;\r
 }\r
+\r
+\r
+\r
+Double_t* DetectorK::PrepareEffFakeKombinations(TMatrixD *probKomb, TMatrixD *probLay) {\r
+\r
+  if (!probLay) {  \r
+    printf("Error: Layer tracking efficiencies not set \n");\r
+    return 0;\r
+  }\r
+\r
+  TMatrixD &tProbKomb = *probKomb;\r
+  TMatrixD &tProbLay = *probLay;\r
+\r
+\r
+  //  Int_t base = tProbLay.GetNcols(); // 3? null, fake, correct\r
+  Int_t nLayer = tProbKomb.GetNcols(); // nlayer? - number of ITS layers\r
+  Int_t komb = tProbKomb.GetNrows(); // 3^nlayer? - number of kombinations\r
+\r
+  // Fill probabilities \r
+\r
+  Double_t probEff =0;\r
+  Double_t probFake =0;\r
+  for (Int_t num=0; num<komb; num++) {\r
+    Int_t flCorr=0, flFake=0, flNull=0; \r
+     for (Int_t l=0; l<nLayer; l++)  {\r
+      if (tProbKomb(num,l)==0) \r
+       flNull++;\r
+      else if (tProbKomb(num,l)==1)\r
+       flFake++;\r
+      else if (tProbKomb(num,l)==2)\r
+       flCorr++;\r
+      else \r
+       printf("Error: unexpected values in combinatorics table\n");\r
+    }\r
+\r
+    Int_t fkAtLeastCorr = fAtLeastCorr;\r
+    if (fAtLeastCorr == -1) fkAtLeastCorr = nLayer; // all hits are "correct"\r
+    \r
+    if (flCorr>=fkAtLeastCorr && flFake==0) { // at least correct but zero fake\r
+      Double_t probEffLayer = 1;\r
+      for (Int_t l=0; l<nLayer; l++) {\r
+       probEffLayer *=  tProbLay(tProbKomb(num,l),l);\r
+       //      cout<<a(num,l)<<" ";\r
+      }\r
+      //      cout<<endl;\r
+      probEff+=probEffLayer;\r
+    }\r
\r
+    if (flFake>=fAtLeastFake) {\r
+      Double_t probFakeLayer = 1;\r
+      for (Int_t l=0; l<nLayer; l++) {\r
+       probFakeLayer *=  tProbLay(tProbKomb(num,l),l);\r
+       //      cout<<a(num,l)<<" ";\r
+      }\r
+      //      cout<<endl;\r
+      probFake+=probFakeLayer;\r
+    }\r
+\r
+  }\r
+  Double_t *probs = new Double_t(2);\r
+  probs[0] = probEff; probs[1] = probFake;\r
+  return probs;\r
+\r
+}\r
index baa482f..e7b1d50 100644 (file)
-#ifndef DETECTORK_H
-#define DETECTORK_H
-
-#include <TNamed.h>
-#include <TList.h>
-#include <TGraph.h>
-#include <Riostream.h>
-
-/***********************************************************
-
-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\r
+#define DETECTORK_H\r
+\r
+#include <TNamed.h>\r
+#include <TList.h>\r
+#include <TGraph.h>\r
+#include <Riostream.h>\r
+\r
+/***********************************************************\r
+\r
+Fast Simulation tool for Inner Tracker Systems\r
+\r
+original code of using the billoir technique was developed\r
+for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov\r
+http://rnc.lbl.gov/~jhthomas\r
+\r
+Changes by S. Rossegger\r
+\r
+July 2011  - Adding the possibility of "fake calculation" and "efficiency calculation" \r
+             with a number of "at least correct clusters on the track"\r
+             Done using the complete combinatorics table with 3^nLayer track outcomes.\r
+\r
+April 2011 - Now uses the Kalman method (aliroot implementation) instead of the Billoir\r
+             technique ... (changes by Ruben Shahoyan)\r
+\r
+March 2011 - Changes to comply with the Alice Offline coding conventions\r
+\r
+Feb. 2011 - Improvement in "lowest pt allowed" -> now uses helix param. for calc. of a,b\r
+\r
+          - Adding a more sophisticaed efficiency calculation which includes\r
+            the possibility to make chi2 cuts via Confidence Levels (method of Ruben Shahoyan)\r
+            plus adding 'basic' detection efficiencies per layer ...\r
+\r
+          - Adding "ITS Stand alone" tracking capabilities via \r
+            forward+backward tracking -> Kalman smoothing is then \r
+            used for the parameter estimates (important for efficiencies)\r
+\r
+Jan. 2011 - Inclusion of ImpactParameter Plots (with vtx estimates)\r
+            to allow comparison with ITS performances in pp data\r
+\r
+Dec. 2010 - Translation into C++ class format \r
+          - Adding various Setters and Getters to build the geometry \r
+           (based on cylinders) plus handling of the layer properties \r
+\r
+\r
+***********************************************************/\r
+\r
+class AliExternalTrackParam; \r
+#include <TMatrixD.h>\r
+\r
+class DetectorK : public TNamed {\r
+\r
+ public:\r
+  \r
+  DetectorK();\r
+  DetectorK(char *name,char *title);\r
+  virtual ~DetectorK();\r
+\r
+  enum {kNptBins = 100}; // less then 400 !!\r
\r
+  void AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes=999999, Float_t zRes=999999, Float_t eff=0.99);\r
+\r
+  void KillLayer(char *name);\r
+  void SetRadius(char *name, Float_t radius);\r
+  void SetRadiationLength(char *name, Float_t radL);\r
+  void SetResolution(char *name, Float_t phiRes=999999, Float_t zRes=999999);\r
+  void SetLayerEfficiency(char *name, Float_t eff=1.0);\r
+  void RemoveLayer(char *name);\r
+\r
+  Float_t GetRadius(char *name);\r
+  Float_t GetRadiationLength(char *name);\r
+  Float_t GetResolution(char *name, Int_t axis=0);\r
+  Float_t GetLayerEfficiency(char *name);\r
+\r
+  void PrintLayout(); \r
+  void PlotLayout(Int_t plotDead = kTRUE);\r
+  \r
+  void MakeAliceAllNew(Bool_t flagTPC =1,Bool_t flagMon=1);\r
+  void MakeAliceCurrent(Int_t AlignResiduals = 0, Bool_t flagTPC =1);\r
+  void AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip=1);\r
+  void RemoveTPC();\r
+\r
+  void SetBField(Float_t bfield) {fBField = bfield; }\r
+  Float_t GetBField() const {return fBField; }\r
+  void SetLhcUPCscale(Float_t lhcUPCscale) {fLhcUPCscale = lhcUPCscale; }\r
+  Float_t GetLhcUPCscale() const { return fLhcUPCscale; }\r
+  void SetParticleMass(Float_t particleMass) {fParticleMass = particleMass; }\r
+  Float_t GetParticleMass() const { return fParticleMass; }\r
+  void SetIntegrationTime(Float_t integrationTime) {fIntegrationTime = integrationTime; }\r
+  Float_t GetIntegrationTime() const { return fIntegrationTime; }\r
+  void SetMaxRadiusOfSlowDetectors(Float_t maxRadiusSlowDet) {fMaxRadiusSlowDet =  maxRadiusSlowDet; }\r
+  Float_t GetMaxRadiusOfSlowDetectors() const { return fMaxRadiusSlowDet; }\r
+  void SetAvgRapidity(Float_t avgRapidity) {fAvgRapidity = avgRapidity; }\r
+  Float_t GetAvgRapidity() const { return fAvgRapidity; }\r
+  void SetConfidenceLevel(Float_t confLevel) {fConfLevel = confLevel; }\r
+  Float_t GetConfidenceLevel() const { return fConfLevel; }\r
+  void SetAtLeastCorr(Int_t atLeastCorr ) {fAtLeastCorr = atLeastCorr; }\r
+  Int_t GetAtLeastCorr() const { return fAtLeastCorr; }\r
+  void SetAtLeastFake(Int_t atLeastFake ) {fAtLeastFake = atLeastFake; }\r
+  Int_t GetAtLeastFake() const { return fAtLeastFake; }\r
+\r
+  void SetdNdEtaCent(Int_t dNdEtaCent ) {fdNdEtaCent = dNdEtaCent; }\r
+  Float_t GetdNdEtaCent() const { return fdNdEtaCent; }\r
+\r
+  \r
+  Float_t GetNumberOfActiveLayers() const {return fNumberOfActiveLayers; }\r
+  Float_t GetNumberOfActiveITSLayers() const {return fNumberOfActiveITSLayers; }\r
+\r
+  void SolveViaBilloir(Int_t flagD0=1,Int_t print=1, Bool_t allPt=1, Double_t meanPt =0.750);\r
+  void SolveDOFminusOneAverage();\r
+\r
+  // Helper functions\r
+  Double_t ThetaMCS                 ( Double_t mass, Double_t RadLength, Double_t momentum ) const;\r
+  Double_t ProbGoodHit              ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; \r
+  Double_t ProbGoodChiSqHit         ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; \r
+  Double_t ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; \r
+  Double_t ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; \r
\r
+  // Howard W. hit distribution and convolution integral\r
+  Double_t Dist              ( Double_t Z, Double_t radius ) ;  \r
+  Double_t HitDensity        ( Double_t radius )   ;\r
+  Double_t UpcHitDensity     ( Double_t radius )   ;\r
+  Double_t IntegratedHitDensity  ( Double_t multiplicity, Double_t radius )   ;\r
+  Double_t OneEventHitDensity    ( Double_t multiplicity, Double_t radius ) const   ;\r
+  Double_t D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) const ;\r
+  \r
+  TGraph* GetGraphMomentumResolution(Int_t color, Int_t linewidth=1);\r
+  TGraph* GetGraphPointingResolution(Int_t axis,Int_t color, Int_t linewidth=1);\r
+  TGraph* GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth=1);\r
+\r
+  TGraph* GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth=1);\r
+\r
+  TGraph* GetGraphRecoEfficiency(Int_t particle, Int_t color, Int_t linewidth=1); \r
+  TGraph* GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth);\r
+  TGraph* GetGraphRecoPurity(Int_t particle,Int_t color, Int_t linewidth);\r
+\r
+  TGraph* GetGraph(Int_t number, Int_t color, Int_t linewidth=1);\r
+\r
+  void MakeStandardPlots(Bool_t add =0, Int_t color=1, Int_t linewidth=1,Bool_t onlyPionEff=0);\r
+\r
+  // method to extend AliExternalTrackParam functionality\r
+  Bool_t GetXatLabR(AliExternalTrackParam* tr,Double_t r,Double_t &x, Double_t bz, Int_t dir=0) const;\r
+\r
+  Double_t* PrepareEffFakeKombinations(TMatrixD *probKomb, TMatrixD *probLay);\r
+\r
+ protected:\r
\r
+  Int_t fNumberOfLayers;        // total number of layers in the model\r
+  Int_t fNumberOfActiveLayers;  // number of active layers in the model\r
+  Int_t fNumberOfActiveITSLayers;  // number of active ITS layers in the model\r
+  TList fLayers;                // List of layer pointers\r
+  Float_t fBField;              // Magnetic Field in Tesla\r
+  Float_t fLhcUPCscale;         // UltraPeripheralElectrons: scale from RHIC to LHC\r
+  Float_t fIntegrationTime;     // electronics integration time\r
+  Float_t fConfLevel;           // Confidence Level for the tracking\r
+  Float_t fAvgRapidity;         // rapidity of the track (= mean)\r
+  Float_t fParticleMass;        // Particle used for tracking. Standard: mass of pion\r
+  Double_t fMaxRadiusSlowDet;   // Maximum radius for slow detectors.  Fast detectors \r
+                                // and only fast detectors reside outside this radius.\r
+  Int_t fAtLeastCorr;     // min. number of correct hits for the track to be "good"\r
+  Int_t fAtLeastFake;     // min. number of fake hits for the track to be "fake"\r
+\r
+  Int_t fdNdEtaCent;       // Multiplicity\r
+\r
+  enum {kMaxNumberOfDetectors = 200};\r
\r
+  Double_t fTransMomenta[kNptBins];                          // array of transverse momenta\r
+  Double_t fMomentumRes[kNptBins];                           // array of momentum resolution\r
+  Double_t fResolutionRPhi[kNptBins];                        // array of rphi resolution\r
+  Double_t fResolutionZ[kNptBins];                           // array of z resolution\r
+  Double_t fDetPointRes[kMaxNumberOfDetectors][kNptBins];    // array of rphi resolution per layer\r
+  Double_t fDetPointZRes[kMaxNumberOfDetectors][kNptBins];   // array of z resolution per layer\r
+  Double_t fEfficiency[3][kNptBins];                         // efficiency for different particles\r
+  Double_t fFake[3][kNptBins];                               // fake prob for different particles\r
+\r
+  ClassDef(DetectorK,1);\r
+};\r
+\r
+#endif\r
diff --git a/ITS/UPGRADE/macros/FastVsSlowSim.C b/ITS/UPGRADE/macros/FastVsSlowSim.C
new file mode 100644 (file)
index 0000000..fc0cad1
--- /dev/null
@@ -0,0 +1,1367 @@
+#ifndef __CINT__
+#include <TTree.h>
+#include <TFile.h>
+#include <TClonesArray.h>
+#include <TStyle.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TProfile.h>
+#include <TGraph.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
+#include "AliRawReader.h"
+#include <AliTPCRawStreamV3.h>
+#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 <TROOT.h>
+#include <TStyle.h>
+#include "TGraphErrors.h"
+#include <TF1.h>
+#include <TGeoManager.h>
+
+#endif
+#include <iostream>
+#include <iomanip>
+#include <cstdio>
+
+/*
+ .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; i<nlayers; i++) {
+    rmin->AddAt(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; i<rmin.GetSize(); i++) {
+    cout<<i<<": (rmin,rmax)=("<<rmin.At(i)<<","<<rmax.At(i)<<")cm"<<endl;
+  }
+}
+
+
+void CountTrackableMCs(TH1F *hAllMC=0, Bool_t onlyPrims=0,Bool_t onlyPion=0);
+void CountPrimaries(TH1F *hMultCount);
+
+
+void ExtractOutputHistos(Bool_t onlyPrims=0,Bool_t onlyPion=0,Int_t plotFlag=0) {
+
+  //  gROOT->SetStyle("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<<xbins[i]<<endl;
+  }
+  //  TH1F *h = new TH1F("h","hist with log x axis",nbins,xbins);
+
+  TH1F *hMultCount = new TH1F("mult","averaged multiplicity (charg. prim)",80,-4.,4.);
+  hMultCount->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; iEv<tree->GetEntries(); 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 "<<label<<"<0) !! -> 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 "<<ncl<<endl;
+
+         nTrackFound++;
+         hAllFound->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"<<endl;
+      }
+      
+    } //entries track esd
+  
+
+  }//entries tree
+  runLoader->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; iEv<runLoader->GetNumberOfEvents(); 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; iMC<stack->GetNtrack(); 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; iref<nref; iref++){
+       AliTrackReference *trR = (AliTrackReference*)trackRef->At(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.GetSize();il++) {
+         if (radPos>=rmin.At(il)-0.1 && radPos<=rmax.At(il)+0.1) {
+           //      cout<<"  in Layer "<<il<<" "<<radPos;
+           nLaySigs.AddAt(1.,il);
+           //      cout<<" "<<nLaySigs.At(il)<<endl;
+         }
+       }
+      }
+
+      if (nLaySigs.GetSum()>=3) {
+       isTrackable =1;
+       //      cout<<nLaySigs.GetSum()<<endl;
+      }
+      
+      if (isTrackable) {
+       Double_t ptMC = particle->Pt();
+       //      Double_t etaMC = particle->Eta();
+       //      if (ptMC>ptmin&&ptMC<ptmax) {nTrackableTracks++;hAllMC->Fill(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 "<<nEvents<<endl;
+  for(Int_t iEv=0; iEv<nEvents; iEv++){
+    rl->GetEvent(iEv);
+    AliStack *s = rl->Stack();
+    for(Int_t iP=0; iP<s->GetNtrack(); 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();
+
+
+
+}
+