Another histos for lumi
[u/mrichter/AliRoot.git] / ITS / AliITSMeanVertexer.cxx
index a346ee9..caa5e59 100644 (file)
@@ -21,7 +21,6 @@
 
 const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
 ClassImp(AliITSMeanVertexer)
-
 ///////////////////////////////////////////////////////////////////////
 //                                                                   //
 // Class to compute vertex position using SPD local reconstruction   //
@@ -32,9 +31,7 @@ ClassImp(AliITSMeanVertexer)
 // Usage:                                                            //
 // Class used by the ITSSPDVertexDiamondda.cxx detector algorithm    //
 ///////////////////////////////////////////////////////////////////////
-
 /* $Id$ */
-
 //______________________________________________________________________
 AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
 fDetTypeRec(NULL),
@@ -52,20 +49,24 @@ fClu0(NULL),
 fIndex(0),
 fErrXCut(0.),
 fRCut(0.),
+fZCutDiamond(0.),
 fLowSPD0(0),
 fHighSPD0(0),
 fMultH(NULL),
 fErrXH(NULL),
 fMultHa(NULL),
 fErrXHa(NULL),
-fDistH(NULL)
+fDistH(NULL),
+fContrH(NULL),
+fContrHa(NULL)
 {
   // Default Constructor
-  Reset();
+  SetZFiducialRegion();   //  sets fZCutDiamond to its default value
+  Reset(kFALSE,kFALSE);
 
   // Histograms initialization
-  const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
-  const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
+  const Float_t xLimit = 0.6, yLimit = 0.6, zLimit = 50.0;
+  const Float_t xDelta = 0.003, yDelta = 0.003, zDelta = 0.2;
   fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
                       2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
                       2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
@@ -75,21 +76,64 @@ fDistH(NULL)
   fVertexZ  = new TH1F("VertexZ"," Longitudinal Vertex Profile",
                       2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
   fVertexZ->SetXTitle("Z , cm");
-  fErrXH = new TH1F("errorX","Error X - before cuts",50,0.,49.);
-  fMultH = new TH1F("multh","mult on layer 1 before cuts",100,0.,300.);
-  fErrXHa = new TH1F("errorXa","Error X - after Filter",50,0.,49.);
+  fErrXH = new TH1F("errorX","Error X - before cuts",100,0.,99.);
+  fMultH = new TH1F("multh","mult on layer 1 before cuts",2400,0.,7200.);
+  fErrXHa = new TH1F("errorXa","Error X - after Filter",100,0.,99.);
   fErrXHa->SetLineColor(kRed);
-  fMultHa = new TH1F("multha","mult on layer 1 - after Filter",100,0.,300.);
+  fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
   fMultHa->SetLineColor(kRed);
   fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
+  fContrH = new TH1F("contrib","Number of contributors - before cuts",300,0.5,300.5);
+  fContrHa = new TH1F("contriba","Number of contributors - after cuts",300,0.5,300.5);
+  fContrHa->SetLineColor(kRed);
+
   fClu0 = new UInt_t [fgkMaxNumOfEvents];
   for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
   fAccEvents.ResetAllBits(kTRUE);
 }
 
 //______________________________________________________________________
-void AliITSMeanVertexer::Reset(){
+void AliITSMeanVertexer::Reset(Bool_t redefine2D,Bool_t complete){
   // Reset averages 
+  // Both arguments are kFALSE when the method is called by the constructor
+  // redefine2D is TRUE when the method is called by ComputeMean at the second
+  //            pass.
+  // redefine2D is FALSE and complete is TRUE when the method is called
+  //            after a failure cof ComputeMean(kTRUE)
+
+  static Int_t counter = 0;
+  if(fVertexXY && redefine2D){
+    counter++;
+    Double_t rangeX=fVertexXY->GetXaxis()->GetXmax()-fVertexXY->GetXaxis()->GetXmin();
+    Double_t rangeY=fVertexXY->GetYaxis()->GetXmax()-fVertexXY->GetYaxis()->GetXmin();
+    Int_t nx=fVertexXY->GetNbinsX();
+    Int_t ny=fVertexXY->GetNbinsY();
+    delete fVertexXY;
+    Double_t xmi=fWeighPos[0]-rangeX/2.;
+    Double_t xma=fWeighPos[0]+rangeX/2.;
+    Double_t ymi=fWeighPos[1]-rangeY/2.;
+    Double_t yma=fWeighPos[1]+rangeY/2.;
+    fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",nx,xmi,xma,ny,ymi,yma);
+    fVertexXY->SetXTitle("X , cm");
+    fVertexXY->SetYTitle("Y , cm");
+    fVertexXY->SetOption("colz");
+    fVertexXY->SetDirectory(0);
+  }
+  else if(fVertexXY && !redefine2D){
+    fVertexXY->Reset();
+  }
+  if(fVertexZ){
+    fVertexZ->Reset();
+    fDistH->Reset();
+    if(complete){
+      fErrXH->Reset();
+      fMultH->Reset();
+      fErrXHa->Reset(); 
+      fMultHa->Reset();
+      fContrH->Reset();
+      fContrHa->Reset();
+    }
+  }
   for(Int_t i=0;i<3;i++){
     fWeighPosSum[i] = 0.;
     fWeighSigSum[i] = 0.;
@@ -101,8 +145,6 @@ void AliITSMeanVertexer::Reset(){
     for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
     fNoEventsContr=0;
     fTotContributors = 0.;
-    if(fVertexXY)fVertexXY->Reset();
-    if(fVertexZ)fVertexZ->Reset();
   }
 }
 //______________________________________________________________________
@@ -139,7 +181,7 @@ Bool_t AliITSMeanVertexer::Init() {
     fVertexer = new AliITSVertexer3D();
     fVertexer->SetDetTypeRec(fDetTypeRec);
     AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
-    alias->SetWideFiducialRegion(40.,0.5);
+    alias->SetWideFiducialRegion(fZCutDiamond,0.5);
     alias->SetNarrowFiducialRegion(0.5,0.5);
     alias->SetDeltaPhiCuts(0.5,0.025);
     alias->SetDCACut(0.1);
@@ -160,6 +202,8 @@ AliITSMeanVertexer::~AliITSMeanVertexer() {
   delete fMultHa;
   delete fErrXHa;
   delete fDistH;
+  delete fContrH;
+  delete fContrHa;
   delete fVertexer;
 }
 
@@ -176,7 +220,7 @@ Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
     AliInfo(Form("Low and high cuts on SPD L0 clusters %d , %d \n",fLowSPD0,fHighSPD0));
     AliInfo(Form("Reconstruct: cut on errX %f \n",fErrXCut));
   }
-  //  AliCodeTimerAuto("",0);
+//  AliCodeTimerAuto("",0);
   AliITSRecPointContainer::Instance()->PrepareToRead();
   TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
   fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
@@ -238,11 +282,16 @@ void AliITSMeanVertexer::WriteVertices(const char *filename){
       cov[3] =  fAverPosSq[0][2];  // cov xz
       cov[4] =  fAverPosSq[1][2];  // cov yz
       cov[5] =  fAverPosSq[2][2];  // variance z
-      AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
+      // We use standard average and not weighed averag now
+      // AliMeanVertex is apparently not taken by the preprocessor; only
+      // the AliESDVertex object is retrieved
+      //      AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
+      AliMeanVertex mv(fAverPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
       mv.SetTitle("Mean Vertex");
       mv.SetName("MeanVertex");
       // we have to add chi2 here
-      AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
+      //      AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
+      AliESDVertex vtx(fAverPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
 
       mv.Write(mv.GetName(),TObject::kOverwrite);
       vtx.Write(vtx.GetName(),TObject::kOverwrite);
@@ -261,6 +310,8 @@ void AliITSMeanVertexer::WriteVertices(const char *filename){
   fMultHa->Write(fMultHa->GetName(),TObject::kOverwrite);
   fErrXHa->Write(fErrXHa->GetName(),TObject::kOverwrite);
   fDistH->Write(fDistH->GetName(),TObject::kOverwrite);
+  fContrH->Write(fContrH->GetName(),TObject::kOverwrite);
+  fContrHa->Write(fContrHa->GetName(),TObject::kOverwrite);
   fmv.Close();
 }
 
@@ -276,11 +327,13 @@ Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
   Double_t ex = vert->GetXRes();
   fMultH->Fill(mult);
   fErrXH->Fill(ex*1000.);
+  fContrH->Fill((Float_t)ncontr);
   if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) && 
      ((ex*1000.)<fErrXCut)) {
     status = kTRUE;
     fMultHa->Fill(mult);
     fErrXHa->Fill(ex*1000.);
+    fContrHa->Fill((Float_t)ncontr);
   }
   return status;
 }
@@ -318,18 +371,21 @@ void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
    // called once with killOutliers = kFALSE and then again with
    // killOutliers = kTRUE
 
-   static Bool_t firstime = kTRUE;
-   if(fNoEventsContr>2 && !killOutliers)return kTRUE;  //ComputeMean is
+   static Bool_t preliminary = kTRUE;
+   static Int_t oldnumbevt = 0;
+   if(!(preliminary || killOutliers))return kTRUE;  //ComputeMean is
                              // executed only once with argument kFALSE
    Double_t wpos[3];
    for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
 
-   if(killOutliers && firstime){
-     firstime = kFALSE;
-     Reset();
+   Int_t istart = oldnumbevt;
+   if(killOutliers)istart = 0;
+   if(killOutliers && preliminary){
+     preliminary = kFALSE;
+     Reset(kTRUE,kFALSE);
    }
-
-   for(Int_t i=0;i<fVertArray.GetEntries();i++){
+   oldnumbevt = fVertArray.GetEntries();
+   for(Int_t i=istart;i<fVertArray.GetEntries();i++){
      AliESDVertex* vert = NULL;
 
      // second pass
@@ -346,7 +402,19 @@ void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
      vert=(AliESDVertex*)fVertArray[i];
      AddToMean(vert);
    }
-   if(fNoEventsContr < 2) return kFALSE;
+   Bool_t bad = ((!killOutliers) && fNoEventsContr < 5) || (killOutliers && fNoEventsContr <2);
+   if(bad) {
+     if(killOutliers){  
+// when the control reachs this point, the preliminary evaluation of the
+// diamond region has to be redone. So a general reset must be issued
+// if this occurs, it is likely that the cuts are badly defined
+       ResetArray();
+       Reset(kFALSE,kTRUE);
+       preliminary = kTRUE;
+       oldnumbevt = 0;
+     }
+     return kFALSE;
+   }
    Double_t nevents = fNoEventsContr;
    for(Int_t i=0;i<3;i++){
      fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];