[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 <string.h>

#include "AliMultiplicity.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliVertexerTracks.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "AliGRPRecoParam.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 char* beamType = esdE->GetBeamType();
  Printf("beam type = %s", beamType);

  Bool_t kLowFlux = kTRUE, kHighFlux = kFALSE;
  // TString pp= "p-p";
  //TString pA= "p-A";
  TString AA= "A-A";

  if (beamType == AA){ 
    kHighFlux = kTRUE;
    kLowFlux = kFALSE;
    Printf ("high flux setting");
    }
  
  AliCDBManager* man = AliCDBManager::Instance();
  man->SetDefaultStorage("raw://");
  Int_t runNb = esdE->GetRunNumber();
  man->SetRun(runNb);
  Printf("runNb = %d", runNb);
  AliCDBEntry *entry = (AliCDBEntry*)man->Get("GRP/Calib/RecoParam/");
  Printf("entry = %p", entry);

  TObjArray *arrayRecoParam = (TObjArray*)entry->GetObject();
  Printf("arrayRecoParam = %p", arrayRecoParam);
  
  AliGRPRecoParam *grpRecoParam=0x0;
  if (kLowFlux) grpRecoParam= (AliGRPRecoParam*)arrayRecoParam->At(1);
  else if (kHighFlux) grpRecoParam = (AliGRPRecoParam*)arrayRecoParam->At(2);
  
  grpRecoParam->Dump();

  AliVertexerTracks *vertexer= new AliVertexerTracks(esdE->GetMagneticField());
  vertexer->SetITSMode();
  vertexer->SetConstraintOff();
    
  if (grpRecoParam) {
    Int_t nCutsVertexer = grpRecoParam->GetVertexerTracksNCuts();
    Double_t *cutsVertexer = new Double_t[nCutsVertexer];
    grpRecoParam->GetVertexerTracksCutsITS(cutsVertexer,nCutsVertexer);
    vertexer->SetCuts(cutsVertexer,nCutsVertexer);
    delete [] cutsVertexer; cutsVertexer = NULL; 
   }
   
  vertexer->SetConstraintOff();

  //track vertex recomputed from the vertexer
  AliESDVertex *trkv = vertexer->FindPrimaryVertex(esdE);
  Printf ("trkv = %p", trkv); 
  //const AliESDVertex *trkv = esdE->GetPrimaryVertexTracks();
  
  //SPD vertex taken from the ESD 
  const AliESDVertex *spdv = esdE->GetPrimaryVertexSPD();

  //ITSSA vertex computed from the vertexer
  vertexer->SetITSpureSA();
  AliESDVertex *itsSAv = vertexer->FindPrimaryVertex(esdE);

  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;
  delete vertexer;
  
  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>
29be8a66	28	#include <string.h>
bbb57828	29
	30	#include "AliMultiplicity.h"
	31	#include "AliESDVertex.h"
	32	#include "AliESDEvent.h"
	33	#include "AliVertexerTracks.h"
29be8a66	34	#include "AliCDBManager.h"
	35	#include "AliCDBEntry.h"
	36	#include "AliGRPRecoParam.h"
bbb57828	37
	38	#include "AliMeanVertexCalibTask.h"
	39
	40
	41	ClassImp(AliMeanVertexCalibTask)
	42
	43	//_____________________________________________________________________
	44	AliMeanVertexCalibTask::AliMeanVertexCalibTask(const char *name):
	45	AliAnalysisTaskSE(name),
	46	fESD(0),
	47	fOutput(0),
	48	fOnlyITSTPCTracks(kFALSE),
	49	fOnlyITSSATracks(kTRUE)
	50	{
	51
	52	// Constructor
	53
	54	// Define input and output slots here
	55	// Output slot #0 writes into a TList container
	56	DefineOutput(1, TList::Class()); //My private output
	57	}
	58
	59	//________________________________________________________________________
	60	AliMeanVertexCalibTask::~AliMeanVertexCalibTask()
	61	{
	62	// Destructor
	63
	64	if (fOutput) {
	65	delete fOutput;
	66	fOutput = 0;
	67	}
	68
	69	}
	70	//________________________________________________________________________
	71	void AliMeanVertexCalibTask::UserCreateOutputObjects()
	72	{
	73
	74	// Several histograms are more conveniently managed in a TList
	75	fOutput = new TList;
	76	fOutput->SetOwner();
	77
	78
	79	TH1F* hSPDVertexX = new TH1F("hSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
	80	fOutput->Add(hSPDVertexX);
	81	TH1F* hSPDVertexY = new TH1F("hSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
	82	fOutput->Add(hSPDVertexY);
	83	TH1F* hSPDVertexZ = new TH1F("hSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
	84	fOutput->Add(hSPDVertexZ);
	85	TH1F* hTRKVertexX = new TH1F("hTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
	86	fOutput->Add(hTRKVertexX);
	87	TH1F* hTRKVertexY = new TH1F("hTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
	88	fOutput->Add(hTRKVertexY);
	89	TH1F* hTRKVertexZ = new TH1F("hTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
	90	fOutput->Add(hTRKVertexZ);
	91
	92	TH2F *hTRKVertexXZ = new TH2F("hTRKVertexXZ", "TRKVertex XZ corr", 200, -1, 1, 200, -20, 20);
	93	fOutput->Add(hTRKVertexXZ);
	94
	95	TH2F *hTRKVertexYZ = new TH2F("hTRKVertexYZ", "TRKVertex YZ corr", 200, -1, 1, 200, -20, 20);
	96	fOutput->Add(hTRKVertexYZ);
	97
	98	TH1F* hTRKVertexXdefMult = new TH1F("hTRKVertexXdefMult","TRKVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
	99	fOutput->Add(hTRKVertexXdefMult);
	100	TH1F* hTRKVertexYdefMult = new TH1F("hTRKVertexYdefMult","TRKVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
101	fOutput->Add(hTRKVertexYdefMult);
a680415a	102
	103	TH1F* hTRKVertexXHighMult = new TH1F("hTRKVertexXHighMult","TRKVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
	104	fOutput->Add(hTRKVertexXHighMult);
	105	TH1F* hTRKVertexYHighMult = new TH1F("hTRKVertexYHighMult","TRKVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
	106	fOutput->Add(hTRKVertexYHighMult);
	107
bbb57828	108
	109	TH1F* hITSSAVertexX = new TH1F("hITSSAVertexX","ITSSAVertex x; x vertex [cm]; events",200,-1,1);
	110	fOutput->Add(hITSSAVertexX);
	111	TH1F* hITSSAVertexY = new TH1F("hITSSAVertexY","ITSSAVertex y; y vertex [cm]; events",200,-1,1);
	112	fOutput->Add(hITSSAVertexY);
	113	TH1F* hITSSAVertexZ = new TH1F("hITSSAVertexZ","ITSSAVertex z; z vertex [cm]; events",200,-20,20);
	114	fOutput->Add(hITSSAVertexZ);
	115
	116	TH2F *hITSSAVertexXZ = new TH2F("hITSSAVertexXZ", "ITSSAVertex XZ corr", 200, -1, 1, 200, -20, 20);
	117	fOutput->Add(hITSSAVertexXZ);
	118
	119	TH2F *hITSSAVertexYZ = new TH2F("hITSSAVertexYZ", "ITSSAVertex YZ corr", 200, -1, 1, 200, -20, 20);
	120	fOutput->Add(hITSSAVertexYZ);
	121
	122	TH1F* hITSSAVertexXdefMult = new TH1F("hITSSAVertexXdefMult","ITSSAVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
	123	fOutput->Add(hITSSAVertexXdefMult);
	124	TH1F* hITSSAVertexYdefMult = new TH1F("hITSSAVertexYdefMult","ITSSAVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
	125	fOutput->Add(hITSSAVertexYdefMult);
	126
a680415a	127
	128	TH1F* hITSSAVertexXHighMult = new TH1F("hITSSAVertexXHighMult","ITSSAVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
	129	fOutput->Add(hITSSAVertexXHighMult);
	130	TH1F* hITSSAVertexYHighMult = new TH1F("hITSSAVertexYHighMult","ITSSAVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
	131	fOutput->Add(hITSSAVertexYHighMult);
	132
	133
bbb57828	134
	135	PostData(1, fOutput);
	136
	137	return;
	138	}
	139
	140
	141	//________________________________________________________________________
	142	void AliMeanVertexCalibTask::UserExec(Option_t *)
	143	{
	144	// Main loop
	145	// Called for each event
	146
	147	if (!InputEvent()) {
	148	Printf("ERROR: fESD not available");
	149	return;
	150	}
	151
	152	AliESDEvent* esdE = (AliESDEvent*) InputEvent();
	153
bbb57828	154	const AliMultiplicity *alimult = esdE->GetMultiplicity();
	155	Int_t ntrklets=0;
	156	if(alimult) ntrklets = alimult->GetNumberOfTracklets();
29be8a66	157
	158	const char* beamType = esdE->GetBeamType();
	159	Printf("beam type = %s", beamType);
	160
	161	Bool_t kLowFlux = kTRUE, kHighFlux = kFALSE;
	162	// TString pp= "p-p";
	163	//TString pA= "p-A";
	164	TString AA= "A-A";
	165
	166	if (beamType == AA){
	167	kHighFlux = kTRUE;
	168	kLowFlux = kFALSE;
	169	Printf ("high flux setting");
	170	}
	171
	172	AliCDBManager* man = AliCDBManager::Instance();
	173	man->SetDefaultStorage("raw://");
	174	Int_t runNb = esdE->GetRunNumber();
	175	man->SetRun(runNb);
	176	Printf("runNb = %d", runNb);
	177	AliCDBEntry entry = (AliCDBEntry)man->Get("GRP/Calib/RecoParam/");
	178	Printf("entry = %p", entry);
	179
	180	TObjArray arrayRecoParam = (TObjArray)entry->GetObject();
	181	Printf("arrayRecoParam = %p", arrayRecoParam);
	182
	183	AliGRPRecoParam *grpRecoParam=0x0;
	184	if (kLowFlux) grpRecoParam= (AliGRPRecoParam*)arrayRecoParam->At(1);
	185	else if (kHighFlux) grpRecoParam = (AliGRPRecoParam*)arrayRecoParam->At(2);
	186
	187	grpRecoParam->Dump();
	188
	189	AliVertexerTracks *vertexer= new AliVertexerTracks(esdE->GetMagneticField());
	190	vertexer->SetITSMode();
	191	vertexer->SetConstraintOff();
	192
	193	if (grpRecoParam) {
	194	Int_t nCutsVertexer = grpRecoParam->GetVertexerTracksNCuts();
	195	Double_t *cutsVertexer = new Double_t[nCutsVertexer];
	196	grpRecoParam->GetVertexerTracksCutsITS(cutsVertexer,nCutsVertexer);
	197	vertexer->SetCuts(cutsVertexer,nCutsVertexer);
	198	delete [] cutsVertexer; cutsVertexer = NULL;
	199	}
bbb57828	200
29be8a66	201	vertexer->SetConstraintOff();
bbb57828	202
29be8a66	203	//track vertex recomputed from the vertexer
	204	AliESDVertex *trkv = vertexer->FindPrimaryVertex(esdE);
	205	Printf ("trkv = %p", trkv);
	206	//const AliESDVertex *trkv = esdE->GetPrimaryVertexTracks();
	207
	208	//SPD vertex taken from the ESD
bbb57828	209	const AliESDVertex *spdv = esdE->GetPrimaryVertexSPD();
bbb57828	210
29be8a66	211	//ITSSA vertex computed from the vertexer
	212	vertexer->SetITSpureSA();
	213	AliESDVertex *itsSAv = vertexer->FindPrimaryVertex(esdE);
bbb57828	214
	215	if(spdv) {
	216	if(spdv->GetNContributors()>0) {
	217	TString title=spdv->GetTitle();
	218	if(title.Contains("3D")) {
	219	((TH1F*)fOutput->FindObject("hSPDVertexX"))->Fill(spdv->GetXv());
	220	((TH1F*)fOutput->FindObject("hSPDVertexY"))->Fill(spdv->GetYv());
	221	}
	222	((TH1F*)fOutput->FindObject("hSPDVertexZ"))->Fill(spdv->GetZv());
	223	}
	224	}
	225
	226
	227	if(trkv) {
	228	if(trkv->GetNContributors()>0) {
	229	((TH1F*)fOutput->FindObject("hTRKVertexX"))->Fill(trkv->GetXv());
	230	((TH1F*)fOutput->FindObject("hTRKVertexY"))->Fill(trkv->GetYv());
	231	((TH1F*)fOutput->FindObject("hTRKVertexZ"))->Fill(trkv->GetZv());
	232
	233	if (ntrklets>30 && ntrklets<45){
	234	((TH1F*)fOutput->FindObject("hTRKVertexXdefMult"))->Fill(trkv->GetXv());
	235	((TH1F*)fOutput->FindObject("hTRKVertexYdefMult"))->Fill(trkv->GetYv());
	236	}
29be8a66	237
	238	if (ntrklets>1500){
	239	((TH1F*)fOutput->FindObject("hTRKVertexXHighMult"))->Fill(trkv->GetXv());
	240	((TH1F*)fOutput->FindObject("hTRKVertexYHighMult"))->Fill(trkv->GetYv());
	241	}
bbb57828	242
	243	((TH2F*)fOutput->FindObject("hTRKVertexXZ"))->Fill(trkv->GetXv(),trkv->GetZv());
	244	((TH2F*)fOutput->FindObject("hTRKVertexYZ"))->Fill(trkv->GetYv(),trkv->GetZv());
	245
	246	}
	247	}
29be8a66	248
bbb57828	249	if (itsSAv){
bbb57828	250	if (itsSAv->GetNContributors()>0){
29be8a66	251
bbb57828	252	((TH1F*)fOutput->FindObject("hITSSAVertexX"))->Fill(itsSAv->GetXv());
	253	((TH1F*)fOutput->FindObject("hITSSAVertexY"))->Fill(itsSAv->GetYv());
	254	((TH1F*)fOutput->FindObject("hITSSAVertexZ"))->Fill(itsSAv->GetZv());
	255
29be8a66	256	if (ntrklets>30 && ntrklets<45){
	257	((TH1F*)fOutput->FindObject("hITSSAVertexXdefMult"))->Fill(itsSAv->GetXv());
	258	((TH1F*)fOutput->FindObject("hITSSAVertexYdefMult"))->Fill(itsSAv->GetYv());
	259	}
	260
	261	if (ntrklets>1500){
	262	((TH1F*)fOutput->FindObject("hITSSAVertexXHighMult"))->Fill(itsSAv->GetXv());
	263	((TH1F*)fOutput->FindObject("hITSSAVertexYHighMult"))->Fill(itsSAv->GetYv());
	264	}
bbb57828	265
	266	((TH2F*)fOutput->FindObject("hITSSAVertexXZ"))->Fill(itsSAv->GetXv(),itsSAv->GetZv());
	267	((TH2F*)fOutput->FindObject("hITSSAVertexYZ"))->Fill(itsSAv->GetYv(),itsSAv->GetZv());
29be8a66	268
bbb57828	269	}
	270	}
	271
	272	delete itsSAv;
29be8a66	273	delete vertexer;
bbb57828	274
	275	PostData(1, fOutput);
	276
	277	return;
	278
	279	}
	280
	281
	282	//________________________________________________________________________
	283	void AliMeanVertexCalibTask::Terminate(Option_t *)
	284	{
	285	// Draw result to the screen
	286	// Called once at the end of the query
	287	fOutput = dynamic_cast<TList*> (GetOutputData(1));
	288	if (!fOutput) {
	289	Printf("ERROR: fOutput not available");
	290	return;
	291	}
	292
	293
	294	return;
	295
	296	}
	297
	298
	299	//__________________________________________________________________________
29be8a66	300	// AliESDVertex* AliMeanVertexCalibTask::ReconstructPrimaryVertex(Bool_t constr,Int_t mode) const {
	301	// // On the fly reco of ITS+TPC vertex from ESD
	302	// // mode=0 use all tracks
	303	// // mode=1 use odd-number tracks
	304	// // mode=2 use even-number tracks
	305
	306	// AliESDEvent* evt = (AliESDEvent*) fInputEvent;
	307	// AliVertexerTracks vertexer(evt->GetMagneticField());
	308	// if(evt->GetNumberOfTracks()<500) {
	309	// vertexer.SetITSMode(); // defaults
	310	// vertexer.SetMinClusters(4); // default is 5
	311	// } else {
	312	// vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb
	313	// }
	314	// if (fOnlyITSSATracks) vertexer.SetITSpureSA();
bbb57828	315
29be8a66	316	// Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
	317	// Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
	318	// Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
	319	// AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
	320	// vertexer.SetVtxStart(initVertex);
	321	// delete initVertex;
	322	// if(!constr) vertexer.SetConstraintOff();
	323
	324	// if(fOnlyITSTPCTracks \|\| fOnlyITSSATracks \|\| mode>0) {
	325	// Int_t iskip=0;
	326	// Int_t ntracks = evt->GetNumberOfTracks();
	327	// Int_t *skip = new Int_t[ntracks];
	328	// for(Int_t i=0;i<ntracks;i++) skip[i]=-1;
	329	// for(Int_t itr=0;itr<ntracks; itr++) {
	330	// AliESDtrack* track = evt->GetTrack(itr);
	331	// if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
	332	// skip[iskip++]=itr;
	333	// continue;
	334	// }
	335	// if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
	336	// skip[iskip++]=itr;
	337	// continue;
	338	// }
	339	// if(mode==1 && itr%2==0) skip[iskip++]=itr;
	340	// if(mode==2 && itr%2==1) skip[iskip++]=itr;
	341	// }
	342	// vertexer.SetSkipTracks(iskip,skip);
	343	// delete [] skip; skip=NULL;
	344	// }
bbb57828	345
29be8a66	346	// return vertexer.FindPrimaryVertex(evt);
29be8a66	347	// }