[u/mrichter/AliRoot.git] / PWG1 / ITS / AliAnalysisTaskITSAlignQA.cxx

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisDataContainer.h"
#include "AliITSRecPoint.h"
#include "AliESDEvent.h"
#include "AliTrackPointArray.h"
#include "AliITSgeomTGeo.h"
#include "AliITSTPArrayFit.h"
#include "AliESDfriend.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "AliITSCalibrationSDD.h"
#include "AliITSresponseSDD.h"
#include "AliGeomManager.h"
#include <TSystem.h>
#include <TTree.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TChain.h>
#include <TGeoGlobalMagField.h>
#include "AliESDInputHandlerRP.h"

/**************************************************************************
 * Copyright(c) 1998-2010, 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.                  *
 **************************************************************************/

//*************************************************************************
// Implementation of class AliAnalysiTaskITSAlignQA
// AliAnalysisTaskSE to extract from ESD + ESDfriends 
// the track-to-point residuals and dE/dx vs, time for SDD modules
//
// Author: F. Prino, prino@to.infn.it
//*************************************************************************


#include "AliAnalysisTaskITSAlignQA.h"

ClassImp(AliAnalysisTaskITSAlignQA)
//______________________________________________________________________________
AliAnalysisTaskITSAlignQA::AliAnalysisTaskITSAlignQA() : AliAnalysisTaskSE("SDD Calibration"), 
  fOutput(0),
  fHistNEvents(0),
  fHistPtAccept(0),
  fDoSPDResiduals(kTRUE),
  fDoSDDResiduals(kTRUE),
  fDoSSDResiduals(kTRUE),
  fDoSDDdEdxCalib(kTRUE),
  fUseITSsaTracks(kFALSE),
  fMinITSpts(3),
  fMinTPCpts(70),
  fMinPt(0.5),
  fNPtBins(8),
  fFitter(0),
  fRunNb(0),
  fOCDBLocation("local://$ALICE_ROOT/OCDB")
{
  //
  fFitter = new AliITSTPArrayFit(5);
  Double_t xbins[9]={0.3,0.5,0.75,1.,1.5,2.,3.,5.,10.};
  SetPtBinLimits(8,xbins);
  DefineOutput(1, TList::Class());
}


//___________________________________________________________________________
AliAnalysisTaskITSAlignQA::~AliAnalysisTaskITSAlignQA(){
  //
  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
  if(fHistNEvents){
    delete fHistNEvents;
    fHistNEvents=0;    
  }
  for(Int_t i=0;i<kNSDDmods;i++){
    if(fHistSDDResidXvsX[i]){ delete fHistSDDResidXvsX[i]; fHistSDDResidXvsX[i]=0;}
    if(fHistSDDResidXvsZ[i]){ delete fHistSDDResidXvsX[i]; fHistSDDResidXvsX[i]=0;}
    if(fHistSDDResidZvsX[i]){ delete fHistSDDResidXvsX[i]; fHistSDDResidXvsX[i]=0;}
    if(fHistSDDResidZvsZ[i]){ delete fHistSDDResidXvsX[i]; fHistSDDResidXvsX[i]=0;}
    if(fHistSDDdEdxvsDrTime[i]){ delete fHistSDDdEdxvsDrTime[i]; fHistSDDdEdxvsDrTime[i]=0;}
    if(fHistSDDDrTimeAll[i]){ delete fHistSDDDrTimeAll[i]; fHistSDDDrTimeAll[i]=0;}
    if(fHistSDDDrTimeExtra[i]){ delete fHistSDDDrTimeExtra[i]; fHistSDDDrTimeExtra[i]=0;}
    if(fHistSDDDrTimeAttac[i]){ delete fHistSDDDrTimeAttac[i]; fHistSDDDrTimeAttac[i]=0;}
  }
  if(fFitter){
    delete fFitter;
    fFitter=0;
  }
}
//___________________________________________________________________________
void AliAnalysisTaskITSAlignQA::UserCreateOutputObjects() {
  //
  LoadGeometryFromOCDB();

  fOutput = new TList();
  fOutput->SetOwner();
  fOutput->SetName("OutputHistos");

  fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
  fHistNEvents->Sumw2();
  fHistNEvents->SetMinimum(0);
  fOutput->Add(fHistNEvents);

  fHistPtAccept = new TH1F("hPtAccept","Pt distrib of accepted tracks",50,0.,5.);
  fHistPtAccept->Sumw2();
  fHistPtAccept->SetMinimum(0);
  fOutput->Add(fHistPtAccept);

  if(fDoSPDResiduals) CreateSPDHistos();
  if(fDoSDDResiduals || fDoSDDdEdxCalib) CreateSDDHistos();
  if(fDoSSDResiduals) CreateSSDHistos();

}

//___________________________________________________________________________
void AliAnalysisTaskITSAlignQA::CreateSPDHistos(){
  // Histos for SPD
  

  for(Int_t iMod=0; iMod<kNSPDmods; iMod++){
    fHistSPDResidX[iMod] = new TH2F(Form("hSPDResidX%d",iMod),
				    Form("hSPDResidX%d",iMod),
				    fNPtBins,fPtBinLimits,
				    250,-0.05,0.05);
    fHistSPDResidX[iMod]->Sumw2();
    fOutput->Add(fHistSPDResidX[iMod]);

    fHistSPDResidZ[iMod] = new TH2F(Form("hSPDResidZ%d",iMod),
				    Form("hSPDResidZ%d",iMod),
				    fNPtBins,fPtBinLimits,
				    250,-0.1,0.1);
    fHistSPDResidZ[iMod]->Sumw2();
    fOutput->Add(fHistSPDResidZ[iMod]);
  }
  return;
}

//___________________________________________________________________________
void AliAnalysisTaskITSAlignQA::CreateSDDHistos(){
  // Histos for SDD

  for(Int_t iMod=0; iMod<kNSDDmods; iMod++){
    if(fDoSDDResiduals){
      fHistSDDResidX[iMod] = new TH2F(Form("hSDDResidX%d",iMod),
				      Form("hSDDResidX%d",iMod),
				      fNPtBins,fPtBinLimits,
				      300,-0.15,0.15);
      fHistSDDResidX[iMod]->Sumw2();
      fOutput->Add(fHistSDDResidX[iMod]);
      
      fHistSDDResidZ[iMod] = new TH2F(Form("hSDDResidZ%d",iMod),
				      Form("hSDDResidZ%d",iMod),
				      fNPtBins,fPtBinLimits,
				      200,-0.1,0.1);
      fHistSDDResidZ[iMod]->Sumw2();
      fOutput->Add(fHistSDDResidZ[iMod]);
      
      fHistSDDResidXvsX[iMod] = new TH2F(Form("hSDDResidXvsX%d",iMod+kNSPDmods),
					 Form("hSDDResidXvsX%d",iMod+kNSPDmods),
					 40,-3.5,3.5,300,-0.15,0.15);   
      fHistSDDResidXvsX[iMod]->Sumw2();
      fOutput->Add(fHistSDDResidXvsX[iMod]);
      
      fHistSDDResidXvsZ[iMod] = new TH2F(Form("hSDDResidXvsZ%d",iMod+kNSPDmods),
					 Form("hSDDResidXvsZ%d",iMod+kNSPDmods),
					 10,-3.8,3.8,300,-0.15,0.15);   
      fHistSDDResidXvsZ[iMod]->Sumw2();
      fOutput->Add(fHistSDDResidXvsZ[iMod]);
      
      fHistSDDResidZvsX[iMod] = new TH2F(Form("hSDDResidZvsX%d",iMod+kNSPDmods),
				      Form("hSDDResidZvsX%d",iMod+kNSPDmods),
				      40,-3.5,3.5,200,-0.1,0.1);   
      fHistSDDResidZvsX[iMod]->Sumw2();
      fOutput->Add(fHistSDDResidZvsX[iMod]);
      
      fHistSDDResidZvsZ[iMod] = new TH2F(Form("hSDDResidZvsZ%d",iMod+kNSPDmods),
					 Form("hSDDResidZvsZ%d",iMod+kNSPDmods),
					 10,-3.8,3.8,200,-0.1,0.1);   
      fHistSDDResidZvsZ[iMod]->Sumw2();
      fOutput->Add(fHistSDDResidZvsZ[iMod]);
    }

    if(fDoSDDdEdxCalib){
      fHistSDDdEdxvsDrTime[iMod] = new TH2F(Form("hSDDdEdxvsDrTime%d",iMod+kNSPDmods),
					    Form("hSDDdEdxvsDrTime%d",iMod+kNSPDmods),
					    16,0.,6400.,100,0.,300.);
      fHistSDDdEdxvsDrTime[iMod]->Sumw2();
      fOutput->Add(fHistSDDdEdxvsDrTime[iMod]);
    }

    fHistSDDDrTimeAll[iMod]=new TH1F(Form("hSDDDrTimeAll%d",iMod+kNSPDmods),
				     Form("hSDDDrTimeAll%d",iMod+kNSPDmods),
				     3200,0.,6400.);
    fHistSDDDrTimeAll[iMod]->Sumw2();
    fHistSDDDrTimeAll[iMod]->SetMinimum(0.);
    fOutput->Add(fHistSDDDrTimeAll[iMod]);

    fHistSDDDrTimeExtra[iMod]=new TH1F(Form("hSDDDrTimeExtra%d",iMod+kNSPDmods),
				       Form("hSDDDrTimeExtra%d",iMod+kNSPDmods),
				       3200,0.,6400.);
    fHistSDDDrTimeExtra[iMod]->Sumw2();
    fHistSDDDrTimeExtra[iMod]->SetMinimum(0.);
    fOutput->Add(fHistSDDDrTimeExtra[iMod]);

    fHistSDDDrTimeAttac[iMod]=new TH1F(Form("hSDDDrTimeAttac%d",iMod+kNSPDmods),
				       Form("hSDDDrTimeAttac%d",iMod+kNSPDmods),
				       3200,0.,6400.);
    fHistSDDDrTimeAttac[iMod]->Sumw2();
    fHistSDDDrTimeAttac[iMod]->SetMinimum(0.);
    fOutput->Add(fHistSDDDrTimeAttac[iMod]);

  }
  return;

}
//___________________________________________________________________________
void AliAnalysisTaskITSAlignQA::CreateSSDHistos(){
  // Histos for SSD
  for(Int_t iMod=0; iMod<kNSSDmods; iMod++){
    fHistSSDResidX[iMod] = new TH2F(Form("hSSDResidX%d",iMod+kNSPDmods+kNSDDmods),
				    Form("hSSDResidX%d",iMod+kNSPDmods+kNSDDmods),
				    fNPtBins,fPtBinLimits,
				    250,-0.1,0.1);
    fHistSSDResidX[iMod]->Sumw2();
    fOutput->Add(fHistSSDResidX[iMod]);

    fHistSSDResidZ[iMod] = new TH2F(Form("hSSDResidZ%d",iMod+kNSPDmods+kNSDDmods),
				    Form("hSSDResidZ%d",iMod+kNSPDmods+kNSDDmods),
				    fNPtBins,fPtBinLimits,
				    250,-1.,1.);
    fHistSSDResidZ[iMod]->Sumw2();
    fOutput->Add(fHistSSDResidZ[iMod]);
  }
  return;
}
//______________________________________________________________________________
void AliAnalysisTaskITSAlignQA::UserExec(Option_t *)
{
  //
  AliESDEvent *esd = (AliESDEvent*) (InputEvent());

  if(!esd) {
    printf("AliAnalysisTaskITSAlignQA::Exec(): bad ESD\n");
    return;
  } 


  if(!ESDfriend()) {
    printf("AliAnalysisTaskITSAlignQA::Exec(): bad ESDfriend\n");
    return;
  }

  fHistNEvents->Fill(0);
  fFitter->SetBz(esd->GetMagneticField());

  const AliTrackPointArray *array = 0;
  Int_t ntracks = esd->GetNumberOfTracks();

  for (Int_t itrack=0; itrack < ntracks; itrack++) {
    AliESDtrack * track = esd->GetTrack(itrack);
    if(!track) continue;
    if(!AcceptTrack(track)) continue;
    array = track->GetTrackPointArray();
    if(!array) continue;

    Int_t npts=array->GetNPoints();
    if(fDoSPDResiduals){ 
      FitAndFillSPD(1,array,npts,track);
      FitAndFillSPD(2,array,npts,track);
    }
    if(fDoSDDResiduals || fDoSDDdEdxCalib) FitAndFillSDD(array,npts,track);
    if(fDoSSDResiduals){ 
      FitAndFillSSD(5,array,npts,track);
      FitAndFillSSD(6,array,npts,track);
    }
  }

  PostData(1,fOutput);
  
}
//___________________________________________________________________________
Bool_t AliAnalysisTaskITSAlignQA::AcceptTrack(AliESDtrack * track){
  // track selection cuts
  Bool_t accept=kTRUE;
  if(fUseITSsaTracks){ 
      if(track->GetNcls(1)>0) accept=kFALSE;
  }else{
    if(track->GetNcls(1)<fMinTPCpts) accept=kFALSE;
  }
  if(track->GetNcls(0) < fMinITSpts) accept=kFALSE;
  Int_t trstatus=track->GetStatus();
  if(!(trstatus&AliESDtrack::kITSrefit)) accept=kFALSE;
  Float_t pt=track->Pt();
  if(pt<fMinPt) accept=kFALSE;
  if(accept) fHistPtAccept->Fill(pt);
  return accept;
}
//___________________________________________________________________________
void AliAnalysisTaskITSAlignQA::FitAndFillSPD(Int_t iLayer, const AliTrackPointArray *array, Int_t npts,AliESDtrack * track){
  // fit track and fills histos for SPD
  fFitter->AttachPoints(array,0, npts-1); 
  Int_t iPtSPD[4],modIdSPD[4];
  Int_t nPtSPD=0;
  Double_t resGlo[3],resLoc[3];
  for(Int_t ipt=0; ipt<npts; ipt++) {
    AliTrackPoint point;
    Int_t modId;
    array->GetPoint(point,ipt);
    Int_t volId = point.GetVolumeID();
    Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
    if(layerId==iLayer){
      modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
      iPtSPD[nPtSPD] = ipt;
      modIdSPD[nPtSPD] = modId;
      ++nPtSPD;
      fFitter->SetCovIScale(ipt,1e-4); 
    }
  }
  if(nPtSPD>0){
    fFitter->Fit(track->Charge(),track->Pt(),0.);
    Double_t chi2=fFitter->GetChi2NDF();
    if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
    for (Int_t ip=0; ip<nPtSPD;ip++) {
      fFitter->GetResiduals(resGlo,iPtSPD[ip]);
      TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSPD[ip]);
      mcurr->MasterToLocalVect(resGlo,resLoc);
      Int_t index=modIdSPD[ip];
      fHistSPDResidX[index]->Fill(track->Pt(),resLoc[0]);
      fHistSPDResidZ[index]->Fill(track->Pt(),resLoc[2]);
    }
  }    
}
//___________________________________________________________________________
void AliAnalysisTaskITSAlignQA::FitAndFillSDD(const AliTrackPointArray *array, Int_t npts, AliESDtrack * track){
  // fit track and fills histos for SDD
  Double_t dedx[4];
  track->GetITSdEdxSamples(dedx);

  fFitter->AttachPoints(array,0, npts-1); 
  Int_t iPtSDD[4],modIdSDD[4];
  Double_t xLocSDD[4],zLocSDD[4];
  Int_t nPtSDD=0;
  Int_t nPtSSDSPD=0;
  Double_t resGlo[3],resLoc[3];
  Float_t  posGloF[3];
  Double_t posGlo[3],posLoc[3];
  for(Int_t ipt=0; ipt<npts; ipt++) {
    AliTrackPoint point;
    Int_t modId;
    array->GetPoint(point,ipt);
    Int_t volId = point.GetVolumeID();
    Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
    if(layerId==3 || layerId==4){
      Float_t drTime=point.GetDriftTime();
      modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
      Int_t index=modId-kNSPDmods;
      fHistSDDDrTimeAll[index]->Fill(drTime);
      if(point.IsExtra()) fHistSDDDrTimeExtra[index]->Fill(drTime);
      else fHistSDDDrTimeAttac[index]->Fill(drTime);
      Float_t dedxLay=dedx[layerId-3];
      if(dedxLay>1.) fHistSDDdEdxvsDrTime[index]->Fill(drTime,dedxLay);
      iPtSDD[nPtSDD] = ipt;
      modIdSDD[nPtSDD] = modId;
      point.GetXYZ(posGloF);
      for(Int_t icoor=0;icoor<3;icoor++) posGlo[icoor]=posGloF[icoor];
      AliITSgeomTGeo::GlobalToLocal(modId,posGlo,posLoc);
      xLocSDD[nPtSDD]=posLoc[0];
      zLocSDD[nPtSDD]=posLoc[2];
      ++nPtSDD;
      fFitter->SetCovIScale(ipt,1e-4); // scaling for inverted errors of SDD
    }else{
      ++nPtSSDSPD;
    }
  }
  if(nPtSDD>0 && nPtSSDSPD>=2){
    fFitter->Fit(track->Charge(),track->Pt(),0.);
    Double_t chi2=fFitter->GetChi2NDF();
    if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
    for (Int_t ip=0; ip<nPtSDD;ip++) {
      fFitter->GetResiduals(resGlo,iPtSDD[ip]);
      TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSDD[ip]);
      mcurr->MasterToLocalVect(resGlo,resLoc);
      Int_t index=modIdSDD[ip]-kNSPDmods;
      fHistSDDResidX[index]->Fill(track->Pt(),resLoc[0]);
      fHistSDDResidZ[index]->Fill(track->Pt(),resLoc[2]);
      fHistSDDResidXvsX[index]->Fill(xLocSDD[ip],resLoc[0]);
      fHistSDDResidXvsZ[index]->Fill(zLocSDD[ip],resLoc[0]);
      fHistSDDResidZvsX[index]->Fill(xLocSDD[ip],resLoc[2]);
      fHistSDDResidZvsZ[index]->Fill(zLocSDD[ip],resLoc[2]);
    }
  }
}
//___________________________________________________________________________
void AliAnalysisTaskITSAlignQA::FitAndFillSSD(Int_t iLayer, const AliTrackPointArray *array, Int_t npts,AliESDtrack * track){
  // fit track and fills histos for SSD
  fFitter->AttachPoints(array,0, npts-1); 
  Int_t iPtSSD[4],modIdSSD[4];
  Int_t nPtSSD=0;
  Double_t resGlo[3],resLoc[3];
  for(Int_t ipt=0; ipt<npts; ipt++) {
    AliTrackPoint point;
    Int_t modId;
    array->GetPoint(point,ipt);
    Int_t volId = point.GetVolumeID();
    Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
    if(layerId==iLayer){
      modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
      iPtSSD[nPtSSD] = ipt;
      modIdSSD[nPtSSD] = modId;
      ++nPtSSD;
      fFitter->SetCovIScale(ipt,1e-4); 
    }  
  }
  if(nPtSSD>0){
    fFitter->Fit(track->Charge(),track->Pt(),0.);
    Double_t chi2=fFitter->GetChi2NDF();
    if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
    for (Int_t ip=0; ip<nPtSSD;ip++) {
      fFitter->GetResiduals(resGlo,iPtSSD[ip]);
      TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSSD[ip]);
      mcurr->MasterToLocalVect(resGlo,resLoc);
      Int_t index=modIdSSD[ip]-kNSPDmods-kNSDDmods;
      fHistSSDResidX[index]->Fill(track->Pt(),resLoc[0]);
      fHistSSDResidZ[index]->Fill(track->Pt(),resLoc[2]);
    }
  }
}
//______________________________________________________________________________
void AliAnalysisTaskITSAlignQA::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  fOutput = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutput) {     
    printf("ERROR: fOutput not available\n");
    return;
  }

  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
  return;
}


//___________________________________________________________________________
void AliAnalysisTaskITSAlignQA::LoadGeometryFromOCDB(){
  //method to get the gGeomanager
  // it is called at the CreatedOutputObject stage
  // to comply with the CAF environment
  AliInfo("Loading geometry");

  AliCDBManager *man = AliCDBManager::Instance();
  man->SetDefaultStorage(fOCDBLocation.Data());
  man->SetRun(fRunNb);
  AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
  AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
  AliGeomManager::GetNalignable("ITS");
  AliGeomManager::ApplyAlignObjsFromCDB("ITS");
}
Commit	Line	Data
89b48e9e	1	#include "AliAnalysisTask.h"
	2	#include "AliAnalysisManager.h"
	3	#include "AliAnalysisDataContainer.h"
	4	#include "AliITSRecPoint.h"
	5	#include "AliESDEvent.h"
	6	#include "AliTrackPointArray.h"
	7	#include "AliITSgeomTGeo.h"
	8	#include "AliITSTPArrayFit.h"
	9	#include "AliESDfriend.h"
	10	#include "AliCDBManager.h"
	11	#include "AliCDBEntry.h"
	12	#include "AliITSCalibrationSDD.h"
	13	#include "AliITSresponseSDD.h"
	14	#include "AliGeomManager.h"
	15	#include <TSystem.h>
	16	#include <TTree.h>
	17	#include <TH1F.h>
	18	#include <TH2F.h>
	19	#include <TChain.h>
	20	#include <TGeoGlobalMagField.h>
	21	#include "AliESDInputHandlerRP.h"
	22
	23	/**************************************************************************
	24	* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
	25	* *
	26	* Author: The ALICE Off-line Project. *
	27	* Contributors are mentioned in the code where appropriate. *
	28	* *
	29	* Permission to use, copy, modify and distribute this software and its *
	30	* documentation strictly for non-commercial purposes is hereby granted *
	31	* without fee, provided that the above copyright notice appears in all *
	32	* copies and that both the copyright notice and this permission notice *
	33	* appear in the supporting documentation. The authors make no claims *
	34	* about the suitability of this software for any purpose. It is *
	35	* provided "as is" without express or implied warranty. *
	36	**************************************************************************/
	37
	38	//*************************************************************************
	39	// Implementation of class AliAnalysiTaskITSAlignQA
	40	// AliAnalysisTaskSE to extract from ESD + ESDfriends
	41	// the track-to-point residuals and dE/dx vs, time for SDD modules
	42	//
	43	// Author: F. Prino, prino@to.infn.it
	44	//*************************************************************************
	45
	46
	47	#include "AliAnalysisTaskITSAlignQA.h"
	48
	49	ClassImp(AliAnalysisTaskITSAlignQA)
	50	//______________________________________________________________________________
	51	AliAnalysisTaskITSAlignQA::AliAnalysisTaskITSAlignQA() : AliAnalysisTaskSE("SDD Calibration"),
	52	fOutput(0),
	53	fHistNEvents(0),
	54	fHistPtAccept(0),
	55	fDoSPDResiduals(kTRUE),
	56	fDoSDDResiduals(kTRUE),
	57	fDoSSDResiduals(kTRUE),
	58	fDoSDDdEdxCalib(kTRUE),
	59	fUseITSsaTracks(kFALSE),
	60	fMinITSpts(3),
	61	fMinTPCpts(70),
	62	fMinPt(0.5),
	63	fNPtBins(8),
	64	fFitter(0),
65	fRunNb(0),
66	fOCDBLocation("local://$ALICE_ROOT/OCDB")
67	{
68	//
69	fFitter = new AliITSTPArrayFit(5);
70	Double_t xbins[9]={0.3,0.5,0.75,1.,1.5,2.,3.,5.,10.};
71	SetPtBinLimits(8,xbins);
72	DefineOutput(1, TList::Class());
73	}
74
75
76	//___________________________________________________________________________
77	AliAnalysisTaskITSAlignQA::~AliAnalysisTaskITSAlignQA(){
78	//
79	if (fOutput) {
80	delete fOutput;
81	fOutput = 0;
82	}
83	if(fHistNEvents){
84	delete fHistNEvents;
85	fHistNEvents=0;
86	}
87	for(Int_t i=0;i<kNSDDmods;i++){
88	if(fHistSDDResidXvsX[i]){ delete fHistSDDResidXvsX[i]; fHistSDDResidXvsX[i]=0;}
89	if(fHistSDDResidXvsZ[i]){ delete fHistSDDResidXvsX[i]; fHistSDDResidXvsX[i]=0;}
90	if(fHistSDDResidZvsX[i]){ delete fHistSDDResidXvsX[i]; fHistSDDResidXvsX[i]=0;}
91	if(fHistSDDResidZvsZ[i]){ delete fHistSDDResidXvsX[i]; fHistSDDResidXvsX[i]=0;}
92	if(fHistSDDdEdxvsDrTime[i]){ delete fHistSDDdEdxvsDrTime[i]; fHistSDDdEdxvsDrTime[i]=0;}
93	if(fHistSDDDrTimeAll[i]){ delete fHistSDDDrTimeAll[i]; fHistSDDDrTimeAll[i]=0;}
94	if(fHistSDDDrTimeExtra[i]){ delete fHistSDDDrTimeExtra[i]; fHistSDDDrTimeExtra[i]=0;}
95	if(fHistSDDDrTimeAttac[i]){ delete fHistSDDDrTimeAttac[i]; fHistSDDDrTimeAttac[i]=0;}
96	}
97	if(fFitter){
98	delete fFitter;
99	fFitter=0;
100	}
101	}
102	//___________________________________________________________________________
103	void AliAnalysisTaskITSAlignQA::UserCreateOutputObjects() {
104	//
105	LoadGeometryFromOCDB();
106
107	fOutput = new TList();
108	fOutput->SetOwner();
109	fOutput->SetName("OutputHistos");
110
111	fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
112	fHistNEvents->Sumw2();
113	fHistNEvents->SetMinimum(0);
114	fOutput->Add(fHistNEvents);
115
116	fHistPtAccept = new TH1F("hPtAccept","Pt distrib of accepted tracks",50,0.,5.);
117	fHistPtAccept->Sumw2();
118	fHistPtAccept->SetMinimum(0);
119	fOutput->Add(fHistPtAccept);
120
121	if(fDoSPDResiduals) CreateSPDHistos();
122	if(fDoSDDResiduals \|\| fDoSDDdEdxCalib) CreateSDDHistos();
123	if(fDoSSDResiduals) CreateSSDHistos();
124
125	}
126
127	//___________________________________________________________________________
128	void AliAnalysisTaskITSAlignQA::CreateSPDHistos(){
129	// Histos for SPD
130
131
132	for(Int_t iMod=0; iMod<kNSPDmods; iMod++){
133	fHistSPDResidX[iMod] = new TH2F(Form("hSPDResidX%d",iMod),
134	Form("hSPDResidX%d",iMod),
135	fNPtBins,fPtBinLimits,
136	250,-0.05,0.05);
137	fHistSPDResidX[iMod]->Sumw2();
138	fOutput->Add(fHistSPDResidX[iMod]);
139
140	fHistSPDResidZ[iMod] = new TH2F(Form("hSPDResidZ%d",iMod),
141	Form("hSPDResidZ%d",iMod),
142	fNPtBins,fPtBinLimits,
143	250,-0.1,0.1);
144	fHistSPDResidZ[iMod]->Sumw2();
145	fOutput->Add(fHistSPDResidZ[iMod]);
146	}
147	return;
148	}
149
150	//___________________________________________________________________________
151	void AliAnalysisTaskITSAlignQA::CreateSDDHistos(){
152	// Histos for SDD
153
154	for(Int_t iMod=0; iMod<kNSDDmods; iMod++){
155	if(fDoSDDResiduals){
156	fHistSDDResidX[iMod] = new TH2F(Form("hSDDResidX%d",iMod),
157	Form("hSDDResidX%d",iMod),
158	fNPtBins,fPtBinLimits,
159	300,-0.15,0.15);
160	fHistSDDResidX[iMod]->Sumw2();
161	fOutput->Add(fHistSDDResidX[iMod]);
162
163	fHistSDDResidZ[iMod] = new TH2F(Form("hSDDResidZ%d",iMod),
164	Form("hSDDResidZ%d",iMod),
165	fNPtBins,fPtBinLimits,
166	200,-0.1,0.1);
167	fHistSDDResidZ[iMod]->Sumw2();
168	fOutput->Add(fHistSDDResidZ[iMod]);
169
170	fHistSDDResidXvsX[iMod] = new TH2F(Form("hSDDResidXvsX%d",iMod+kNSPDmods),
171	Form("hSDDResidXvsX%d",iMod+kNSPDmods),
172	40,-3.5,3.5,300,-0.15,0.15);
173	fHistSDDResidXvsX[iMod]->Sumw2();
174	fOutput->Add(fHistSDDResidXvsX[iMod]);
175
176	fHistSDDResidXvsZ[iMod] = new TH2F(Form("hSDDResidXvsZ%d",iMod+kNSPDmods),
177	Form("hSDDResidXvsZ%d",iMod+kNSPDmods),
178	10,-3.8,3.8,300,-0.15,0.15);
179	fHistSDDResidXvsZ[iMod]->Sumw2();
180	fOutput->Add(fHistSDDResidXvsZ[iMod]);
181
182	fHistSDDResidZvsX[iMod] = new TH2F(Form("hSDDResidZvsX%d",iMod+kNSPDmods),
183	Form("hSDDResidZvsX%d",iMod+kNSPDmods),
184	40,-3.5,3.5,200,-0.1,0.1);
185	fHistSDDResidZvsX[iMod]->Sumw2();
186	fOutput->Add(fHistSDDResidZvsX[iMod]);
187
188	fHistSDDResidZvsZ[iMod] = new TH2F(Form("hSDDResidZvsZ%d",iMod+kNSPDmods),
189	Form("hSDDResidZvsZ%d",iMod+kNSPDmods),
190	10,-3.8,3.8,200,-0.1,0.1);
191	fHistSDDResidZvsZ[iMod]->Sumw2();
192	fOutput->Add(fHistSDDResidZvsZ[iMod]);
193	}
194
195	if(fDoSDDdEdxCalib){
196	fHistSDDdEdxvsDrTime[iMod] = new TH2F(Form("hSDDdEdxvsDrTime%d",iMod+kNSPDmods),
197	Form("hSDDdEdxvsDrTime%d",iMod+kNSPDmods),
198	16,0.,6400.,100,0.,300.);
199	fHistSDDdEdxvsDrTime[iMod]->Sumw2();
200	fOutput->Add(fHistSDDdEdxvsDrTime[iMod]);
201	}
202
203	fHistSDDDrTimeAll[iMod]=new TH1F(Form("hSDDDrTimeAll%d",iMod+kNSPDmods),
204	Form("hSDDDrTimeAll%d",iMod+kNSPDmods),
205	3200,0.,6400.);
206	fHistSDDDrTimeAll[iMod]->Sumw2();
207	fHistSDDDrTimeAll[iMod]->SetMinimum(0.);
208	fOutput->Add(fHistSDDDrTimeAll[iMod]);
209
210	fHistSDDDrTimeExtra[iMod]=new TH1F(Form("hSDDDrTimeExtra%d",iMod+kNSPDmods),
211	Form("hSDDDrTimeExtra%d",iMod+kNSPDmods),
212	3200,0.,6400.);
213	fHistSDDDrTimeExtra[iMod]->Sumw2();
214	fHistSDDDrTimeExtra[iMod]->SetMinimum(0.);
215	fOutput->Add(fHistSDDDrTimeExtra[iMod]);
216
217	fHistSDDDrTimeAttac[iMod]=new TH1F(Form("hSDDDrTimeAttac%d",iMod+kNSPDmods),
218	Form("hSDDDrTimeAttac%d",iMod+kNSPDmods),
219	3200,0.,6400.);
220	fHistSDDDrTimeAttac[iMod]->Sumw2();
221	fHistSDDDrTimeAttac[iMod]->SetMinimum(0.);
222	fOutput->Add(fHistSDDDrTimeAttac[iMod]);
223
224	}
225	return;
226
227	}
228	//___________________________________________________________________________
229	void AliAnalysisTaskITSAlignQA::CreateSSDHistos(){
230	// Histos for SSD
231	for(Int_t iMod=0; iMod<kNSSDmods; iMod++){
232	fHistSSDResidX[iMod] = new TH2F(Form("hSSDResidX%d",iMod+kNSPDmods+kNSDDmods),
233	Form("hSSDResidX%d",iMod+kNSPDmods+kNSDDmods),
234	fNPtBins,fPtBinLimits,
235	250,-0.1,0.1);
236	fHistSSDResidX[iMod]->Sumw2();
237	fOutput->Add(fHistSSDResidX[iMod]);
238
239	fHistSSDResidZ[iMod] = new TH2F(Form("hSSDResidZ%d",iMod+kNSPDmods+kNSDDmods),
240	Form("hSSDResidZ%d",iMod+kNSPDmods+kNSDDmods),
241	fNPtBins,fPtBinLimits,
242	250,-1.,1.);
243	fHistSSDResidZ[iMod]->Sumw2();
244	fOutput->Add(fHistSSDResidZ[iMod]);
245	}
246	return;
247	}
248	//______________________________________________________________________________
249	void AliAnalysisTaskITSAlignQA::UserExec(Option_t *)
250	{
251	//
252	AliESDEvent esd = (AliESDEvent) (InputEvent());
253
254	if(!esd) {
255	printf("AliAnalysisTaskITSAlignQA::Exec(): bad ESD\n");
256	return;
257	}
258
259
260	if(!ESDfriend()) {
261	printf("AliAnalysisTaskITSAlignQA::Exec(): bad ESDfriend\n");
262	return;
263	}
264
265	fHistNEvents->Fill(0);
266	fFitter->SetBz(esd->GetMagneticField());
267
268	const AliTrackPointArray *array = 0;
269	Int_t ntracks = esd->GetNumberOfTracks();
270
271	for (Int_t itrack=0; itrack < ntracks; itrack++) {
272	AliESDtrack * track = esd->GetTrack(itrack);
273	if(!track) continue;
274	if(!AcceptTrack(track)) continue;
275	array = track->GetTrackPointArray();
276	if(!array) continue;
277
278	Int_t npts=array->GetNPoints();
279	if(fDoSPDResiduals){
280	FitAndFillSPD(1,array,npts,track);
281	FitAndFillSPD(2,array,npts,track);
282	}
283	if(fDoSDDResiduals \|\| fDoSDDdEdxCalib) FitAndFillSDD(array,npts,track);
284	if(fDoSSDResiduals){
285	FitAndFillSSD(5,array,npts,track);
286	FitAndFillSSD(6,array,npts,track);
287	}
288	}
289
290	PostData(1,fOutput);
291
292	}
293	//___________________________________________________________________________
294	Bool_t AliAnalysisTaskITSAlignQA::AcceptTrack(AliESDtrack * track){
295	// track selection cuts
296	Bool_t accept=kTRUE;
297	if(fUseITSsaTracks){
298	if(track->GetNcls(1)>0) accept=kFALSE;
299	}else{
300	if(track->GetNcls(1)<fMinTPCpts) accept=kFALSE;
301	}
302	if(track->GetNcls(0) < fMinITSpts) accept=kFALSE;
303	Int_t trstatus=track->GetStatus();
304	if(!(trstatus&AliESDtrack::kITSrefit)) accept=kFALSE;
305	Float_t pt=track->Pt();
306	if(pt<fMinPt) accept=kFALSE;
307	if(accept) fHistPtAccept->Fill(pt);
308	return accept;
309	}
310	//___________________________________________________________________________
311	void AliAnalysisTaskITSAlignQA::FitAndFillSPD(Int_t iLayer, const AliTrackPointArray array, Int_t npts,AliESDtrack track){
312	// fit track and fills histos for SPD
313	fFitter->AttachPoints(array,0, npts-1);
314	Int_t iPtSPD[4],modIdSPD[4];
315	Int_t nPtSPD=0;
316	Double_t resGlo[3],resLoc[3];
317	for(Int_t ipt=0; ipt<npts; ipt++) {
318	AliTrackPoint point;
319	Int_t modId;
320	array->GetPoint(point,ipt);
321	Int_t volId = point.GetVolumeID();
322	Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
323	if(layerId==iLayer){
324	modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
325	iPtSPD[nPtSPD] = ipt;
326	modIdSPD[nPtSPD] = modId;
327	++nPtSPD;
328	fFitter->SetCovIScale(ipt,1e-4);
329	}
330	}
331	if(nPtSPD>0){
332	fFitter->Fit(track->Charge(),track->Pt(),0.);
333	Double_t chi2=fFitter->GetChi2NDF();
334	if ( chi2<0 \|\| chi2>1e4 ) return; // fit failed, abandon this track
335	for (Int_t ip=0; ip<nPtSPD;ip++) {
336	fFitter->GetResiduals(resGlo,iPtSPD[ip]);
337	TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSPD[ip]);
338	mcurr->MasterToLocalVect(resGlo,resLoc);
339	Int_t index=modIdSPD[ip];
340	fHistSPDResidX[index]->Fill(track->Pt(),resLoc[0]);
341	fHistSPDResidZ[index]->Fill(track->Pt(),resLoc[2]);
342	}
343	}
344	}
345	//___________________________________________________________________________
346	void AliAnalysisTaskITSAlignQA::FitAndFillSDD(const AliTrackPointArray array, Int_t npts, AliESDtrack track){
347	// fit track and fills histos for SDD
348	Double_t dedx[4];
349	track->GetITSdEdxSamples(dedx);
350
351	fFitter->AttachPoints(array,0, npts-1);
352	Int_t iPtSDD[4],modIdSDD[4];
353	Double_t xLocSDD[4],zLocSDD[4];
354	Int_t nPtSDD=0;
355	Int_t nPtSSDSPD=0;
356	Double_t resGlo[3],resLoc[3];
357	Float_t posGloF[3];
358	Double_t posGlo[3],posLoc[3];
359	for(Int_t ipt=0; ipt<npts; ipt++) {
360	AliTrackPoint point;
361	Int_t modId;
362	array->GetPoint(point,ipt);
363	Int_t volId = point.GetVolumeID();
364	Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
365	if(layerId==3 \|\| layerId==4){
366	Float_t drTime=point.GetDriftTime();
367	modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
368	Int_t index=modId-kNSPDmods;
369	fHistSDDDrTimeAll[index]->Fill(drTime);
370	if(point.IsExtra()) fHistSDDDrTimeExtra[index]->Fill(drTime);
371	else fHistSDDDrTimeAttac[index]->Fill(drTime);
372	Float_t dedxLay=dedx[layerId-3];
373	if(dedxLay>1.) fHistSDDdEdxvsDrTime[index]->Fill(drTime,dedxLay);
374	iPtSDD[nPtSDD] = ipt;
375	modIdSDD[nPtSDD] = modId;
376	point.GetXYZ(posGloF);
377	for(Int_t icoor=0;icoor<3;icoor++) posGlo[icoor]=posGloF[icoor];
378	AliITSgeomTGeo::GlobalToLocal(modId,posGlo,posLoc);
379	xLocSDD[nPtSDD]=posLoc[0];
380	zLocSDD[nPtSDD]=posLoc[2];
381	++nPtSDD;
382	fFitter->SetCovIScale(ipt,1e-4); // scaling for inverted errors of SDD
383	}else{
384	++nPtSSDSPD;
385	}
386	}
387	if(nPtSDD>0 && nPtSSDSPD>=2){
388	fFitter->Fit(track->Charge(),track->Pt(),0.);
389	Double_t chi2=fFitter->GetChi2NDF();
390	if ( chi2<0 \|\| chi2>1e4 ) return; // fit failed, abandon this track
391	for (Int_t ip=0; ip<nPtSDD;ip++) {
392	fFitter->GetResiduals(resGlo,iPtSDD[ip]);
393	TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSDD[ip]);
394	mcurr->MasterToLocalVect(resGlo,resLoc);
395	Int_t index=modIdSDD[ip]-kNSPDmods;
396	fHistSDDResidX[index]->Fill(track->Pt(),resLoc[0]);
397	fHistSDDResidZ[index]->Fill(track->Pt(),resLoc[2]);
398	fHistSDDResidXvsX[index]->Fill(xLocSDD[ip],resLoc[0]);
399	fHistSDDResidXvsZ[index]->Fill(zLocSDD[ip],resLoc[0]);
400	fHistSDDResidZvsX[index]->Fill(xLocSDD[ip],resLoc[2]);
401	fHistSDDResidZvsZ[index]->Fill(zLocSDD[ip],resLoc[2]);
402	}
403	}
404	}
405	//___________________________________________________________________________
406	void AliAnalysisTaskITSAlignQA::FitAndFillSSD(Int_t iLayer, const AliTrackPointArray array, Int_t npts,AliESDtrack track){
407	// fit track and fills histos for SSD
408	fFitter->AttachPoints(array,0, npts-1);
409	Int_t iPtSSD[4],modIdSSD[4];
410	Int_t nPtSSD=0;
411	Double_t resGlo[3],resLoc[3];
412	for(Int_t ipt=0; ipt<npts; ipt++) {
413	AliTrackPoint point;
414	Int_t modId;
415	array->GetPoint(point,ipt);
416	Int_t volId = point.GetVolumeID();
417	Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
418	if(layerId==iLayer){
419	modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
420	iPtSSD[nPtSSD] = ipt;
421	modIdSSD[nPtSSD] = modId;
422	++nPtSSD;
423	fFitter->SetCovIScale(ipt,1e-4);
424	}
425	}
426	if(nPtSSD>0){
427	fFitter->Fit(track->Charge(),track->Pt(),0.);
428	Double_t chi2=fFitter->GetChi2NDF();
429	if ( chi2<0 \|\| chi2>1e4 ) return; // fit failed, abandon this track
430	for (Int_t ip=0; ip<nPtSSD;ip++) {
431	fFitter->GetResiduals(resGlo,iPtSSD[ip]);
432	TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSSD[ip]);
433	mcurr->MasterToLocalVect(resGlo,resLoc);
434	Int_t index=modIdSSD[ip]-kNSPDmods-kNSDDmods;
435	fHistSSDResidX[index]->Fill(track->Pt(),resLoc[0]);
436	fHistSSDResidZ[index]->Fill(track->Pt(),resLoc[2]);
437	}
438	}
439	}
440	//______________________________________________________________________________
441	void AliAnalysisTaskITSAlignQA::Terminate(Option_t /option*/)
442	{
443	// Terminate analysis
444	fOutput = dynamic_cast<TList*> (GetOutputData(1));
445	if (!fOutput) {
446	printf("ERROR: fOutput not available\n");
447	return;
448	}
449
450	fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
451	printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
452	return;
453	}
454
455
456	//___________________________________________________________________________
457	void AliAnalysisTaskITSAlignQA::LoadGeometryFromOCDB(){
458	//method to get the gGeomanager
459	// it is called at the CreatedOutputObject stage
460	// to comply with the CAF environment
461	AliInfo("Loading geometry");
462
463	AliCDBManager *man = AliCDBManager::Instance();
464	man->SetDefaultStorage(fOCDBLocation.Data());
465	man->SetRun(fRunNb);
466	AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
467	AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
468	AliGeomManager::GetNalignable("ITS");
469	AliGeomManager::ApplyAlignObjsFromCDB("ITS");
470	}
471
472
473
474
475
476