PWG1/AliComparisonRes.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliComparisonRes class. It keeps information from
   3 // comparison of reconstructed and MC particle tracks. In addtion,
   4 // it keeps selection cuts used during comparison. The comparison
   5 // information is stored in the ROOT histograms. Analysis of these
   6 // histograms can be done by using Analyse() class function. The result of
   7 // the analysis (histograms/graphs) are stored in the folder which is
   8 // a data member of AliComparisonRes.
   9 //
  10 // Author: J.Otwinowski 04/02/2008
  11 //------------------------------------------------------------------------------
  12
  13 /*
  14
  15   // after running comparison task, read the file, and get component
  16   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
  17   LoadMyLibs();
  18
  19   TFile f("Output.root");
  20   //AliComparisonRes * compObj = (AliComparisonRes*)f.Get("AliComparisonRes");
  21   AliComparisonRes * compObj = (AliComparisonRes*)cOutput->FindObject("AliComparisonRes");
  22
  23   // analyse comparison data
  24   compObj->Analyse();
  25
  26   // the output histograms/graphs will be stored in the folder "folderRes"
  27   compObj->GetAnalysisFolder()->ls("*");
  28
  29   // user can save whole comparison object (or only folder with anlysed histograms)
  30   // in the seperate output file (e.g.)
  31   TFile fout("Analysed_Res.root","recreate");
  32   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  33   fout.Close();
  34
  35 */
  36
  37 #include "TCanvas.h"
  38 #include "TH1.h"
  39 #include "TAxis.h"
  40
  41 #include "AliComparisonRes.h"
  42 #include "AliESDRecInfo.h"
  43 #include "AliESDVertex.h"
  44 #include "AliESDtrack.h"
  45 #include "AliLog.h"
  46 #include "AliMCInfo.h"
  47 #include "AliMCInfoCuts.h"
  48 #include "AliRecInfoCuts.h"
  49 #include "AliTracker.h"
  50 #include "AliTreeDraw.h"
  51
  52 using namespace std;
  53
  54 ClassImp(AliComparisonRes)
  55
  56 //_____________________________________________________________________________
  57 AliComparisonRes::AliComparisonRes():
  58   AliComparisonObject("AliComparisonRes"),
  59   fResolHisto(0),
  60   fPullHisto(0),
  61
  62   // Cuts
  63   fCutsRC(0),
  64   fCutsMC(0),
  65
  66   // histogram folder
  67   fAnalysisFolder(0)
  68 {
  69   //Init();
  70 }
  71
  72 //_____________________________________________________________________________
  73 AliComparisonRes::AliComparisonRes(Char_t* name="AliComparisonRes", Char_t* title="AliComparisonRes",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
  74 //AliComparisonRes::AliComparisonRes(Char_t* name, Char_t* title,Int_t analysisMode,Bool_t hptGenerator):
  75   AliComparisonObject(name,title),
  76   fResolHisto(0),
  77   fPullHisto(0),
  78
  79   // Cuts
  80   fCutsRC(0),
  81   fCutsMC(0),
  82
  83   // histogram folder
  84   fAnalysisFolder(0)
  85 {
  86   // named constructor
  87   //
  88   SetAnalysisMode(analysisMode);
  89   SetHptGenerator(hptGenerator);
  90
  91   Init();
  92 }
  93
  94 //_____________________________________________________________________________
  95 AliComparisonRes::~AliComparisonRes()
  96 {
  97   // destructor
  98
  99   if(fResolHisto) delete fResolHisto; fResolHisto=0;
 100   if(fPullHisto)  delete fPullHisto;  fPullHisto=0;
 101
 102   if(fCutsRC) delete fCutsRC; fCutsRC=0;
 103   if(fCutsMC) delete fCutsMC; fCutsMC=0;
 104
 105   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 106 }
 107
 108 //_____________________________________________________________________________
 109 void AliComparisonRes::Init(){
 110
 111   //
 112   // histogram bining
 113   //
 114
 115   // set pt bins
 116    Int_t nPtBins = 31;
 117    Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
 118   Double_t ptMin = 0., ptMax = 10.;
 119
 120
 121   if(IsHptGenerator() == kTRUE) {
 122     nPtBins = 100;
 123     ptMin = 0.; ptMax = 100.;
 124   }
 125
 126   // res_y:res_z:res_phi,res_lambda:res_1pt:y:z:eta:phi:pt
 127   Int_t binsResolHisto[10]={100,100,100,100,100,50,100,30,144,nPtBins};
 128   Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2,-0.03,-20.,-1.5,0.,ptMin};
 129   Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, 0.03, 20., 1.5,2.*TMath::Pi(), ptMax};
 130
 131
 132   fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_1pt:y:z:eta:phi:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
 133   if(!IsHptGenerator()) fResolHisto->SetBinEdges(9,binsPt);
 134
 135   fResolHisto->GetAxis(0)->SetTitle("res_y (cm)");
 136   fResolHisto->GetAxis(1)->SetTitle("res_z (cm)");
 137   fResolHisto->GetAxis(2)->SetTitle("res_phi (rad)");
 138   fResolHisto->GetAxis(3)->SetTitle("res_lambda (rad)");
 139   fResolHisto->GetAxis(4)->SetTitle("res_1pt (GeV/c)^{-1}");
 140   fResolHisto->GetAxis(5)->SetTitle("y (cm)");
 141   fResolHisto->GetAxis(6)->SetTitle("z (cm)");
 142   fResolHisto->GetAxis(7)->SetTitle("eta");
 143   fResolHisto->GetAxis(8)->SetTitle("phi (rad)");
 144   fResolHisto->GetAxis(9)->SetTitle("pt (GeV/c)");
 145   fResolHisto->Sumw2();
 146
 147   //pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt
 148   Int_t binsPullHisto[10]={100,100,100,100,100,50,100,30,144,nPtBins};
 149   Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,-0.03, -20.,-1.5, 0., ptMin};
 150   Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., 0.03, 20., 1.5, 2.*TMath::Pi(),ptMax};
 151
 152   fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto);
 153   if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,binsPt);
 154
 155   fPullHisto->GetAxis(0)->SetTitle("res_y (cm)");
 156   fPullHisto->GetAxis(1)->SetTitle("res_z (cm)");
 157   fPullHisto->GetAxis(2)->SetTitle("res_phi (rad)");
 158   fPullHisto->GetAxis(3)->SetTitle("res_lambda (rad)");
 159   fPullHisto->GetAxis(4)->SetTitle("res_1pt (GeV/c)^{-1}");
 160   fPullHisto->GetAxis(5)->SetTitle("y (rad)");
 161   fPullHisto->GetAxis(6)->SetTitle("z (rad)");
 162   fPullHisto->GetAxis(7)->SetTitle("eta");
 163   fPullHisto->GetAxis(8)->SetTitle("phi (rad)");
 164   fPullHisto->GetAxis(9)->SetTitle("pt (GeV/c)");
 165   fPullHisto->Sumw2();
 166
 167   // Init cuts
 168   if(!fCutsMC)
 169     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 170   if(!fCutsRC)
 171     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 172
 173   // init folder
 174   fAnalysisFolder = CreateFolder("folderRes","Analysis Resolution Folder");
 175 }
 176
 177 //_____________________________________________________________________________
 178 void AliComparisonRes::ProcessTPC(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 179 {
 180   // Fill resolution comparison information
 181   AliExternalTrackParam *track = 0;
 182   Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]
 183   Double_t kMaxD      = 123456.0; // max distance
 184   AliESDVertex vertexMC;  // MC primary vertex
 185
 186   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 187
 188   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 189   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 190   //Float_t delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC;
 191
 192   Float_t mceta =  infoMC->GetParticle().Eta();
 193   Float_t mcphi =  infoMC->GetParticle().Phi();
 194   if(mcphi<0) mcphi += 2.*TMath::Pi();
 195   Float_t mcpt = infoMC->GetParticle().Pt();
 196
 197   // distance to Prim. vertex
 198   const Double_t* dv = infoMC->GetVDist();
 199   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
 200
 201   // Only 5 charged particle species (e,mu,pi,K,p)
 202   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 203   if (!isPrim) return;
 204   //if (infoRC->GetStatus(1)!=3) return; // TPC refit
 205   if (!infoRC->GetESDtrack()) return;
 206   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 207
 208   // calculate and set prim. vertex
 209   vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] );
 210   vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] );
 211   vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] );
 212
 213   // calculate track parameters at vertex
 214   const AliExternalTrackParam *innerTPC =  0;
 215   if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0)
 216   {
 217     if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
 218     {
 219       Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov);
 220
 221       // Fill parametrisation histograms (only TPC track)
 222       if(bDCAStatus)
 223       {
 224           // select primaries
 225           if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 226           {
 227
 228             deltaYTPC= track->GetY()-infoMC->GetParticle().Vy();
 229             deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz();
 230             deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
 231             deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
 232             delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 233
 234             pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2());
 235             pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 236             pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 237             pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 238             pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 239
 240             Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
 241             fResolHisto->Fill(vResolHisto);
 242
 243             Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
 244             fPullHisto->Fill(vPullHisto);
 245           }
 246
 247           /*
 248             delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2);
 249             deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt);
 250             deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt);
 251             deltaPhi1PtTPC = deltaPhiTPC   / (0.1+1/mcpt);
 252             deltaTheta1PtTPC = deltaLambdaTPC / (0.1+1/mcpt);
 253
 254             fPhiEtaPtTPC->Fill(TMath::ATan2(innerTPC->Py(),innerTPC->Px()), innerTPC->Eta(), innerTPC->Pt());
 255
 256             f1Pt2ResolS1PtTPC->Fill(s1mcpt,delta1Pt2TPC);
 257             fYResolS1PtTPC->Fill(s1mcpt,deltaY1PtTPC);
 258             fZResolS1PtTPC->Fill(s1mcpt,deltaZ1PtTPC);
 259             fPhiResolS1PtTPC->Fill(s1mcpt,deltaPhi1PtTPC);
 260             fThetaResolS1PtTPC->Fill(s1mcpt,deltaTheta1PtTPC);
 261
 262             fPtResolPtTPC->Fill(mcpt,deltaPtTPC);
 263             fPtPullPtTPC->Fill(mcpt,pullPtTPC);
 264             fPhiResolEtaTPC->Fill(eta,deltaPhiTPC);
 265             fPhiPullEtaTPC->Fill(eta,pullPhiTPC);
 266             fLambdaResolEtaTPC->Fill(eta,deltaLambdaTPC);
 267             fLambdaPullEtaTPC->Fill(eta,pullLambdaTPC);
 268           */
 269         }
 270     delete track;
 271     }
 272   }
 273 }
 274
 275 //_____________________________________________________________________________
 276 void AliComparisonRes::ProcessTPCITS(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 277 {
 278   Int_t clusterITS[200];
 279   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 280
 281   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 282   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 283
 284   Float_t mceta =  infoMC->GetParticle().Eta();
 285   Float_t mcphi =  infoMC->GetParticle().Phi();
 286   if(mcphi<0) mcphi += 2.*TMath::Pi();
 287   Float_t mcpt = infoMC->GetParticle().Pt();
 288
 289   // distance to Prim. vertex
 290   const Double_t* dv = infoMC->GetVDist();
 291   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
 292
 293   // Only 5 charged particle species (e,mu,pi,K,p)
 294   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 295   if (!isPrim) return;
 296   if (infoRC->GetStatus(1)!=3) return; // TPC refit
 297   if (!infoRC->GetESDtrack()) return;
 298   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 299
 300   infoRC->GetESDtrack()->GetImpactParameters(dca,cov);
 301
 302   // select primaries
 303   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 304   {
 305     deltaYTPC= infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy();
 306     deltaZTPC = infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz();
 307     deltaLambdaTPC = TMath::ATan2(infoRC->GetESDtrack()->Pz(),infoRC->GetESDtrack()->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
 308     deltaPhiTPC = TMath::ATan2(infoRC->GetESDtrack()->Py(),infoRC->GetESDtrack()->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
 309     delta1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt)*mcpt;
 310
 311     pullYTPC= (infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaY2());
 312     pullZTPC = (infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaZ2());
 313     pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaTgl2());
 314     pullPhiTPC = deltaPhiTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaSnp2());
 315     pull1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2());
 316
 317     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
 318     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
 319
 320     // TPC and ITS clusters in the system
 321     if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>fCutsRC->GetMinNClustersITS())
 322     {
 323       fResolHisto->Fill(vResolHisto);
 324       fPullHisto->Fill(vPullHisto);
 325     }
 326   }
 327 }
 328
 329 //_____________________________________________________________________________
 330 void AliComparisonRes::ProcessConstrained(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 331 {
 332   // Fill resolution comparison information (constarained parameters)
 333   //
 334   AliExternalTrackParam *track = 0;
 335   Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]
 336   Double_t kMaxD      = 123456.0; // max distance
 337   AliESDVertex vertexMC;  // MC primary vertex
 338
 339   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 340
 341   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 342   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 343
 344   Float_t mceta =  infoMC->GetParticle().Eta();
 345   Float_t mcphi =  infoMC->GetParticle().Phi();
 346   if(mcphi<0) mcphi += 2.*TMath::Pi();
 347   Float_t mcpt = infoMC->GetParticle().Pt();
 348
 349   // distance to Prim. vertex
 350   const Double_t* dv = infoMC->GetVDist();
 351   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
 352
 353   // calculate and set prim. vertex
 354   vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] );
 355   vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] );
 356   vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] );
 357
 358   // Check selection cuts
 359   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 360   if (!isPrim) return;
 361   if (infoRC->GetStatus(1)!=3) return; // TPC refit
 362   if (!infoRC->GetESDtrack()) return;
 363   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 364
 365   // constrained parameters resolution
 366   const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam();
 367   if(!cparam) return;
 368
 369   if ((track = new AliExternalTrackParam(*cparam)) != 0)
 370   {
 371     Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov);
 372     if(bDCAStatus) {
 373       if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 374       {
 375         deltaYTPC= track->GetY()-infoMC->GetParticle().Vy();
 376         deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz();
 377         deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
 378         deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
 379         delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 380
 381         pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2());
 382         pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 383         pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 384         pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 385         pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 386
 387         Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
 388         fResolHisto->Fill(vResolHisto);
 389
 390         Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
 391         fPullHisto->Fill(vPullHisto);
 392       }
 393     }
 394   delete track;
 395   }
 396 }
 397
 398 //_____________________________________________________________________________
 399 void AliComparisonRes::Exec(AliMCInfo* const infoMC, AliESDRecInfo* const infoRC){
 400
 401   // Process comparison information
 402   if(GetAnalysisMode() == 0) ProcessTPC(infoMC,infoRC);
 403   else if(GetAnalysisMode() == 1) ProcessTPCITS(infoMC,infoRC);
 404   else if(GetAnalysisMode() == 2) ProcessConstrained(infoMC,infoRC);
 405   else {
 406     printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 407     return;
 408   }
 409 }
 410
 411 //_____________________________________________________________________________
 412 TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){
 413   // Create resolution histograms
 414
 415   TH1F *hisr, *hism;
 416   if (!gPad) new TCanvas;
 417   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
 418   if (type) return hism;
 419   else
 420     return hisr;
 421 }
 422
 423 //_____________________________________________________________________________
 424 void AliComparisonRes::Analyse(){
 425   // Analyse comparison information and store output histograms
 426   // in the folder "folderRes"
 427   //
 428   TH1::AddDirectory(kFALSE);
 429   TH1F *h=0;
 430   TH2F *h2D=0;
 431   TObjArray *aFolderObj = new TObjArray;
 432
 433   // write results in the folder
 434   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
 435   c->cd();
 436
 437   char name[256];
 438   char res_title[256] = {"res_y:res_z:res_lambda:res_phi:res_1pt:y:z:eta:phi:pt"} ;
 439   char pull_title[256] = {"pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt"};
 440
 441   for(Int_t i=0; i<5; i++)
 442   {
 443     for(Int_t j=5; j<10; j++)
 444     {
 445       if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.2,10.); // pt threshold
 446       if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window
 447
 448       h2D = (TH2F*)fResolHisto->Projection(i,j);
 449       h = AliComparisonRes::MakeResol(h2D,1,0);
 450       sprintf(name,"h_res_%d_vs_%d",i,j);
 451       h->SetName(name);
 452       h->SetTitle(res_title);
 453
 454       aFolderObj->Add(h);
 455
 456       h = AliComparisonRes::MakeResol(h2D,1,1);
 457       sprintf(name,"h_mean_res_%d_vs_%d",i,j);
 458       h->SetName(name);
 459       h->SetTitle(res_title);
 460
 461       aFolderObj->Add(h);
 462
 463       if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
 464       if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
 465
 466       //
 467       if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.2,10.);
 468       if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-0.9,0.9);
 469
 470       h2D = (TH2F*)fPullHisto->Projection(i,j);
 471       h = AliComparisonRes::MakeResol(h2D,1,0);
 472       sprintf(name,"h_pull_%d_vs_%d",i,j);
 473       h->SetName(name);
 474       h->SetTitle(pull_title);
 475
 476       aFolderObj->Add(h);
 477
 478       h = AliComparisonRes::MakeResol(h2D,1,1);
 479       sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
 480       h->SetName(name);
 481       h->SetTitle(pull_title);
 482
 483       aFolderObj->Add(h);
 484
 485       if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.);
 486       if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
 487     }
 488   }
 489
 490   // export objects to analysis folder
 491   fAnalysisFolder = ExportToFolder(aFolderObj);
 492
 493   // delete only TObjArray
 494   if(aFolderObj) delete aFolderObj;
 495 }
 496
 497 //_____________________________________________________________________________
 498 TFolder* AliComparisonRes::ExportToFolder(TObjArray * array)
 499 {
 500   // recreate folder avery time and export objects to new one
 501   //
 502   AliComparisonRes * comp=this;
 503   TFolder *folder = comp->GetAnalysisFolder();
 504
 505   TString name, title;
 506   TFolder *newFolder = 0;
 507   Int_t i = 0;
 508   Int_t size = array->GetSize();
 509
 510   if(folder) {
 511      // get name and title from old folder
 512      name = folder->GetName();
 513      title = folder->GetTitle();
 514
 515          // delete old one
 516      delete folder;
 517
 518          // create new one
 519      newFolder = CreateFolder(name.Data(),title.Data());
 520      newFolder->SetOwner();
 521
 522          // add objects to folder
 523      while(i < size) {
 524            newFolder->Add(array->At(i));
 525            i++;
 526          }
 527   }
 528
 529 return newFolder;
 530 }
 531
 532 //_____________________________________________________________________________
 533 Long64_t AliComparisonRes::Merge(TCollection* const list)
 534 {
 535   // Merge list of objects (needed by PROOF)
 536
 537   if (!list)
 538   return 0;
 539
 540   if (list->IsEmpty())
 541   return 1;
 542
 543   TIterator* iter = list->MakeIterator();
 544   TObject* obj = 0;
 545
 546   // collection of generated histograms
 547   Int_t count=0;
 548   while((obj = iter->Next()) != 0)
 549   {
 550   AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj);
 551   if (entry == 0) continue;
 552
 553   fResolHisto->Add(entry->fResolHisto);
 554   fPullHisto->Add(entry->fPullHisto);
 555
 556   count++;
 557   }
 558
 559 return count;
 560 }
 561
 562 //_____________________________________________________________________________
 563 TFolder* AliComparisonRes::CreateFolder(TString name,TString title) {
 564 // create folder for analysed histograms
 565 //
 566 TFolder *folder = 0;
 567   folder = new TFolder(name.Data(),title.Data());
 568
 569   return folder;
 570 }