[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 *hTRKVertexXvsMult = new TH2F("hTRKVertexXvsMult", "TRKVertex X vs mult", 200, -1, 1, 300, 0, 3000);
  fOutput->Add(hTRKVertexXvsMult);
  
  TH2F *hTRKVertexYvsMult = new TH2F("hTRKVertexYvsMult", "TRKVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
  fOutput->Add(hTRKVertexYvsMult);
  
  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);

  TH2F *hITSSAVertexXvsMult = new TH2F("hITSSAVertexXvsMult", "ITSSAVertex X vs mult", 200, -1, 1, 300, 0, 3000);
  fOutput->Add(hITSSAVertexXvsMult);
  
  TH2F *hITSSAVertexYvsMult = new TH2F("hITSSAVertexYvsMult", "ITSSAVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
  fOutput->Add(hITSSAVertexYvsMult);
  
  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();
  if (runNb > 0) {
    man->SetRun(runNb);
    // Printf("runNb = %d", runNb);
  }
  
  AliCDBEntry *entry = (AliCDBEntry*)man->Get("GRP/Calib/RecoParam/");
  // Printf("entry = %p", entry);
  TObjArray *arrayRecoParam=0x0;
  if (entry) {
    arrayRecoParam = (TObjArray*)entry->GetObject();
    // Printf("arrayRecoParam = %p", arrayRecoParam);
  }
  else { 
    Printf("CDBEntry not found");
    return;
  }
  AliGRPRecoParam *grpRecoParam=0x0;
  if (kLowFlux) grpRecoParam= (AliGRPRecoParam*)arrayRecoParam->At(1);
  else if (kHighFlux) grpRecoParam = (AliGRPRecoParam*)arrayRecoParam->At(2);
  
  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);
  
  //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->GetX());
	((TH1F*)fOutput->FindObject("hSPDVertexY"))->Fill(spdv->GetY());
      }
      ((TH1F*)fOutput->FindObject("hSPDVertexZ"))->Fill(spdv->GetZ());
    }
  }
  
  
  if(trkv) {
    if(trkv->GetNContributors()>0) {
      ((TH1F*)fOutput->FindObject("hTRKVertexX"))->Fill(trkv->GetX());
      ((TH1F*)fOutput->FindObject("hTRKVertexY"))->Fill(trkv->GetY());
      ((TH1F*)fOutput->FindObject("hTRKVertexZ"))->Fill(trkv->GetZ());

      ((TH2F*)fOutput->FindObject("hTRKVertexXvsMult"))->Fill(trkv->GetX(), ntrklets);
      ((TH2F*)fOutput->FindObject("hTRKVertexYvsMult"))->Fill(trkv->GetY(), ntrklets);
      
      if (ntrklets>30 && ntrklets<45){
	((TH1F*)fOutput->FindObject("hTRKVertexXdefMult"))->Fill(trkv->GetX());
	((TH1F*)fOutput->FindObject("hTRKVertexYdefMult"))->Fill(trkv->GetY());
      }
      
      if (ntrklets>1500){
	((TH1F*)fOutput->FindObject("hTRKVertexXHighMult"))->Fill(trkv->GetX());
	((TH1F*)fOutput->FindObject("hTRKVertexYHighMult"))->Fill(trkv->GetY());
      }
      
      ((TH2F*)fOutput->FindObject("hTRKVertexXZ"))->Fill(trkv->GetX(),trkv->GetZ());
      ((TH2F*)fOutput->FindObject("hTRKVertexYZ"))->Fill(trkv->GetY(),trkv->GetZ());
      
    }
  }
  
  if (itsSAv){
    if (itsSAv->GetNContributors()>0){
      
      ((TH1F*)fOutput->FindObject("hITSSAVertexX"))->Fill(itsSAv->GetX());
      ((TH1F*)fOutput->FindObject("hITSSAVertexY"))->Fill(itsSAv->GetY());
      ((TH1F*)fOutput->FindObject("hITSSAVertexZ"))->Fill(itsSAv->GetZ());

      ((TH2F*)fOutput->FindObject("hITSSAVertexXvsMult"))->Fill(itsSAv->GetX(), ntrklets);
      ((TH2F*)fOutput->FindObject("hITSSAVertexYvsMult"))->Fill(itsSAv->GetY(), ntrklets);
      
      if (ntrklets>30 && ntrklets<45){
	((TH1F*)fOutput->FindObject("hITSSAVertexXdefMult"))->Fill(itsSAv->GetX());
	((TH1F*)fOutput->FindObject("hITSSAVertexYdefMult"))->Fill(itsSAv->GetY());
      }
      
      if (ntrklets>1500){
	((TH1F*)fOutput->FindObject("hITSSAVertexXHighMult"))->Fill(itsSAv->GetX());
	((TH1F*)fOutput->FindObject("hITSSAVertexYHighMult"))->Fill(itsSAv->GetY());
      }
      
      ((TH2F*)fOutput->FindObject("hITSSAVertexXZ"))->Fill(itsSAv->GetX(),itsSAv->GetZ());
      ((TH2F*)fOutput->FindObject("hITSSAVertexYZ"))->Fill(itsSAv->GetY(),itsSAv->GetZ());
      
    }
  }

  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);
791a58dc	91
	92	TH2F *hTRKVertexXvsMult = new TH2F("hTRKVertexXvsMult", "TRKVertex X vs mult", 200, -1, 1, 300, 0, 3000);
	93	fOutput->Add(hTRKVertexXvsMult);
	94
	95	TH2F *hTRKVertexYvsMult = new TH2F("hTRKVertexYvsMult", "TRKVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
	96	fOutput->Add(hTRKVertexYvsMult);
bbb57828	97
	98	TH2F *hTRKVertexXZ = new TH2F("hTRKVertexXZ", "TRKVertex XZ corr", 200, -1, 1, 200, -20, 20);
	99	fOutput->Add(hTRKVertexXZ);
791a58dc	100
bbb57828	101	TH2F *hTRKVertexYZ = new TH2F("hTRKVertexYZ", "TRKVertex YZ corr", 200, -1, 1, 200, -20, 20);
	102	fOutput->Add(hTRKVertexYZ);
	103
	104	TH1F* hTRKVertexXdefMult = new TH1F("hTRKVertexXdefMult","TRKVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
	105	fOutput->Add(hTRKVertexXdefMult);
	106	TH1F* hTRKVertexYdefMult = new TH1F("hTRKVertexYdefMult","TRKVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
	107	fOutput->Add(hTRKVertexYdefMult);
a680415a	108
791a58dc	109	TH1F* hTRKVertexXHighMult = new TH1F("hTRKVertexXHighMult","TRKVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
	110	fOutput->Add(hTRKVertexXHighMult);
	111	TH1F* hTRKVertexYHighMult = new TH1F("hTRKVertexYHighMult","TRKVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
	112	fOutput->Add(hTRKVertexYHighMult);
	113
bbb57828	114	TH1F* hITSSAVertexX = new TH1F("hITSSAVertexX","ITSSAVertex x; x vertex [cm]; events",200,-1,1);
	115	fOutput->Add(hITSSAVertexX);
	116	TH1F* hITSSAVertexY = new TH1F("hITSSAVertexY","ITSSAVertex y; y vertex [cm]; events",200,-1,1);
	117	fOutput->Add(hITSSAVertexY);
	118	TH1F* hITSSAVertexZ = new TH1F("hITSSAVertexZ","ITSSAVertex z; z vertex [cm]; events",200,-20,20);
	119	fOutput->Add(hITSSAVertexZ);
	120
	121	TH2F *hITSSAVertexXZ = new TH2F("hITSSAVertexXZ", "ITSSAVertex XZ corr", 200, -1, 1, 200, -20, 20);
	122	fOutput->Add(hITSSAVertexXZ);
	123
	124	TH2F *hITSSAVertexYZ = new TH2F("hITSSAVertexYZ", "ITSSAVertex YZ corr", 200, -1, 1, 200, -20, 20);
	125	fOutput->Add(hITSSAVertexYZ);
791a58dc	126
	127	TH2F *hITSSAVertexXvsMult = new TH2F("hITSSAVertexXvsMult", "ITSSAVertex X vs mult", 200, -1, 1, 300, 0, 3000);
	128	fOutput->Add(hITSSAVertexXvsMult);
	129
	130	TH2F *hITSSAVertexYvsMult = new TH2F("hITSSAVertexYvsMult", "ITSSAVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
	131	fOutput->Add(hITSSAVertexYvsMult);
bbb57828	132
	133	TH1F* hITSSAVertexXdefMult = new TH1F("hITSSAVertexXdefMult","ITSSAVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
	134	fOutput->Add(hITSSAVertexXdefMult);
	135	TH1F* hITSSAVertexYdefMult = new TH1F("hITSSAVertexYdefMult","ITSSAVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
	136	fOutput->Add(hITSSAVertexYdefMult);
	137
791a58dc	138
	139	TH1F* hITSSAVertexXHighMult = new TH1F("hITSSAVertexXHighMult","ITSSAVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
	140	fOutput->Add(hITSSAVertexXHighMult);
	141	TH1F* hITSSAVertexYHighMult = new TH1F("hITSSAVertexYHighMult","ITSSAVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
	142	fOutput->Add(hITSSAVertexYHighMult);
	143
bbb57828	144	PostData(1, fOutput);
	145
	146	return;
	147	}
	148
	149
	150	//________________________________________________________________________
	151	void AliMeanVertexCalibTask::UserExec(Option_t *)
	152	{
	153	// Main loop
	154	// Called for each event
	155
	156	if (!InputEvent()) {
	157	Printf("ERROR: fESD not available");
	158	return;
	159	}
	160
	161	AliESDEvent* esdE = (AliESDEvent*) InputEvent();
	162
bbb57828	163	const AliMultiplicity *alimult = esdE->GetMultiplicity();
	164	Int_t ntrklets=0;
	165	if(alimult) ntrklets = alimult->GetNumberOfTracklets();
29be8a66	166
29be8a66	167	const char* beamType = esdE->GetBeamType();
87b471fa	168	// Printf("beam type = %s", beamType);
29be8a66	169
	170	Bool_t kLowFlux = kTRUE, kHighFlux = kFALSE;
	171	// TString pp= "p-p";
	172	//TString pA= "p-A";
	173	TString AA= "A-A";
	174
	175	if (beamType == AA){
	176	kHighFlux = kTRUE;
	177	kLowFlux = kFALSE;
87b471fa	178	// Printf ("high flux setting");
29be8a66	179	}
	180
	181	AliCDBManager* man = AliCDBManager::Instance();
9fa1a9e8	182	//man->SetDefaultStorage("raw://");
29be8a66	183	Int_t runNb = esdE->GetRunNumber();
791a58dc	184	if (runNb > 0) {
791a58dc	185	man->SetRun(runNb);
87b471fa	186	// Printf("runNb = %d", runNb);
791a58dc	187	}
791a58dc	188
29be8a66	189	AliCDBEntry entry = (AliCDBEntry)man->Get("GRP/Calib/RecoParam/");
87b471fa	190	// Printf("entry = %p", entry);
791a58dc	191	TObjArray *arrayRecoParam=0x0;
	192	if (entry) {
	193	arrayRecoParam = (TObjArray*)entry->GetObject();
87b471fa	194	// Printf("arrayRecoParam = %p", arrayRecoParam);
791a58dc	195	}
	196	else {
	197	Printf("CDBEntry not found");
	198	return;
	199	}
29be8a66	200	AliGRPRecoParam *grpRecoParam=0x0;
	201	if (kLowFlux) grpRecoParam= (AliGRPRecoParam*)arrayRecoParam->At(1);
	202	else if (kHighFlux) grpRecoParam = (AliGRPRecoParam*)arrayRecoParam->At(2);
	203
29be8a66	204	AliVertexerTracks *vertexer= new AliVertexerTracks(esdE->GetMagneticField());
	205	vertexer->SetITSMode();
	206	vertexer->SetConstraintOff();
	207
	208	if (grpRecoParam) {
	209	Int_t nCutsVertexer = grpRecoParam->GetVertexerTracksNCuts();
	210	Double_t *cutsVertexer = new Double_t[nCutsVertexer];
	211	grpRecoParam->GetVertexerTracksCutsITS(cutsVertexer,nCutsVertexer);
	212	vertexer->SetCuts(cutsVertexer,nCutsVertexer);
	213	delete [] cutsVertexer; cutsVertexer = NULL;
	214	}
bbb57828	215
29be8a66	216	vertexer->SetConstraintOff();
bbb57828	217
29be8a66	218	//track vertex recomputed from the vertexer
29be8a66	219	AliESDVertex *trkv = vertexer->FindPrimaryVertex(esdE);
87b471fa	220
29be8a66	221	//const AliESDVertex *trkv = esdE->GetPrimaryVertexTracks();
	222
	223	//SPD vertex taken from the ESD
bbb57828	224	const AliESDVertex *spdv = esdE->GetPrimaryVertexSPD();
bbb57828	225
29be8a66	226	//ITSSA vertex computed from the vertexer
	227	vertexer->SetITSpureSA();
	228	AliESDVertex *itsSAv = vertexer->FindPrimaryVertex(esdE);
bbb57828	229
	230	if(spdv) {
	231	if(spdv->GetNContributors()>0) {
	232	TString title=spdv->GetTitle();
	233	if(title.Contains("3D")) {
e690d4d0	234	((TH1F*)fOutput->FindObject("hSPDVertexX"))->Fill(spdv->GetX());
e690d4d0	235	((TH1F*)fOutput->FindObject("hSPDVertexY"))->Fill(spdv->GetY());
bbb57828	236	}
e690d4d0	237	((TH1F*)fOutput->FindObject("hSPDVertexZ"))->Fill(spdv->GetZ());
bbb57828	238	}
	239	}
	240
	241
	242	if(trkv) {
	243	if(trkv->GetNContributors()>0) {
e690d4d0	244	((TH1F*)fOutput->FindObject("hTRKVertexX"))->Fill(trkv->GetX());
	245	((TH1F*)fOutput->FindObject("hTRKVertexY"))->Fill(trkv->GetY());
	246	((TH1F*)fOutput->FindObject("hTRKVertexZ"))->Fill(trkv->GetZ());
791a58dc	247
e690d4d0	248	((TH2F*)fOutput->FindObject("hTRKVertexXvsMult"))->Fill(trkv->GetX(), ntrklets);
e690d4d0	249	((TH2F*)fOutput->FindObject("hTRKVertexYvsMult"))->Fill(trkv->GetY(), ntrklets);
bbb57828	250
bbb57828	251	if (ntrklets>30 && ntrklets<45){
e690d4d0	252	((TH1F*)fOutput->FindObject("hTRKVertexXdefMult"))->Fill(trkv->GetX());
e690d4d0	253	((TH1F*)fOutput->FindObject("hTRKVertexYdefMult"))->Fill(trkv->GetY());
bbb57828	254	}
29be8a66	255
29be8a66	256	if (ntrklets>1500){
e690d4d0	257	((TH1F*)fOutput->FindObject("hTRKVertexXHighMult"))->Fill(trkv->GetX());
e690d4d0	258	((TH1F*)fOutput->FindObject("hTRKVertexYHighMult"))->Fill(trkv->GetY());
29be8a66	259	}
bbb57828	260
e690d4d0	261	((TH2F*)fOutput->FindObject("hTRKVertexXZ"))->Fill(trkv->GetX(),trkv->GetZ());
e690d4d0	262	((TH2F*)fOutput->FindObject("hTRKVertexYZ"))->Fill(trkv->GetY(),trkv->GetZ());
bbb57828	263
	264	}
	265	}
29be8a66	266
bbb57828	267	if (itsSAv){
bbb57828	268	if (itsSAv->GetNContributors()>0){
29be8a66	269
e690d4d0	270	((TH1F*)fOutput->FindObject("hITSSAVertexX"))->Fill(itsSAv->GetX());
	271	((TH1F*)fOutput->FindObject("hITSSAVertexY"))->Fill(itsSAv->GetY());
	272	((TH1F*)fOutput->FindObject("hITSSAVertexZ"))->Fill(itsSAv->GetZ());
791a58dc	273
e690d4d0	274	((TH2F*)fOutput->FindObject("hITSSAVertexXvsMult"))->Fill(itsSAv->GetX(), ntrklets);
e690d4d0	275	((TH2F*)fOutput->FindObject("hITSSAVertexYvsMult"))->Fill(itsSAv->GetY(), ntrklets);
bbb57828	276
29be8a66	277	if (ntrklets>30 && ntrklets<45){
e690d4d0	278	((TH1F*)fOutput->FindObject("hITSSAVertexXdefMult"))->Fill(itsSAv->GetX());
e690d4d0	279	((TH1F*)fOutput->FindObject("hITSSAVertexYdefMult"))->Fill(itsSAv->GetY());
29be8a66	280	}
	281
	282	if (ntrklets>1500){
e690d4d0	283	((TH1F*)fOutput->FindObject("hITSSAVertexXHighMult"))->Fill(itsSAv->GetX());
e690d4d0	284	((TH1F*)fOutput->FindObject("hITSSAVertexYHighMult"))->Fill(itsSAv->GetY());
29be8a66	285	}
bbb57828	286
e690d4d0	287	((TH2F*)fOutput->FindObject("hITSSAVertexXZ"))->Fill(itsSAv->GetX(),itsSAv->GetZ());
e690d4d0	288	((TH2F*)fOutput->FindObject("hITSSAVertexYZ"))->Fill(itsSAv->GetY(),itsSAv->GetZ());
29be8a66	289
bbb57828	290	}
	291	}
	292
	293	delete itsSAv;
29be8a66	294	delete vertexer;
bbb57828	295
	296	PostData(1, fOutput);
	297
	298	return;
	299
	300	}
	301
	302
	303	//________________________________________________________________________
	304	void AliMeanVertexCalibTask::Terminate(Option_t *)
	305	{
	306	// Draw result to the screen
	307	// Called once at the end of the query
	308	fOutput = dynamic_cast<TList*> (GetOutputData(1));
	309	if (!fOutput) {
	310	Printf("ERROR: fOutput not available");
	311	return;
	312	}
	313
	314
	315	return;
	316
	317	}
	318
	319
	320	//__________________________________________________________________________
29be8a66	321	// AliESDVertex* AliMeanVertexCalibTask::ReconstructPrimaryVertex(Bool_t constr,Int_t mode) const {
	322	// // On the fly reco of ITS+TPC vertex from ESD
	323	// // mode=0 use all tracks
	324	// // mode=1 use odd-number tracks
	325	// // mode=2 use even-number tracks
	326
	327	// AliESDEvent* evt = (AliESDEvent*) fInputEvent;
	328	// AliVertexerTracks vertexer(evt->GetMagneticField());
	329	// if(evt->GetNumberOfTracks()<500) {
	330	// vertexer.SetITSMode(); // defaults
	331	// vertexer.SetMinClusters(4); // default is 5
	332	// } else {
	333	// vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb
	334	// }
	335	// if (fOnlyITSSATracks) vertexer.SetITSpureSA();
bbb57828	336
29be8a66	337	// Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
	338	// Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
	339	// Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
	340	// AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
	341	// vertexer.SetVtxStart(initVertex);
	342	// delete initVertex;
	343	// if(!constr) vertexer.SetConstraintOff();
	344
	345	// if(fOnlyITSTPCTracks \|\| fOnlyITSSATracks \|\| mode>0) {
	346	// Int_t iskip=0;
	347	// Int_t ntracks = evt->GetNumberOfTracks();
	348	// Int_t *skip = new Int_t[ntracks];
	349	// for(Int_t i=0;i<ntracks;i++) skip[i]=-1;
	350	// for(Int_t itr=0;itr<ntracks; itr++) {
	351	// AliESDtrack* track = evt->GetTrack(itr);
	352	// if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
	353	// skip[iskip++]=itr;
	354	// continue;
	355	// }
	356	// if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
	357	// skip[iskip++]=itr;
	358	// continue;
	359	// }
	360	// if(mode==1 && itr%2==0) skip[iskip++]=itr;
	361	// if(mode==2 && itr%2==1) skip[iskip++]=itr;
	362	// }
	363	// vertexer.SetSkipTracks(iskip,skip);
	364	// delete [] skip; skip=NULL;
	365	// }
bbb57828	366
29be8a66	367	// return vertexer.FindPrimaryVertex(evt);
29be8a66	368	// }