PWG1/AliPerformanceMC.cxx

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