[u/mrichter/AliRoot.git] / PWGPP / ITS / AliMeanVertexCalibTask.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//*************************************************************************
// Class AliMeanVertexCalibTask
// AliAnalysisTask to extract from ESD the information on primary vertex
// reconstruction in order to compute the MeanVertex object
//
// Author:  D.Caffarri, davide.caffarri@pd.infn.it  
//          A.Dainese, andrea.dainese@pd.infn.it
//*************************************************************************


#include <TH1F.h>
#include <TH2F.h>

#include "AliMultiplicity.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliVertexerTracks.h"

#include "AliMeanVertexCalibTask.h"


ClassImp(AliMeanVertexCalibTask)

//_____________________________________________________________________
AliMeanVertexCalibTask::AliMeanVertexCalibTask(const char *name):
  AliAnalysisTaskSE(name),
  fESD(0), 
  fOutput(0),
  fOnlyITSTPCTracks(kFALSE),
  fOnlyITSSATracks(kTRUE)
{

  // Constructor
  
  // Define input and output slots here
  // Output slot #0 writes into a TList container
  DefineOutput(1, TList::Class());  //My private output
}

//________________________________________________________________________
AliMeanVertexCalibTask::~AliMeanVertexCalibTask()
{
  // Destructor

  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
  
}
//________________________________________________________________________
void AliMeanVertexCalibTask::UserCreateOutputObjects()
{

  // Several histograms are more conveniently managed in a TList
  fOutput = new TList;
  fOutput->SetOwner();

  
  TH1F* hSPDVertexX = new TH1F("hSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
  fOutput->Add(hSPDVertexX);
  TH1F* hSPDVertexY = new TH1F("hSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
  fOutput->Add(hSPDVertexY);
  TH1F* hSPDVertexZ = new TH1F("hSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
  fOutput->Add(hSPDVertexZ);
  TH1F* hTRKVertexX = new TH1F("hTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
  fOutput->Add(hTRKVertexX);
  TH1F* hTRKVertexY = new TH1F("hTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
  fOutput->Add(hTRKVertexY);
  TH1F* hTRKVertexZ = new TH1F("hTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
  fOutput->Add(hTRKVertexZ);
  
  TH2F *hTRKVertexXZ = new TH2F("hTRKVertexXZ", "TRKVertex XZ corr", 200, -1, 1, 200, -20, 20);
  fOutput->Add(hTRKVertexXZ);

  TH2F *hTRKVertexYZ = new TH2F("hTRKVertexYZ", "TRKVertex YZ corr", 200, -1, 1, 200, -20, 20);
  fOutput->Add(hTRKVertexYZ);
  
  TH1F* hTRKVertexXdefMult = new TH1F("hTRKVertexXdefMult","TRKVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);  
  fOutput->Add(hTRKVertexXdefMult);
  TH1F* hTRKVertexYdefMult = new TH1F("hTRKVertexYdefMult","TRKVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
  fOutput->Add(hTRKVertexYdefMult);
	
	TH1F* hTRKVertexXHighMult = new TH1F("hTRKVertexXHighMult","TRKVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);  
	fOutput->Add(hTRKVertexXHighMult);
	TH1F* hTRKVertexYHighMult = new TH1F("hTRKVertexYHighMult","TRKVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
	fOutput->Add(hTRKVertexYHighMult);	
	

  TH1F* hITSSAVertexX = new TH1F("hITSSAVertexX","ITSSAVertex x; x vertex [cm]; events",200,-1,1);
  fOutput->Add(hITSSAVertexX);
  TH1F* hITSSAVertexY = new TH1F("hITSSAVertexY","ITSSAVertex y; y vertex [cm]; events",200,-1,1);
  fOutput->Add(hITSSAVertexY);
  TH1F* hITSSAVertexZ = new TH1F("hITSSAVertexZ","ITSSAVertex z; z vertex [cm]; events",200,-20,20);
  fOutput->Add(hITSSAVertexZ);
  
  TH2F *hITSSAVertexXZ = new TH2F("hITSSAVertexXZ", "ITSSAVertex XZ corr", 200, -1, 1, 200, -20, 20);
  fOutput->Add(hITSSAVertexXZ);

  TH2F *hITSSAVertexYZ = new TH2F("hITSSAVertexYZ", "ITSSAVertex YZ corr", 200, -1, 1, 200, -20, 20);
  fOutput->Add(hITSSAVertexYZ);
  
  TH1F* hITSSAVertexXdefMult = new TH1F("hITSSAVertexXdefMult","ITSSAVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);  
  fOutput->Add(hITSSAVertexXdefMult);
  TH1F* hITSSAVertexYdefMult = new TH1F("hITSSAVertexYdefMult","ITSSAVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
  fOutput->Add(hITSSAVertexYdefMult);

	
   TH1F* hITSSAVertexXHighMult = new TH1F("hITSSAVertexXHighMult","ITSSAVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);  
   fOutput->Add(hITSSAVertexXHighMult);
   TH1F* hITSSAVertexYHighMult = new TH1F("hITSSAVertexYHighMult","ITSSAVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
   fOutput->Add(hITSSAVertexYHighMult);
	
	
  PostData(1, fOutput);
  
  return;
}


//________________________________________________________________________
void AliMeanVertexCalibTask::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
  
  if (!InputEvent()) {
    Printf("ERROR: fESD not available");
    return;
  }
  
  AliESDEvent* esdE = (AliESDEvent*) InputEvent();

  const AliMultiplicity *alimult = esdE->GetMultiplicity();
  Int_t ntrklets=0;
  if(alimult) ntrklets = alimult->GetNumberOfTracklets();
   

  const AliESDVertex *trkv = esdE->GetPrimaryVertexTracks();
  const AliESDVertex *spdv = esdE->GetPrimaryVertexSPD();

  AliESDVertex *itsSAv = ReconstructPrimaryVertex(kFALSE); 

  if(spdv) {
    if(spdv->GetNContributors()>0) {
      TString title=spdv->GetTitle();
      if(title.Contains("3D")) {
	((TH1F*)fOutput->FindObject("hSPDVertexX"))->Fill(spdv->GetXv());
	((TH1F*)fOutput->FindObject("hSPDVertexY"))->Fill(spdv->GetYv());
      }
      ((TH1F*)fOutput->FindObject("hSPDVertexZ"))->Fill(spdv->GetZv());
    }
  }
  
  
  if(trkv) {
    if(trkv->GetNContributors()>0) {
      ((TH1F*)fOutput->FindObject("hTRKVertexX"))->Fill(trkv->GetXv());
      ((TH1F*)fOutput->FindObject("hTRKVertexY"))->Fill(trkv->GetYv());
      ((TH1F*)fOutput->FindObject("hTRKVertexZ"))->Fill(trkv->GetZv());
      
      if (ntrklets>30 && ntrklets<45){
	((TH1F*)fOutput->FindObject("hTRKVertexXdefMult"))->Fill(trkv->GetXv());
	((TH1F*)fOutput->FindObject("hTRKVertexYdefMult"))->Fill(trkv->GetYv());
      }
		
	  if (ntrklets>1500){
			((TH1F*)fOutput->FindObject("hTRKVertexXHighMult"))->Fill(trkv->GetXv());
			((TH1F*)fOutput->FindObject("hTRKVertexYHighMult"))->Fill(trkv->GetYv());
		}
      
      ((TH2F*)fOutput->FindObject("hTRKVertexXZ"))->Fill(trkv->GetXv(),trkv->GetZv());
      ((TH2F*)fOutput->FindObject("hTRKVertexYZ"))->Fill(trkv->GetYv(),trkv->GetZv());
      
    }
  }

  if (itsSAv){
    if (itsSAv->GetNContributors()>0){
  
      ((TH1F*)fOutput->FindObject("hITSSAVertexX"))->Fill(itsSAv->GetXv());
      ((TH1F*)fOutput->FindObject("hITSSAVertexY"))->Fill(itsSAv->GetYv());
      ((TH1F*)fOutput->FindObject("hITSSAVertexZ"))->Fill(itsSAv->GetZv());
      
		if (ntrklets>30 && ntrklets<45){
			((TH1F*)fOutput->FindObject("hITSSAVertexXdefMult"))->Fill(itsSAv->GetXv());
		    ((TH1F*)fOutput->FindObject("hITSSAVertexYdefMult"))->Fill(itsSAv->GetYv());
		}
		
		if (ntrklets>1500){
			((TH1F*)fOutput->FindObject("hITSSAVertexXHighMult"))->Fill(itsSAv->GetXv());
			((TH1F*)fOutput->FindObject("hITSSAVertexYHighMult"))->Fill(itsSAv->GetYv());
		}
      
      ((TH2F*)fOutput->FindObject("hITSSAVertexXZ"))->Fill(itsSAv->GetXv(),itsSAv->GetZv());
      ((TH2F*)fOutput->FindObject("hITSSAVertexYZ"))->Fill(itsSAv->GetYv(),itsSAv->GetZv());
     
    }
  }

  delete itsSAv;
  
  PostData(1, fOutput);
  
  return;
  
}


//________________________________________________________________________
void AliMeanVertexCalibTask::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query
  fOutput = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutput) {     
    Printf("ERROR: fOutput not available");
    return;
  }


  return;

}


//__________________________________________________________________________
AliESDVertex* AliMeanVertexCalibTask::ReconstructPrimaryVertex(Bool_t constr,Int_t mode) const {
  // On the fly reco of ITS+TPC vertex from ESD
  // mode=0 use all tracks
  // mode=1 use odd-number tracks
  // mode=2 use even-number tracks

  AliESDEvent* evt = (AliESDEvent*) fInputEvent;
  AliVertexerTracks vertexer(evt->GetMagneticField());
  if(evt->GetNumberOfTracks()<500) {
    vertexer.SetITSMode(); // defaults
    vertexer.SetMinClusters(4); // default is 5
  } else { 
    vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb
  } 
  if (fOnlyITSSATracks) vertexer.SetITSpureSA();
  
  Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
  Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0}; 
  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
  AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
  vertexer.SetVtxStart(initVertex);
  delete initVertex;
  if(!constr) vertexer.SetConstraintOff();

  if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>0) {
    Int_t iskip=0;
    Int_t ntracks = evt->GetNumberOfTracks();
    Int_t *skip = new Int_t[ntracks];
    for(Int_t i=0;i<ntracks;i++) skip[i]=-1;
    for(Int_t itr=0;itr<ntracks; itr++) {
      AliESDtrack* track = evt->GetTrack(itr);
      if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
	skip[iskip++]=itr;
	continue;
      }
      if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
	skip[iskip++]=itr;
	continue;
      }
      if(mode==1 && itr%2==0) skip[iskip++]=itr;
      if(mode==2 && itr%2==1) skip[iskip++]=itr;
    }
    vertexer.SetSkipTracks(iskip,skip);
    delete [] skip; skip=NULL;
  }
  
  return vertexer.FindPrimaryVertex(evt);
}
Commit	Line	Data
bbb57828	1	/**************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//*************************************************************************
	17	// Class AliMeanVertexCalibTask
	18	// AliAnalysisTask to extract from ESD the information on primary vertex
	19	// reconstruction in order to compute the MeanVertex object
	20	//
	21	// Author: D.Caffarri, davide.caffarri@pd.infn.it
	22	// A.Dainese, andrea.dainese@pd.infn.it
	23	//*************************************************************************
	24
	25
	26	#include <TH1F.h>
	27	#include <TH2F.h>
	28
	29	#include "AliMultiplicity.h"
	30	#include "AliESDVertex.h"
	31	#include "AliESDEvent.h"
	32	#include "AliVertexerTracks.h"
	33
	34	#include "AliMeanVertexCalibTask.h"
	35
	36
	37	ClassImp(AliMeanVertexCalibTask)
	38
	39	//_____________________________________________________________________
	40	AliMeanVertexCalibTask::AliMeanVertexCalibTask(const char *name):
	41	AliAnalysisTaskSE(name),
	42	fESD(0),
	43	fOutput(0),
	44	fOnlyITSTPCTracks(kFALSE),
	45	fOnlyITSSATracks(kTRUE)
	46	{
	47
	48	// Constructor
	49
	50	// Define input and output slots here
	51	// Output slot #0 writes into a TList container
	52	DefineOutput(1, TList::Class()); //My private output
	53	}
	54
	55	//________________________________________________________________________
	56	AliMeanVertexCalibTask::~AliMeanVertexCalibTask()
	57	{
	58	// Destructor
	59
	60	if (fOutput) {
	61	delete fOutput;
	62	fOutput = 0;
	63	}
	64
65	}
66	//________________________________________________________________________
67	void AliMeanVertexCalibTask::UserCreateOutputObjects()
68	{
69
70	// Several histograms are more conveniently managed in a TList
71	fOutput = new TList;
72	fOutput->SetOwner();
73
74
75	TH1F* hSPDVertexX = new TH1F("hSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
76	fOutput->Add(hSPDVertexX);
77	TH1F* hSPDVertexY = new TH1F("hSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
78	fOutput->Add(hSPDVertexY);
79	TH1F* hSPDVertexZ = new TH1F("hSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
80	fOutput->Add(hSPDVertexZ);
81	TH1F* hTRKVertexX = new TH1F("hTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
82	fOutput->Add(hTRKVertexX);
83	TH1F* hTRKVertexY = new TH1F("hTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
84	fOutput->Add(hTRKVertexY);
85	TH1F* hTRKVertexZ = new TH1F("hTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
86	fOutput->Add(hTRKVertexZ);
87
88	TH2F *hTRKVertexXZ = new TH2F("hTRKVertexXZ", "TRKVertex XZ corr", 200, -1, 1, 200, -20, 20);
89	fOutput->Add(hTRKVertexXZ);
90
91	TH2F *hTRKVertexYZ = new TH2F("hTRKVertexYZ", "TRKVertex YZ corr", 200, -1, 1, 200, -20, 20);
92	fOutput->Add(hTRKVertexYZ);
93
94	TH1F* hTRKVertexXdefMult = new TH1F("hTRKVertexXdefMult","TRKVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
95	fOutput->Add(hTRKVertexXdefMult);
96	TH1F* hTRKVertexYdefMult = new TH1F("hTRKVertexYdefMult","TRKVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
97	fOutput->Add(hTRKVertexYdefMult);
a680415a	98
	99	TH1F* hTRKVertexXHighMult = new TH1F("hTRKVertexXHighMult","TRKVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
	100	fOutput->Add(hTRKVertexXHighMult);
	101	TH1F* hTRKVertexYHighMult = new TH1F("hTRKVertexYHighMult","TRKVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
	102	fOutput->Add(hTRKVertexYHighMult);
	103
bbb57828	104
	105	TH1F* hITSSAVertexX = new TH1F("hITSSAVertexX","ITSSAVertex x; x vertex [cm]; events",200,-1,1);
	106	fOutput->Add(hITSSAVertexX);
	107	TH1F* hITSSAVertexY = new TH1F("hITSSAVertexY","ITSSAVertex y; y vertex [cm]; events",200,-1,1);
	108	fOutput->Add(hITSSAVertexY);
	109	TH1F* hITSSAVertexZ = new TH1F("hITSSAVertexZ","ITSSAVertex z; z vertex [cm]; events",200,-20,20);
	110	fOutput->Add(hITSSAVertexZ);
	111
	112	TH2F *hITSSAVertexXZ = new TH2F("hITSSAVertexXZ", "ITSSAVertex XZ corr", 200, -1, 1, 200, -20, 20);
	113	fOutput->Add(hITSSAVertexXZ);
	114
	115	TH2F *hITSSAVertexYZ = new TH2F("hITSSAVertexYZ", "ITSSAVertex YZ corr", 200, -1, 1, 200, -20, 20);
	116	fOutput->Add(hITSSAVertexYZ);
	117
	118	TH1F* hITSSAVertexXdefMult = new TH1F("hITSSAVertexXdefMult","ITSSAVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
	119	fOutput->Add(hITSSAVertexXdefMult);
	120	TH1F* hITSSAVertexYdefMult = new TH1F("hITSSAVertexYdefMult","ITSSAVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
	121	fOutput->Add(hITSSAVertexYdefMult);
	122
a680415a	123
	124	TH1F* hITSSAVertexXHighMult = new TH1F("hITSSAVertexXHighMult","ITSSAVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
	125	fOutput->Add(hITSSAVertexXHighMult);
	126	TH1F* hITSSAVertexYHighMult = new TH1F("hITSSAVertexYHighMult","ITSSAVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
	127	fOutput->Add(hITSSAVertexYHighMult);
	128
	129
bbb57828	130
	131	PostData(1, fOutput);
	132
	133	return;
	134	}
	135
	136
	137	//________________________________________________________________________
	138	void AliMeanVertexCalibTask::UserExec(Option_t *)
	139	{
	140	// Main loop
	141	// Called for each event
	142
	143	if (!InputEvent()) {
	144	Printf("ERROR: fESD not available");
	145	return;
	146	}
	147
	148	AliESDEvent* esdE = (AliESDEvent*) InputEvent();
	149
bbb57828	150	const AliMultiplicity *alimult = esdE->GetMultiplicity();
	151	Int_t ntrklets=0;
	152	if(alimult) ntrklets = alimult->GetNumberOfTracklets();
	153
	154
	155	const AliESDVertex *trkv = esdE->GetPrimaryVertexTracks();
	156	const AliESDVertex *spdv = esdE->GetPrimaryVertexSPD();
	157
	158	AliESDVertex *itsSAv = ReconstructPrimaryVertex(kFALSE);
	159
	160	if(spdv) {
	161	if(spdv->GetNContributors()>0) {
	162	TString title=spdv->GetTitle();
	163	if(title.Contains("3D")) {
	164	((TH1F*)fOutput->FindObject("hSPDVertexX"))->Fill(spdv->GetXv());
	165	((TH1F*)fOutput->FindObject("hSPDVertexY"))->Fill(spdv->GetYv());
	166	}
	167	((TH1F*)fOutput->FindObject("hSPDVertexZ"))->Fill(spdv->GetZv());
	168	}
	169	}
	170
	171
	172	if(trkv) {
	173	if(trkv->GetNContributors()>0) {
	174	((TH1F*)fOutput->FindObject("hTRKVertexX"))->Fill(trkv->GetXv());
	175	((TH1F*)fOutput->FindObject("hTRKVertexY"))->Fill(trkv->GetYv());
	176	((TH1F*)fOutput->FindObject("hTRKVertexZ"))->Fill(trkv->GetZv());
	177
	178	if (ntrklets>30 && ntrklets<45){
	179	((TH1F*)fOutput->FindObject("hTRKVertexXdefMult"))->Fill(trkv->GetXv());
	180	((TH1F*)fOutput->FindObject("hTRKVertexYdefMult"))->Fill(trkv->GetYv());
	181	}
a680415a	182
	183	if (ntrklets>1500){
	184	((TH1F*)fOutput->FindObject("hTRKVertexXHighMult"))->Fill(trkv->GetXv());
	185	((TH1F*)fOutput->FindObject("hTRKVertexYHighMult"))->Fill(trkv->GetYv());
	186	}
bbb57828	187
	188	((TH2F*)fOutput->FindObject("hTRKVertexXZ"))->Fill(trkv->GetXv(),trkv->GetZv());
	189	((TH2F*)fOutput->FindObject("hTRKVertexYZ"))->Fill(trkv->GetYv(),trkv->GetZv());
	190
	191	}
	192	}
	193
	194	if (itsSAv){
	195	if (itsSAv->GetNContributors()>0){
	196
	197	((TH1F*)fOutput->FindObject("hITSSAVertexX"))->Fill(itsSAv->GetXv());
	198	((TH1F*)fOutput->FindObject("hITSSAVertexY"))->Fill(itsSAv->GetYv());
	199	((TH1F*)fOutput->FindObject("hITSSAVertexZ"))->Fill(itsSAv->GetZv());
	200
a680415a	201	if (ntrklets>30 && ntrklets<45){
	202	((TH1F*)fOutput->FindObject("hITSSAVertexXdefMult"))->Fill(itsSAv->GetXv());
	203	((TH1F*)fOutput->FindObject("hITSSAVertexYdefMult"))->Fill(itsSAv->GetYv());
	204	}
	205
	206	if (ntrklets>1500){
	207	((TH1F*)fOutput->FindObject("hITSSAVertexXHighMult"))->Fill(itsSAv->GetXv());
	208	((TH1F*)fOutput->FindObject("hITSSAVertexYHighMult"))->Fill(itsSAv->GetYv());
	209	}
bbb57828	210
	211	((TH2F*)fOutput->FindObject("hITSSAVertexXZ"))->Fill(itsSAv->GetXv(),itsSAv->GetZv());
	212	((TH2F*)fOutput->FindObject("hITSSAVertexYZ"))->Fill(itsSAv->GetYv(),itsSAv->GetZv());
	213
	214	}
	215	}
	216
	217	delete itsSAv;
	218
	219	PostData(1, fOutput);
	220
	221	return;
	222
	223	}
	224
	225
	226	//________________________________________________________________________
	227	void AliMeanVertexCalibTask::Terminate(Option_t *)
	228	{
	229	// Draw result to the screen
	230	// Called once at the end of the query
	231	fOutput = dynamic_cast<TList*> (GetOutputData(1));
	232	if (!fOutput) {
	233	Printf("ERROR: fOutput not available");
	234	return;
	235	}
	236
	237
	238	return;
	239
	240	}
	241
	242
	243	//__________________________________________________________________________
	244	AliESDVertex* AliMeanVertexCalibTask::ReconstructPrimaryVertex(Bool_t constr,Int_t mode) const {
	245	// On the fly reco of ITS+TPC vertex from ESD
	246	// mode=0 use all tracks
	247	// mode=1 use odd-number tracks
	248	// mode=2 use even-number tracks
	249
	250	AliESDEvent* evt = (AliESDEvent*) fInputEvent;
	251	AliVertexerTracks vertexer(evt->GetMagneticField());
	252	if(evt->GetNumberOfTracks()<500) {
	253	vertexer.SetITSMode(); // defaults
	254	vertexer.SetMinClusters(4); // default is 5
	255	} else {
	256	vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb
	257	}
	258	if (fOnlyITSSATracks) vertexer.SetITSpureSA();
	259
	260	Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
	261	Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
	262	Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
	263	AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
	264	vertexer.SetVtxStart(initVertex);
	265	delete initVertex;
	266	if(!constr) vertexer.SetConstraintOff();
	267
	268	if(fOnlyITSTPCTracks \|\| fOnlyITSSATracks \|\| mode>0) {
	269	Int_t iskip=0;
550dc9e9	270	Int_t ntracks = evt->GetNumberOfTracks();
550dc9e9	271	Int_t *skip = new Int_t[ntracks];
22f3142a	272	for(Int_t i=0;i<ntracks;i++) skip[i]=-1;
550dc9e9	273	for(Int_t itr=0;itr<ntracks; itr++) {
bbb57828	274	AliESDtrack* track = evt->GetTrack(itr);
	275	if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
	276	skip[iskip++]=itr;
	277	continue;
	278	}
	279	if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
	280	skip[iskip++]=itr;
	281	continue;
	282	}
	283	if(mode==1 && itr%2==0) skip[iskip++]=itr;
	284	if(mode==2 && itr%2==1) skip[iskip++]=itr;
	285	}
	286	vertexer.SetSkipTracks(iskip,skip);
	287	delete [] skip; skip=NULL;
	288	}
	289
	290	return vertexer.FindPrimaryVertex(evt);
	291	}