PWG1/AliPerformanceRes.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliPerformanceRes 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 AliPerformanceRes.
   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   //AliPerformanceRes * compObj = (AliPerformanceRes*)f.Get("AliPerformanceRes");
  21   AliPerformanceRes * compObj = (AliPerformanceRes*)cOutput->FindObject("AliPerformanceRes");
  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 "TH2.h"
  40 #include "TAxis.h"
  41 #include "TF1.h"
  42
  43 #include "AliPerformanceRes.h"
  44 #include "AliESDEvent.h"
  45 #include "AliESDVertex.h"
  46 #include "AliESDtrack.h"
  47 #include "AliLog.h"
  48 #include "AliMCEvent.h"
  49 #include "AliMCParticle.h"
  50 #include "AliHeader.h"
  51 #include "AliGenEventHeader.h"
  52 #include "AliStack.h"
  53 #include "AliMCInfoCuts.h"
  54 #include "AliRecInfoCuts.h"
  55 #include "AliTracker.h"
  56 #include "AliTreeDraw.h"
  57
  58 using namespace std;
  59
  60 ClassImp(AliPerformanceRes)
  61
  62 //_____________________________________________________________________________
  63 AliPerformanceRes::AliPerformanceRes():
  64   AliPerformanceObject("AliPerformanceRes"),
  65   fResolHisto(0),
  66   fPullHisto(0),
  67
  68   // Cuts
  69   fCutsRC(0),
  70   fCutsMC(0),
  71
  72   // histogram folder
  73   fAnalysisFolder(0)
  74 {
  75   Init();
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliPerformanceRes::AliPerformanceRes(Char_t* name="AliPerformanceRes", Char_t* title="AliPerformanceRes",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
  80   AliPerformanceObject(name,title),
  81   fResolHisto(0),
  82   fPullHisto(0),
  83
  84   // Cuts
  85   fCutsRC(0),
  86   fCutsMC(0),
  87
  88   // histogram folder
  89   fAnalysisFolder(0)
  90 {
  91   // named constructor
  92   //
  93   SetAnalysisMode(analysisMode);
  94   SetHptGenerator(hptGenerator);
  95
  96   Init();
  97 }
  98
  99 //_____________________________________________________________________________
 100 AliPerformanceRes::~AliPerformanceRes()
 101 {
 102   // destructor
 103
 104   if(fResolHisto) delete fResolHisto; fResolHisto=0;
 105   if(fPullHisto)  delete fPullHisto;  fPullHisto=0;
 106
 107   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 108 }
 109
 110 //_____________________________________________________________________________
 111 void AliPerformanceRes::Init(){
 112
 113   //
 114   // histogram bining
 115   //
 116
 117   // set pt bins
 118   Int_t nPtBins = 50;
 119   Double_t ptMin = 1.e-2, ptMax = 10.;
 120
 121   Double_t *binsPt = 0;
 122   if (IsHptGenerator())  {
 123     nPtBins = 100; ptMax = 100.;
 124     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 125   } else {
 126     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 127   }
 128
 129   /*
 130   Int_t nPtBins = 31;
 131   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.};
 132   Double_t ptMin = 0., ptMax = 10.;
 133
 134   if(IsHptGenerator() == kTRUE) {
 135     nPtBins = 100;
 136     ptMin = 0.; ptMax = 100.;
 137   }
 138   */
 139
 140   Double_t yMin = -0.02, yMax = 0.02;
 141   Double_t zMin = -12.0, zMax = 12.0;
 142   if(GetAnalysisMode() == 3) { // TrackRef coordinate system
 143     yMin = -100.; yMax = 100.;
 144     zMin = -100.; zMax = 100.;
 145   }
 146
 147   // res_y:res_z:res_phi,res_lambda:res_1pt:y:z:eta:phi:pt
 148   Int_t binsResolHisto[10]={100,100,100,100,100,25,50,30,144,nPtBins};
 149   Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2, yMin, zMin,-1.5,0.,ptMin};
 150   Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, yMax, zMax, 1.5,2.*TMath::Pi(), ptMax};
 151
 152   fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_1pt:y:z:eta:phi:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
 153   fResolHisto->SetBinEdges(9,binsPt);
 154
 155   fResolHisto->GetAxis(0)->SetTitle("y-y_{mc} (cm)");
 156   fResolHisto->GetAxis(1)->SetTitle("z-z_{mc} (cm)");
 157   fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{mc} (rad)");
 158   fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{mc} (rad)");
 159   fResolHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)");
 160   fResolHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
 161   fResolHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
 162   fResolHisto->GetAxis(7)->SetTitle("#eta_{mc}");
 163   fResolHisto->GetAxis(8)->SetTitle("#phi_{mc} (rad)");
 164   fResolHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
 165   fResolHisto->Sumw2();
 166
 167   //pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt
 168   Int_t binsPullHisto[10]={100,100,100,100,100,50,50,30,144,nPtBins};
 169   Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1.5, 0., ptMin};
 170   Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1.5, 2.*TMath::Pi(),ptMax};
 171
 172   fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto);
 173   if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,binsPt);
 174
 175   fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
 176   fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
 177   fPullHisto->GetAxis(2)->SetTitle("(#phi-#phi_{mc})/#sigma");
 178   fPullHisto->GetAxis(3)->SetTitle("(#lambda-#lambda_{mc})/#sigma");
 179   fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma");
 180   fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
 181   fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
 182   fPullHisto->GetAxis(7)->SetTitle("#eta_{mc}");
 183   fPullHisto->GetAxis(8)->SetTitle("#phi_{mc} (rad)");
 184   fPullHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
 185   fPullHisto->Sumw2();
 186
 187   // Init cuts
 188   if(!fCutsMC)
 189     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 190   if(!fCutsRC)
 191     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 192
 193   // init folder
 194   fAnalysisFolder = CreateFolder("folderRes","Analysis Resolution Folder");
 195 }
 196
 197 //_____________________________________________________________________________
 198 void AliPerformanceRes::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
 199 {
 200   if(!esdTrack) return;
 201
 202   // Fill TPC only resolution comparison information
 203   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
 204   if(!track) return;
 205
 206   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 207   esdTrack->GetImpactParametersTPC(dca,cov);
 208
 209   //
 210   // Fill rec vs MC information
 211   //
 212   if(!stack) return;
 213   Int_t label = TMath::Abs(esdTrack->GetLabel());
 214   if(label == 0) return;
 215
 216   TParticle* particle = stack->Particle(label);
 217   if(!particle) return;
 218   if(!particle->GetPDG()) return;
 219   if(particle->GetPDG()->Charge()==0) return;
 220   //printf("charge %d \n",particle->GetPDG()->Charge());
 221
 222   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 223   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 224
 225   Float_t mceta =  particle->Eta();
 226   Float_t mcphi =  particle->Phi();
 227   if(mcphi<0) mcphi += 2.*TMath::Pi();
 228   Float_t mcpt = particle->Pt();
 229
 230   // nb. TPC clusters cut
 231   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 232
 233   // select primaries
 234   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 235   {
 236     deltaYTPC= track->GetY()-particle->Vy();
 237     deltaZTPC = track->GetZ()-particle->Vz();
 238     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 239     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
 240     delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 241
 242     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
 243     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 244
 245     Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 246     Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 247     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
 248     pullPhiTPC = deltaPhiTPC / sigma_phi;
 249
 250     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 251     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 252     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 253     else pull1PtTPC = 0.;
 254
 255     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 256     fResolHisto->Fill(vResolHisto);
 257
 258     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 259     fPullHisto->Fill(vPullHisto);
 260   }
 261 }
 262
 263 //_____________________________________________________________________________
 264 void AliPerformanceRes::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack)
 265 {
 266   // Fill resolution comparison information (TPC+ITS)
 267   if(!esdTrack) return;
 268
 269   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 270   esdTrack->GetImpactParameters(dca,cov);
 271
 272   //
 273   // Fill rec vs MC information
 274   //
 275   if(!stack) return;
 276
 277   Int_t label = TMath::Abs(esdTrack->GetLabel());
 278   if(label == 0) return;
 279
 280   TParticle* particle = stack->Particle(label);
 281   if(!particle) return;
 282   if(particle->GetPDG()->Charge()==0) return;
 283   //printf("charge %d \n",particle->GetPDG()->Charge());
 284
 285   Float_t mceta =  particle->Eta();
 286   Float_t mcphi =  particle->Phi();
 287   if(mcphi<0) mcphi += 2.*TMath::Pi();
 288   Float_t mcpt = particle->Pt();
 289
 290   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 291   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 292   Int_t clusterITS[200];
 293   if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return;  // min. nb. ITS clusters
 294
 295   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 296   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 297
 298   // select primaries
 299   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 300   {
 301     deltaYTPC= esdTrack->GetY()-particle->Vy();
 302     deltaZTPC = esdTrack->GetZ()-particle->Vz();
 303     deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 304     deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
 305     delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
 306
 307     pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
 308     pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
 309
 310     Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2());
 311     Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
 312     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
 313     pullPhiTPC = deltaPhiTPC / sigma_phi;
 314
 315     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
 316     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2());
 317     if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
 318     else pull1PtTPC = 0.;
 319
 320     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 321     fResolHisto->Fill(vResolHisto);
 322
 323     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 324     fPullHisto->Fill(vPullHisto);
 325   }
 326 }
 327
 328 //_____________________________________________________________________________
 329 void AliPerformanceRes::ProcessConstrained(AliStack* const stack, AliESDtrack *const esdTrack)
 330 {
 331   // Fill resolution comparison information (constarained parameters)
 332   //
 333   if(!esdTrack) return;
 334
 335   const AliExternalTrackParam * track = esdTrack->GetConstrainedParam();
 336   if(!track) return;
 337
 338   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 339   esdTrack->GetImpactParameters(dca,cov);
 340
 341   //
 342   // Fill rec vs MC information
 343   //
 344   if(!stack) return;
 345
 346   Int_t label = TMath::Abs(esdTrack->GetLabel());
 347   if(label == 0) return;
 348
 349   TParticle* particle = stack->Particle(label);
 350   if(!particle) return;
 351   if(particle->GetPDG()->Charge()==0) return;
 352   //printf("charge %d \n",particle->GetPDG()->Charge());
 353
 354   Float_t mceta =  particle->Eta();
 355   Float_t mcphi =  particle->Phi();
 356   if(mcphi<0) mcphi += 2.*TMath::Pi();
 357   Float_t mcpt = particle->Pt();
 358
 359   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 360   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 361
 362   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 363   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 364
 365   // select primaries
 366   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 367   {
 368     deltaYTPC= track->GetY()-particle->Vy();
 369     deltaZTPC = track->GetZ()-particle->Vz();
 370     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 371     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
 372     delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 373
 374     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
 375     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 376
 377     Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 378     Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 379     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
 380     pullPhiTPC = deltaPhiTPC / sigma_phi;
 381
 382     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 383     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 384
 385     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 386     else pull1PtTPC = 0.;
 387
 388     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 389     fResolHisto->Fill(vResolHisto);
 390
 391     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 392     fPullHisto->Fill(vPullHisto);
 393   }
 394 }
 395
 396 //_____________________________________________________________________________
 397 void AliPerformanceRes::ProcessInnerTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack)
 398 {
 399   //
 400   // Fill resolution comparison information (inner params at TPC reference point)
 401   //
 402   if(!esdTrack) return;
 403
 404   const AliExternalTrackParam * innerParam = esdTrack->GetInnerParam();
 405   if(!innerParam) return;
 406
 407   // create new AliExternalTrackParam
 408   AliExternalTrackParam *track = new AliExternalTrackParam(*innerParam);
 409   if(!track) return;
 410
 411   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 412   esdTrack->GetImpactParameters(dca,cov);
 413
 414   //
 415   // Fill rec vs MC information
 416   //
 417   if(!mcEvent) return;
 418
 419   Int_t label = TMath::Abs(esdTrack->GetLabel());
 420   if(label == 0) return;
 421   AliMCParticle *mcParticle = mcEvent->GetTrack(label);
 422   if(!mcParticle) return;
 423
 424   // get the first TPC track reference
 425   AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle);
 426   if(!ref0) return;
 427
 428   // calculate alpha angle
 429   Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
 430   Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 431   //if (alpha<0) alpha += TMath::TwoPi();
 432
 433   // rotate inner track to local coordinate system
 434   // and propagate to the radius of the first track referenco of TPC
 435   Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
 436   Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
 437   if(!isOK) return;
 438
 439   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
 440   Float_t mcphi =  ref0->Phi();
 441   if(mcphi<0) mcphi += 2.*TMath::Pi();
 442   Float_t mcpt = ref0->Pt();
 443
 444   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 445   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 446
 447   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 448   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 449
 450   // select primaries
 451   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 452   {
 453     //deltaYTPC= track->GetY()-ref0->Y();
 454     deltaYTPC= track->GetY(); // it corresponds to deltaY after alpha rotation
 455     deltaZTPC = track->GetZ()-ref0->Z();
 456     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
 457     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
 458     if(mcpt) delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 459     else delta1PtTPC = 0.;
 460
 461     //pullYTPC= (track->GetY()-ref0->Y()) / TMath::Sqrt(track->GetSigmaY2());
 462     pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2()); //it corresponds to deltaY after alpha rotation
 463     pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
 464
 465     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 466     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 467
 468     Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 469     Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 470     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
 471     pullPhiTPC = deltaPhiTPC / sigma_phi;
 472
 473     if(mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 474     else pull1PtTPC = 0.;
 475
 476     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,ref0->Y(),ref0->Z(),mceta,mcphi,mcpt};
 477     fResolHisto->Fill(vResolHisto);
 478
 479     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mceta,mcphi,mcpt};
 480     fPullHisto->Fill(vPullHisto);
 481   }
 482
 483   if(track) delete track;
 484 }
 485
 486 //_____________________________________________________________________________
 487 AliTrackReference * AliPerformanceRes::GetFirstTPCTrackRef(AliMCParticle *mcParticle)
 488 {
 489   // return first TPC track reference
 490   if(!mcParticle) return 0;
 491
 492   Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
 493   AliTrackReference *ref0 = 0;
 494   for (Int_t iref = 0; iref < nTrackRef; iref++) {
 495     ref0 = mcParticle->GetTrackReference(iref);
 496     if(ref0 && (ref0->DetectorId()==AliTrackReference::kTPC))
 497       break;
 498   }
 499
 500 return ref0;
 501 }
 502
 503 //_____________________________________________________________________________
 504 void AliPerformanceRes::Exec(AliMCEvent* const mcEvent, AliESDEvent* const esdEvent, const Bool_t bUseMC)
 505 {
 506   // Process comparison information
 507   //
 508   if(!esdEvent)
 509   {
 510       AliDebug(AliLog::kError, "esdEvent not available");
 511       return;
 512   }
 513   AliHeader* header = 0;
 514   AliGenEventHeader* genHeader = 0;
 515   AliStack* stack = 0;
 516   TArrayF vtxMC(3);
 517
 518   if(bUseMC)
 519   {
 520     if(!mcEvent) {
 521       AliDebug(AliLog::kError, "mcEvent not available");
 522       return;
 523     }
 524
 525     // get MC event header
 526     header = mcEvent->Header();
 527     if (!header) {
 528       AliDebug(AliLog::kError, "Header not available");
 529       return;
 530     }
 531     // MC particle stack
 532     stack = mcEvent->Stack();
 533     if (!stack) {
 534       AliDebug(AliLog::kError, "Stack not available");
 535       return;
 536     }
 537
 538     // get MC vertex
 539     genHeader = header->GenEventHeader();
 540     if (!genHeader) {
 541       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 542       return;
 543     }
 544     genHeader->PrimaryVertex(vtxMC);
 545
 546   } // end bUseMC
 547
 548   //  Process events
 549   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 550   {
 551     AliESDtrack*track = esdEvent->GetTrack(iTrack);
 552     if(!track) continue;
 553
 554     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
 555     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
 556     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
 557     else if(GetAnalysisMode() == 3) ProcessInnerTPC(mcEvent,track);
 558     else {
 559       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 560       return;
 561     }
 562   }
 563 }
 564
 565 //_____________________________________________________________________________
 566 TH1F* AliPerformanceRes::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
 567   // Create resolution histograms
 568
 569   TH1F *hisr, *hism;
 570   if (!gPad) new TCanvas;
 571   hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
 572   if (type) return hism;
 573   else
 574     return hisr;
 575 }
 576
 577 //_____________________________________________________________________________
 578 void AliPerformanceRes::Analyse(){
 579   // Analyse comparison information and store output histograms
 580   // in the folder "folderRes"
 581   //
 582   TH1::AddDirectory(kFALSE);
 583   TH1F *h=0;
 584   TH2F *h2D=0;
 585   TObjArray *aFolderObj = new TObjArray;
 586
 587   // write results in the folder
 588   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
 589   c->cd();
 590
 591   char name[256];
 592   char title[256];
 593
 594   for(Int_t i=0; i<5; i++)
 595   {
 596     for(Int_t j=5; j<10; j++)
 597     {
 598       if(j!=7) fResolHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window
 599       fResolHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
 600       if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
 601
 602       h2D = (TH2F*)fResolHisto->Projection(i,j);
 603
 604       h = AliPerformanceRes::MakeResol(h2D,1,0,100);
 605       sprintf(name,"h_res_%d_vs_%d",i,j);
 606
 607       /*
 608       if(i<2) {
 609         h->SetMinimum(0.);
 610         h->SetMaximum(1.);
 611       } else {
 612         h->SetMinimum(0.);
 613         h->SetMaximum(0.2);
 614       }
 615       */
 616       h->SetName(name);
 617
 618       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 619       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
 620       h->GetYaxis()->SetTitle(title);
 621       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 622       h->SetTitle(title);
 623
 624       if(j==9) h->SetBit(TH1::kLogX);
 625       aFolderObj->Add(h);
 626
 627       h = AliPerformanceRes::MakeResol(h2D,1,1,100);
 628       //h = (TH1F*)arr->At(1);
 629       sprintf(name,"h_mean_res_%d_vs_%d",i,j);
 630       /*
 631       h->SetMinimum(-0.05);
 632       h->SetMaximum(0.05);
 633       if(i<2) {
 634         h->SetMinimum(-0.05);
 635         h->SetMaximum(0.05);
 636       } else {
 637         h->SetMinimum(-0.02);
 638         h->SetMaximum(0.02);
 639       }
 640       */
 641       h->SetName(name);
 642
 643       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 644       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
 645       h->GetYaxis()->SetTitle(title);
 646
 647       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 648       h->SetTitle(title);
 649
 650       if(j==9) h->SetBit(TH1::kLogX);
 651       aFolderObj->Add(h);
 652
 653       fResolHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
 654       fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
 655
 656       //
 657       if(j!=7) fPullHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window
 658       fPullHisto->GetAxis(9)->SetRangeUser(0.16,10.);  // pt threshold
 659       if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
 660
 661       h2D = (TH2F*)fPullHisto->Projection(i,j);
 662
 663       h = AliPerformanceRes::MakeResol(h2D,1,0,100);
 664       sprintf(name,"h_pull_%d_vs_%d",i,j);
 665       //h->SetMinimum(0.);
 666       //h->SetMaximum(2.);
 667       h->SetName(name);
 668
 669       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 670       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
 671       h->GetYaxis()->SetTitle(title);
 672       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 673       h->SetTitle(title);
 674
 675       if(j==9) h->SetBit(TH1::kLogX);
 676       aFolderObj->Add(h);
 677
 678       h = AliPerformanceRes::MakeResol(h2D,1,1,100);
 679       sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
 680       //h->SetMinimum(-0.2);
 681       //h->SetMaximum(0.2);
 682       h->SetName(name);
 683
 684       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 685       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
 686       h->GetYaxis()->SetTitle(title);
 687       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 688       h->SetTitle(title);
 689
 690       if(j==9) h->SetBit(TH1::kLogX);
 691       aFolderObj->Add(h);
 692
 693       fPullHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
 694       fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.);
 695     }
 696   }
 697
 698   // export objects to analysis folder
 699   fAnalysisFolder = ExportToFolder(aFolderObj);
 700
 701   // delete only TObjArray
 702   if(aFolderObj) delete aFolderObj;
 703 }
 704
 705 //_____________________________________________________________________________
 706 TFolder* AliPerformanceRes::ExportToFolder(TObjArray * array)
 707 {
 708   // recreate folder avery time and export objects to new one
 709   //
 710   AliPerformanceRes * comp=this;
 711   TFolder *folder = comp->GetAnalysisFolder();
 712
 713   TString name, title;
 714   TFolder *newFolder = 0;
 715   Int_t i = 0;
 716   Int_t size = array->GetSize();
 717
 718   if(folder) {
 719      // get name and title from old folder
 720      name = folder->GetName();
 721      title = folder->GetTitle();
 722
 723          // delete old one
 724      delete folder;
 725
 726          // create new one
 727      newFolder = CreateFolder(name.Data(),title.Data());
 728      newFolder->SetOwner();
 729
 730          // add objects to folder
 731      while(i < size) {
 732            newFolder->Add(array->At(i));
 733            i++;
 734          }
 735   }
 736
 737 return newFolder;
 738 }
 739
 740 //_____________________________________________________________________________
 741 Long64_t AliPerformanceRes::Merge(TCollection* const list)
 742 {
 743   // Merge list of objects (needed by PROOF)
 744
 745   if (!list)
 746   return 0;
 747
 748   if (list->IsEmpty())
 749   return 1;
 750
 751   TIterator* iter = list->MakeIterator();
 752   TObject* obj = 0;
 753
 754   // collection of generated histograms
 755   Int_t count=0;
 756   while((obj = iter->Next()) != 0)
 757   {
 758   AliPerformanceRes* entry = dynamic_cast<AliPerformanceRes*>(obj);
 759   if (entry == 0) continue;
 760
 761   fResolHisto->Add(entry->fResolHisto);
 762   fPullHisto->Add(entry->fPullHisto);
 763
 764   count++;
 765   }
 766
 767 return count;
 768 }
 769
 770 //_____________________________________________________________________________
 771 TFolder* AliPerformanceRes::CreateFolder(TString name,TString title) {
 772 // create folder for analysed histograms
 773 //
 774 TFolder *folder = 0;
 775   folder = new TFolder(name.Data(),title.Data());
 776
 777   return folder;
 778 }