PWG1/TPC/AliPerformanceMatch.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliPerformanceMatch class. It checks the track matching
   3 // quality (residuals, pulls) bewteen tracking detectors TPC-ITS and TPC-TRD
   4 // with TPC as the reference detector. In addition, the ITS and TRD track
   5 // reconstruction efficiency with respect to TPC is calculated.
   6 //
   7 // The matchinig quality parameters are stored in the THnSparse histograms.
   8 // The analysis of these histograms to extract reference information is done in
   9 // the Analyse() function (post-processing).
  10 //
  11 // The histograms with reference information can be exported to the ROOT folder.
  12 //
  13 // Author: J.Otwinowski 17/10/2009
  14 //------------------------------------------------------------------------------
  15
  16 /*
  17
  18   // after running comparison task, read the file, and get component
  19   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
  20   LoadMyLibs();
  21
  22   TFile f("Output.root");
  23   AliPerformanceMatch * compObj = (AliPerformanceMatch*)coutput->FindObject("AliPerformanceMatchTPCTRD");
  24
  25   // analyse comparison data
  26   compObj->Analyse();
  27
  28   // the output histograms/graphs will be stored in the folder "folderMatch"
  29   compObj->GetAnalysisFolder()->ls("*");
  30
  31   // user can save whole comparison object (or only folder with anlysed histograms)
  32   // in the seperate output file (e.g.)
  33   TFile fout("Analysed_Match.root","recreate");
  34   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  35   fout.Close();
  36
  37 */
  38
  39 #include "TCanvas.h"
  40 #include "TH1.h"
  41 #include "TH2.h"
  42 #include "TAxis.h"
  43 #include "TF1.h"
  44
  45 #include "AliPerformanceMatch.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDVertex.h"
  48 #include "AliESDtrack.h"
  49 #include "AliESDfriendTrack.h"
  50 #include "AliESDfriend.h"
  51 #include "AliLog.h"
  52 #include "AliMCEvent.h"
  53 #include "AliMCParticle.h"
  54 #include "AliHeader.h"
  55 #include "AliGenEventHeader.h"
  56 #include "AliStack.h"
  57 #include "AliMCInfoCuts.h"
  58 #include "AliRecInfoCuts.h"
  59 #include "AliTracker.h"
  60 #include "AliTRDtrackV1.h"
  61 #include "AliTreeDraw.h"
  62
  63 using namespace std;
  64
  65 ClassImp(AliPerformanceMatch)
  66
  67 //_____________________________________________________________________________
  68 AliPerformanceMatch::AliPerformanceMatch():
  69   AliPerformanceObject("AliPerformanceMatch"),
  70   fResolHisto(0),
  71   fPullHisto(0),
  72   fTrackingEffHisto(0),
  73
  74   // Cuts
  75   fCutsRC(0),
  76   fCutsMC(0),
  77
  78   // histogram folder
  79   fAnalysisFolder(0)
  80 {
  81   Init();
  82 }
  83
  84 //_____________________________________________________________________________
  85 AliPerformanceMatch::AliPerformanceMatch(Char_t* name="AliPerformanceMatch", Char_t* title="AliPerformanceMatch",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
  86   AliPerformanceObject(name,title),
  87   fResolHisto(0),
  88   fPullHisto(0),
  89   fTrackingEffHisto(0),
  90
  91   // Cuts
  92   fCutsRC(0),
  93   fCutsMC(0),
  94
  95   // histogram folder
  96   fAnalysisFolder(0)
  97 {
  98   // named constructor
  99   //
 100   SetAnalysisMode(analysisMode);
 101   SetHptGenerator(hptGenerator);
 102
 103   Init();
 104 }
 105
 106 //_____________________________________________________________________________
 107 AliPerformanceMatch::~AliPerformanceMatch()
 108 {
 109   // destructor
 110
 111   if(fResolHisto) delete fResolHisto; fResolHisto=0;
 112   if(fPullHisto)  delete fPullHisto;  fPullHisto=0;
 113   if(fTrackingEffHisto) delete fTrackingEffHisto; fTrackingEffHisto = 0x0;
 114
 115   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 116 }
 117
 118 //_____________________________________________________________________________
 119 void AliPerformanceMatch::Init(){
 120
 121   //
 122   // Make performance histogrms
 123   //
 124
 125   // set pt bins
 126   Int_t nPtBins = 50;
 127   Double_t ptMin = 1.e-2, ptMax = 10.;
 128
 129   Double_t *binsPt = 0;
 130   if (IsHptGenerator())  {
 131     nPtBins = 100; ptMax = 100.;
 132     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 133   } else {
 134     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 135   }
 136
 137   Double_t yMin = -0.02, yMax = 0.02;
 138   Double_t zMin = -12.0, zMax = 12.0;
 139   if(GetAnalysisMode() == 0 || GetAnalysisMode() == 1) {
 140     yMin = -100.; yMax = 100.;
 141     zMin = -100.; zMax = 100.;
 142   }
 143
 144   // res_y:res_z:res_phi,res_lambda:res_pt:y:z:eta:phi:pt:isRec
 145   Int_t binsResolHisto[11]={100,100,100,100,100,25,50,90,30,nPtBins,2};
 146   Double_t minResolHisto[11]={-1.,-1.,-0.03,-0.03,-0.2, yMin, zMin, 0., -1.5, ptMin,0};
 147   Double_t maxResolHisto[11]={ 1., 1., 0.03, 0.03, 0.2, yMax, zMax, 2.*TMath::Pi(), 1.5, ptMax,2};
 148
 149   fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_pt:y:z:phi:eta:pt:isRec",11,binsResolHisto,minResolHisto,maxResolHisto);
 150   fResolHisto->SetBinEdges(9,binsPt);
 151
 152   fResolHisto->GetAxis(0)->SetTitle("y-y_{ref} (cm)");
 153   fResolHisto->GetAxis(1)->SetTitle("z-z_{ref} (cm)");
 154   fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{ref} (rad)");
 155   fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{ref} (rad)");
 156   fResolHisto->GetAxis(4)->SetTitle("(p_{T}/p_{Tref}-1)");
 157   fResolHisto->GetAxis(5)->SetTitle("y_{ref} (cm)");
 158   fResolHisto->GetAxis(6)->SetTitle("z_{ref} (cm)");
 159   fResolHisto->GetAxis(7)->SetTitle("#phi_{ref} (rad)");
 160   fResolHisto->GetAxis(8)->SetTitle("#eta_{ref}");
 161   fResolHisto->GetAxis(9)->SetTitle("p_{Tref} (GeV/c)");
 162   fResolHisto->GetAxis(10)->SetTitle("isReconstructed");
 163   fResolHisto->Sumw2();
 164
 165   //pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt:isRec
 166   Int_t binsPullHisto[11]={100,100,100,100,100,50,50,50,50,nPtBins,2};
 167   Double_t minPullHisto[11]={-5.,-5.,-5.,-5.,-5., yMin, zMin,-1.,-2.0, ptMin,0};
 168   Double_t maxPullHisto[11]={ 5., 5., 5., 5., 5., yMax, zMax, 1., 2.0, ptMax,2};
 169   fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt:isRec",11,binsPullHisto,minPullHisto,maxPullHisto);
 170
 171   fPullHisto->GetAxis(0)->SetTitle("(y-y_{ref})/#sigma");
 172   fPullHisto->GetAxis(1)->SetTitle("(z-z_{ref})/#sigma");
 173   fPullHisto->GetAxis(2)->SetTitle("(sin#phi-sin#phi_{ref})/#sigma");
 174   fPullHisto->GetAxis(3)->SetTitle("(tan#lambda-tan#lambda_{ref})/#sigma");
 175   fPullHisto->GetAxis(4)->SetTitle("(p_{Tref}/p_{T}-1)/#sigma");
 176   fPullHisto->GetAxis(5)->SetTitle("y_{ref} (cm)");
 177   fPullHisto->GetAxis(6)->SetTitle("z_{ref} (cm)");
 178   fPullHisto->GetAxis(7)->SetTitle("sin#phi_{ref}");
 179   fPullHisto->GetAxis(8)->SetTitle("tan#lambda_{ref}");
 180   fPullHisto->GetAxis(9)->SetTitle("1/p_{Tref} (GeV/c)^{-1}");
 181   fPullHisto->GetAxis(10)->SetTitle("isReconstructed");
 182   fPullHisto->Sumw2();
 183
 184   // -> has match:y:z:snp:tgl:phi:pt:ITSclusters
 185   Int_t binsTrackingEffHisto[8]    = { 2,    50, 100, 50, 50, 90,           100,  7   };
 186   Double_t minTrackingEffHisto[8]  = {-0.5, -25, -50, -1, -2, 0.,            0,   -0.5 };
 187   Double_t maxTrackingEffHisto[8]  = { 1.5,  25,  50,  1,  2, 2*TMath::Pi(), 20,   6.5 };
 188
 189   fTrackingEffHisto = new THnSparseF("fTrackingEffHisto","has match:y:z:snp:tgl:phi:pt:ITSclusters",8,binsTrackingEffHisto,minTrackingEffHisto,maxTrackingEffHisto);
 190   fTrackingEffHisto->GetAxis(0)->SetTitle("IsMatching");
 191   fTrackingEffHisto->GetAxis(1)->SetTitle("local y (cm)");
 192   fTrackingEffHisto->GetAxis(2)->SetTitle("z (cm)");
 193   fTrackingEffHisto->GetAxis(3)->SetTitle("sin(#phi)");
 194   fTrackingEffHisto->GetAxis(4)->SetTitle("tan(#lambda)");
 195   fTrackingEffHisto->GetAxis(5)->SetTitle("phi (rad)");
 196   fTrackingEffHisto->GetAxis(6)->SetTitle("p_{T}");
 197   fTrackingEffHisto->GetAxis(7)->SetTitle("number of ITS clusters");
 198   fTrackingEffHisto->Sumw2();
 199
 200   // Init cuts
 201   if(!fCutsMC)
 202     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 203   if(!fCutsRC)
 204     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 205
 206   // init folder
 207   fAnalysisFolder = CreateFolder("folderMatch","Analysis Matching Folder");
 208 }
 209
 210 //_____________________________________________________________________________
 211 void AliPerformanceMatch::ProcessITSTPC(Int_t iTrack, AliESDEvent *const esdEvent, AliStack* /*const stack*/, AliESDtrack *const esdTrack, AliESDfriendTrack *const esdFriendTrack)
 212 {
 213   //
 214   // addition to standard analysis - check if ITS stand-alone tracks have a match in the TPC
 215   // Origin: A. Kalwait
 216   // Modified: J. Otwinowski
 217   if(!esdTrack) return;
 218   if(!esdFriendTrack) return;
 219
 220   //
 221   if (esdTrack->GetInnerParam()) return; // ITS stand-alone tracks have not TPC inner param
 222   if (esdTrack->GetITSclusters(0) < fCutsRC->GetMinNClustersITS()) return;
 223   //AliTracker::PropagateTrackToBxByBz(esdTrack,80.0,0.1056,1.0,kTRUE); // we propagate the ITS stand-alone to the inner TPC radius
 224   Bool_t hasMatch = kFALSE;
 225     for (Int_t jTrack = 0; jTrack < esdEvent->GetNumberOfTracks(); jTrack++) {
 226       // loop for all those tracks and check if the corresponding TPC track is found
 227       if (jTrack==iTrack) continue;
 228       AliESDtrack *trackTPC = esdEvent->GetTrack(jTrack);
 229       if (!trackTPC) continue;
 230       if (!trackTPC->GetTPCInnerParam()) continue;
 231
 232       AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(trackTPC->GetTPCInnerParam()));
 233       if(!innerTPC) continue;
 234
 235       AliTracker::PropagateTrackToBxByBz(innerTPC,2.8,trackTPC->GetMass(),1.0,kTRUE);
 236       Double_t x[3]; trackTPC->GetXYZ(x);
 237       Double_t b[3]; AliTracker::GetBxByBz(x,b);
 238
 239       Double_t dz[2],cov[3];
 240       innerTPC->PropagateToDCABxByBz(esdEvent->GetPrimaryVertexSPD(),b,kVeryBig,dz,cov);
 241
 242       // rough checking if they match
 243       /*
 244       printf("+++++++++++++++++++++++++++++++++++++++++++++++ ");
 245       printf("esdTrack->GetY() %f, esdTrack->GetSnp() %f, esdTrack->GetTgl() %f \n",
 246               esdTrack->GetY(), esdTrack->GetSnp(), esdTrack->GetTgl());
 247       printf("innerTPC->GetY() %f, innerTPC->GetSnp() %f, innerTPC->GetTgl() %f \n",
 248               innerTPC->GetY() , innerTPC->GetSnp() , innerTPC->GetTgl());
 249       */
 250       if (TMath::Abs(esdTrack->GetY() - innerTPC->GetY()) > 3) { delete innerTPC; continue; }
 251       if (TMath::Abs(esdTrack->GetSnp() - innerTPC->GetSnp()) > 0.2) { delete innerTPC; continue; }
 252       if (TMath::Abs(esdTrack->GetTgl() - innerTPC->GetTgl()) > 0.2) { delete innerTPC; continue; }
 253
 254       hasMatch = kTRUE;
 255       if(innerTPC) delete innerTPC;
 256     }
 257     //has match:y:z:snp:tgl:phi:pt:ITSclusters
 258     Double_t vecTrackingEff[8] = { hasMatch,esdTrack->GetY(),esdTrack->GetZ(),esdTrack->GetSnp(),esdTrack->GetTgl(),esdTrack->Phi(), esdTrack->Pt(),esdTrack->GetITSclusters(0) };
 259     fTrackingEffHisto->Fill(vecTrackingEff);
 260
 261 }
 262
 263 //_____________________________________________________________________________
 264 void AliPerformanceMatch::ProcessTPCITS(AliStack* /*const stack*/, AliESDtrack *const esdTrack, AliESDfriendTrack *const esdFriendTrack)
 265 {
 266   //
 267   // Match TPC and ITS min-bias tracks
 268   // at radius between detectors
 269   //
 270   if(!esdTrack) return;
 271   if(!esdFriendTrack) return;
 272
 273   //
 274   // Propagate tracks to the radius between TPC-ITS
 275   // using B-field and material budget
 276   //
 277   Double_t radius = fCutsRC->GetTPCITSMatchingRadius();
 278   Double_t mass = esdTrack->GetMass();
 279   Double_t step=1.0; // cm
 280
 281   //
 282   // Propagate TPCinner (reference detector)
 283   //
 284   Bool_t isTPCOK=kFALSE;
 285   AliExternalTrackParam *innerTPC=NULL;
 286
 287   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 288   esdTrack->GetImpactParametersTPC(dca,cov);
 289
 290   //
 291   // select primaries
 292   //
 293   Double_t dcaToVertex = -1;
 294   if( fCutsRC->GetDCAToVertex2D() )
 295   {
 296       dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY()                    + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
 297   }
 298   if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
 299   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
 300   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
 301
 302   if( (esdTrack->GetNcls(1)>fCutsRC->GetMinNClustersTPC()) &&
 303       (esdTrack->GetTPCInnerParam()) &&
 304       (innerTPC=new AliExternalTrackParam(*(esdTrack->GetTPCInnerParam()))))
 305   {
 306      isTPCOK = AliTracker::PropagateTrackToBxByBz(innerTPC,radius,mass,step,kTRUE);
 307   }
 308   if(!isTPCOK) {
 309    if(innerTPC) delete innerTPC;
 310    return;
 311   }
 312
 313   //
 314   // Propagate ITSouter
 315   //
 316   Bool_t isITSOK=kFALSE;
 317   AliExternalTrackParam *outerITS=NULL;
 318
 319   if( (esdTrack->GetNcls(0)>fCutsRC->GetMinNClustersITS()) &&
 320       (esdFriendTrack->GetITSOut()) &&
 321       (outerITS=new AliExternalTrackParam(*(esdFriendTrack->GetITSOut()))))
 322   {
 323     isITSOK = AliTracker::PropagateTrackToBxByBz(outerITS,radius,mass,step,kTRUE);
 324   }
 325
 326   //
 327   // Fill histograms (TPC reference detector)
 328   //
 329   if(isTPCOK)
 330     FillHistograms(innerTPC,outerITS,isITSOK);
 331
 332   if(outerITS) delete outerITS;
 333   if(innerTPC) delete innerTPC;
 334 }
 335
 336 //_____________________________________________________________________________
 337 void AliPerformanceMatch::ProcessTPCTRD(AliStack* /*const stack*/, AliESDtrack *const esdTrack, AliESDfriendTrack *const esdFriendTrack)
 338 {
 339   //
 340   // Match TPC and TRD min-bias tracks
 341   // at radius between detectors. TPC is the reference detector.
 342   //
 343   if(!esdTrack) return;
 344   if(!esdFriendTrack) return;
 345
 346   //
 347   // Propagate tracks to the radius between TPC-TRD
 348   // using B-field and material budget
 349   //
 350   Double_t radius = fCutsRC->GetTPCTRDMatchingRadius();
 351   Double_t mass = esdTrack->GetMass();
 352   Double_t step=1.0; // cm
 353
 354   //
 355   // Propagate TPCouter (reference detector)
 356   //
 357   Bool_t isTPCOK=kFALSE;
 358   AliExternalTrackParam *outerTPC=NULL;
 359
 360   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 361   esdTrack->GetImpactParametersTPC(dca,cov);
 362
 363   //
 364   // select primaries
 365   //
 366   Double_t dcaToVertex = -1;
 367   if( fCutsRC->GetDCAToVertex2D() )
 368   {
 369       dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY()                    + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
 370   }
 371   if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
 372   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
 373   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
 374
 375   if( (esdTrack->GetNcls(1)>fCutsRC->GetMinNClustersTPC()) &&
 376       (esdFriendTrack->GetTPCOut()) &&
 377       (outerTPC=new AliExternalTrackParam(*(esdFriendTrack->GetTPCOut()))))
 378   {
 379      isTPCOK = AliTracker::PropagateTrackToBxByBz(outerTPC,radius,mass,step,kTRUE);
 380   }
 381   if(!isTPCOK)  {
 382     if(outerTPC) delete outerTPC;
 383     return;
 384   }
 385
 386   //
 387   // Propagate TRDinner
 388   //
 389   Bool_t isTRDOK = kFALSE;
 390   AliExternalTrackParam *innerTRD=NULL;
 391
 392   // get TRD track
 393   AliTRDtrackV1 *trdTrack=NULL; //esdFriendTrack = fESDfriend->GetTrack(itrk);
 394   TObject *calObject=NULL;
 395   Int_t icalib = 0;
 396   while((calObject = esdFriendTrack->GetCalibObject(icalib++))) {
 397     if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
 398     if(!(trdTrack = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
 399   }
 400
 401   if( (trdTrack) &&
 402       (trdTrack->GetNumberOfTracklets()>fCutsRC->GetMinNTrackletsTRD()) &&
 403       (trdTrack->GetTracklet(0)) &&
 404       (esdFriendTrack->GetTRDIn()) &&
 405       (innerTRD = new AliExternalTrackParam(*(esdFriendTrack->GetTRDIn()))))
 406   {
 407     isTRDOK = AliTracker::PropagateTrackToBxByBz(innerTRD,radius,mass,step,kTRUE);
 408   }
 409
 410   //
 411   // Fill histograms (TPC reference detector)
 412   //
 413   if(isTPCOK)
 414     FillHistograms(outerTPC,innerTRD,isTRDOK);
 415
 416   if(outerTPC) delete outerTPC;
 417   if(innerTRD) delete innerTRD;
 418 }
 419
 420 //_____________________________________________________________________________
 421 void AliPerformanceMatch::FillHistograms(AliExternalTrackParam *const refParam, AliExternalTrackParam *const param, Bool_t isRec)
 422 {
 423   //
 424   // fill performance histograms
 425   // (TPC always as reference)
 426   //
 427   if(!refParam) return;
 428
 429   //
 430   // Deltas (dy,dz,dphi,dtheta,dpt)
 431   //
 432   Float_t delta[5] = {0};
 433   if(param && isRec) {
 434     delta[0] = param->GetY()-refParam->GetY();
 435     delta[1] = param->GetZ()-refParam->GetZ();
 436     delta[2] = TMath::ATan2(param->Py(),param->Px())-TMath::ATan2(refParam->Py(),refParam->Px());
 437     delta[3] = TMath::ATan2(param->Pz(),param->Pt())-TMath::ATan2(refParam->Pz(),refParam->Pt());
 438     if(refParam->Pt()) delta[4] = (param->Pt()-refParam->Pt())/refParam->Pt();
 439   }
 440   //
 441   // Pulls (y,z,snp,tanl,1/pt)
 442   //
 443   Float_t sigma[5] = {0};
 444   Float_t pull[5] = {0};
 445   if(param && isRec) {
 446     sigma[0] = TMath::Sqrt(param->GetSigmaY2()+refParam->GetSigmaY2());
 447     sigma[1] = TMath::Sqrt(param->GetSigmaZ2()+refParam->GetSigmaZ2());
 448     sigma[2] = TMath::Sqrt(param->GetSigmaSnp2()+refParam->GetSigmaSnp2());
 449     sigma[3] = TMath::Sqrt(param->GetSigmaTgl2()+refParam->GetSigmaTgl2());
 450     sigma[4] = TMath::Sqrt(param->GetSigma1Pt2()+refParam->GetSigma1Pt2());
 451     if(sigma[0]) pull[0] = delta[0] / sigma[0];
 452     if(sigma[1]) pull[1] = delta[1] / sigma[1];
 453     if(sigma[2]) pull[2] = (param->GetSnp()-refParam->GetSnp()) / sigma[2];
 454     if(sigma[3]) pull[3] = (param->GetTgl()-refParam->GetTgl()) / sigma[3];
 455     if(sigma[4]) pull[4] = (param->OneOverPt()-refParam->OneOverPt()) / sigma[4];
 456   }
 457
 458   // Fill histograms
 459   Double_t vResolHisto[11] = {delta[0],delta[1],delta[2],delta[3],delta[4],refParam->GetY(),refParam->GetZ(),refParam->Phi(),refParam->Eta(),refParam->Pt(),isRec};
 460   fResolHisto->Fill(vResolHisto);
 461
 462   Double_t vPullHisto[11] = {pull[0],pull[1],pull[2],pull[3],pull[4],refParam->GetY(),refParam->GetZ(),refParam->GetSnp(),refParam->GetTgl(),refParam->OneOverPt(),isRec};
 463   fPullHisto->Fill(vPullHisto);
 464 }
 465
 466 //_____________________________________________________________________________
 467 void AliPerformanceMatch::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
 468 {
 469   // Process comparison information
 470   //
 471   if(!esdEvent)
 472   {
 473     Error("Exec","esdEvent not available");
 474     return;
 475   }
 476   AliHeader* header = 0;
 477   AliGenEventHeader* genHeader = 0;
 478   AliStack* stack = 0;
 479   TArrayF vtxMC(3);
 480
 481   if(bUseMC)
 482   {
 483     if(!mcEvent) {
 484       Error("Exec","mcEvent not available");
 485       return;
 486     }
 487     // get MC event header
 488     header = mcEvent->Header();
 489     if (!header) {
 490       Error("Exec","Header not available");
 491       return;
 492     }
 493     // MC particle stack
 494     stack = mcEvent->Stack();
 495     if (!stack) {
 496       Error("Exec","Stack not available");
 497       return;
 498     }
 499     // get MC vertex
 500     genHeader = header->GenEventHeader();
 501     if (!genHeader) {
 502       Error("Exec","Could not retrieve genHeader from Header");
 503       return;
 504     }
 505     genHeader->PrimaryVertex(vtxMC);
 506   }
 507
 508   // use ESD friends
 509   if(bUseESDfriend) {
 510     if(!esdFriend) {
 511       Error("Exec","esdFriend not available");
 512       return;
 513     }
 514   }
 515
 516   // trigger
 517   if(!bUseMC && GetTriggerClass()) {
 518     Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
 519     if(!isEventTriggered) return;
 520   }
 521
 522   // get TPC event vertex
 523   const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTPC();
 524   if(vtxESD && (vtxESD->GetStatus()<=0)) return;
 525
 526   //  Process events
 527   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 528   {
 529     AliESDtrack *track = esdEvent->GetTrack(iTrack);
 530     if(!track) continue;
 531
 532     AliESDfriendTrack *friendTrack=0;
 533     if(bUseESDfriend) {
 534       friendTrack=esdFriend->GetTrack(iTrack);
 535       if(!friendTrack) continue;
 536     }
 537
 538     if(GetAnalysisMode() == 0) ProcessTPCITS(stack,track,friendTrack);
 539     else if(GetAnalysisMode() == 1) ProcessTPCTRD(stack,track,friendTrack);
 540     else if(GetAnalysisMode() == 2) ProcessITSTPC(iTrack,esdEvent,stack,track,friendTrack);
 541     else {
 542       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 543       return;
 544     }
 545   }
 546 }
 547
 548 //_____________________________________________________________________________
 549 TH1F* AliPerformanceMatch::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
 550   // Create resolution histograms
 551
 552   TH1F *hisr, *hism;
 553   if (!gPad) new TCanvas;
 554   hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
 555   if (type) return hism;
 556   else
 557     return hisr;
 558 }
 559
 560 //_____________________________________________________________________________
 561 void AliPerformanceMatch::Analyse() {
 562   // Analyse comparison information and store output histograms
 563   // in the folder "folderMatch"
 564   //
 565   TH1::AddDirectory(kFALSE);
 566   TH1F *h=0;
 567   TH1F *h2=0;
 568   TH2F *h2D=0;
 569   TObjArray *aFolderObj = new TObjArray;
 570
 571   // write results in the folder
 572   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
 573   c->cd();
 574
 575   char name[256];
 576   char title[256];
 577
 578   if(GetAnalysisMode()==0 || GetAnalysisMode()==1) {
 579
 580   fResolHisto->GetAxis(10)->SetRangeUser(1.0,2.0); // only reconstructed
 581   fPullHisto->GetAxis(10)->SetRangeUser(1.0,2.0);  // only reconstructed
 582   for(Int_t i=0; i<5; i++)
 583   {
 584     for(Int_t j=5; j<10; j++)
 585     {
 586       //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
 587       if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.0,0.89); // eta window
 588       else fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.49);
 589       fResolHisto->GetAxis(9)->SetRangeUser(0.1,100.); // pt threshold
 590
 591       h2D = (TH2F*)fResolHisto->Projection(i,j);
 592
 593       h = AliPerformanceMatch::MakeResol(h2D,1,0,100);
 594       sprintf(name,"h_res_%d_vs_%d",i,j);
 595       h->SetName(name);
 596
 597       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 598       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
 599       h->GetYaxis()->SetTitle(title);
 600       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 601       h->SetTitle(title);
 602
 603       //if(j==9) h->SetBit(TH1::kLogX);
 604       aFolderObj->Add(h);
 605
 606       h = AliPerformanceMatch::MakeResol(h2D,1,1,100);
 607       //h = (TH1F*)arr->At(1);
 608       sprintf(name,"h_mean_res_%d_vs_%d",i,j);
 609       h->SetName(name);
 610
 611       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 612       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
 613       h->GetYaxis()->SetTitle(title);
 614
 615       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 616       h->SetTitle(title);
 617
 618       if(j==9) h->SetBit(TH1::kLogX);
 619       aFolderObj->Add(h);
 620
 621       //
 622       if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(0.0,0.89); // eta window
 623       else  fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.49); // eta window
 624       fPullHisto->GetAxis(9)->SetRangeUser(0.1,100.);  // pt threshold
 625
 626       h2D = (TH2F*)fPullHisto->Projection(i,j);
 627
 628       h = AliPerformanceMatch::MakeResol(h2D,1,0,100);
 629       sprintf(name,"h_pull_%d_vs_%d",i,j);
 630       h->SetName(name);
 631
 632       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 633       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
 634       h->GetYaxis()->SetTitle(title);
 635       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 636       h->SetTitle(title);
 637
 638       if(j==9) h->SetBit(TH1::kLogX);
 639       aFolderObj->Add(h);
 640
 641       h = AliPerformanceMatch::MakeResol(h2D,1,1,100);
 642       sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
 643       h->SetName(name);
 644
 645       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 646       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
 647       h->GetYaxis()->SetTitle(title);
 648       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 649       h->SetTitle(title);
 650
 651       //if(j==9) h->SetBit(TH1::kLogX);
 652       aFolderObj->Add(h);
 653     }
 654   }
 655
 656   //
 657   // Efficiency plots
 658   //
 659   for(Int_t i=5; i<10; i++)
 660   {
 661     if(i!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
 662     else fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.49);
 663     fResolHisto->GetAxis(9)->SetRangeUser(0.1,100.); // pt threshold
 664
 665     fResolHisto->GetAxis(10)->SetRange(1,fResolHisto->GetAxis(10)->GetNbins()); // all
 666     h = (TH1F*)fResolHisto->Projection(i);
 667
 668     fResolHisto->GetAxis(10)->SetRange(2,2); // only reconstructed
 669     h2 = (TH1F*)fResolHisto->Projection(i);
 670
 671     TH1F* h2c = (TH1F*)h2->Clone();
 672     h2c->Divide(h2,h,1,1,"B");
 673
 674     sprintf(name,"h_eff_%d",i);
 675     h2c->SetName(name);
 676
 677     h2c->GetXaxis()->SetTitle(h2c->GetXaxis()->GetTitle());
 678     h2c->GetYaxis()->SetTitle("efficiency");
 679     h2c->SetTitle("matching effciency");
 680
 681     aFolderObj->Add(h2c);
 682   }
 683
 684   }
 685
 686   //
 687   // TPC efficiency wrt ITS
 688   //
 689   if(GetAnalysisMode()==2) {
 690
 691     h = (TH1F*)fTrackingEffHisto->Projection(0);
 692     aFolderObj->Add(h);
 693
 694     for(Int_t i=1; i<7; i++)
 695     {
 696       //
 697       //
 698       // calculate efficiency
 699       //
 700
 701       // all ITS standalone tracks
 702       fTrackingEffHisto->GetAxis(0)->SetRange(1,fTrackingEffHisto->GetAxis(0)->GetNbins());
 703       h = (TH1F*)fTrackingEffHisto->Projection(i);
 704
 705       // TPC tracks which has matching with TPC
 706       fTrackingEffHisto->GetAxis(0)->SetRange(2,2);
 707       h2 = (TH1F*)fTrackingEffHisto->Projection(i);
 708
 709       TH1F* h2c = (TH1F*)h2->Clone();
 710       h2c->Divide(h2,h,1,1,"B");
 711
 712       sprintf(name,"h_TPC_eff_%d",i);
 713       h2c->SetName(name);
 714
 715       h2c->GetXaxis()->SetTitle(h2c->GetXaxis()->GetTitle());
 716       h2c->GetYaxis()->SetTitle("efficiency");
 717       h2c->SetTitle("TPC effciency wrt ITS");
 718
 719       aFolderObj->Add(h2c);
 720     }
 721
 722   }
 723
 724   // export objects to analysis folder
 725   fAnalysisFolder = ExportToFolder(aFolderObj);
 726
 727   // delete only TObjArray
 728   if(aFolderObj) delete aFolderObj;
 729 }
 730
 731 //_____________________________________________________________________________
 732 TFolder* AliPerformanceMatch::ExportToFolder(TObjArray * array)
 733 {
 734   // recreate folder avery time and export objects to new one
 735   //
 736   AliPerformanceMatch * comp=this;
 737   TFolder *folder = comp->GetAnalysisFolder();
 738
 739   TString name, title;
 740   TFolder *newFolder = 0;
 741   Int_t i = 0;
 742   Int_t size = array->GetSize();
 743
 744   if(folder) {
 745      // get name and title from old folder
 746      name = folder->GetName();
 747      title = folder->GetTitle();
 748
 749          // delete old one
 750      delete folder;
 751
 752          // create new one
 753      newFolder = CreateFolder(name.Data(),title.Data());
 754      newFolder->SetOwner();
 755
 756          // add objects to folder
 757      while(i < size) {
 758            newFolder->Add(array->At(i));
 759            i++;
 760          }
 761   }
 762
 763 return newFolder;
 764 }
 765
 766 //_____________________________________________________________________________
 767 Long64_t AliPerformanceMatch::Merge(TCollection* const list)
 768 {
 769   // Merge list of objects (needed by PROOF)
 770
 771   if (!list)
 772   return 0;
 773
 774   if (list->IsEmpty())
 775   return 1;
 776
 777   TIterator* iter = list->MakeIterator();
 778   TObject* obj = 0;
 779
 780   // collection of generated histograms
 781   Int_t count=0;
 782   while((obj = iter->Next()) != 0)
 783   {
 784   AliPerformanceMatch* entry = dynamic_cast<AliPerformanceMatch*>(obj);
 785   if (entry == 0) continue;
 786
 787   fResolHisto->Add(entry->fResolHisto);
 788   fPullHisto->Add(entry->fPullHisto);
 789   fTrackingEffHisto->Add(entry->fTrackingEffHisto);
 790
 791   count++;
 792   }
 793
 794 return count;
 795 }
 796
 797 //_____________________________________________________________________________
 798 TFolder* AliPerformanceMatch::CreateFolder(TString name,TString title) {
 799 // create folder for analysed histograms
 800 //
 801 TFolder *folder = 0;
 802   folder = new TFolder(name.Data(),title.Data());
 803
 804   return folder;
 805 }