PWG1/TPC/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,90,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   // Process pure MC information
 190   //
 191   AliHeader* header = 0;
 192   AliGenEventHeader* genHeader = 0;
 193   AliStack* stack = 0;
 194   TArrayF vtxMC(3);
 195
 196   if(bUseMC)
 197   {
 198     if(!mcEvent) {
 199       Error("Exec","mcEvent not available");
 200       return;
 201     }
 202     // get MC event header
 203     header = mcEvent->Header();
 204     if (!header) {
 205       Error("Exec","Header not available");
 206       return;
 207     }
 208     // MC particle stack
 209     stack = mcEvent->Stack();
 210     if (!stack) {
 211       Error("Exec","Stack not available");
 212       return;
 213     }
 214     // get MC vertex
 215     genHeader = header->GenEventHeader();
 216     if (!genHeader) {
 217       Error("Exec","Could not retrieve genHeader from Header");
 218       return;
 219     }
 220     genHeader->PrimaryVertex(vtxMC);
 221   }
 222   else {
 223     Error("Exec","MC information required!");
 224     return;
 225   }
 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 = (AliMCParticle*) 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::PropagateTrackToBxByBz()
 326       //
 327       if(refTPCIn && refTPCOut)
 328         ProcessInnerTPC(refTPCIn,refTPCOut,particle);
 329     }
 330
 331     if(GetAnalysisMode() == 4) {
 332       // propagate refTPCIn to refTPCOut using AliExternalTrackParam::PropagateToBxByBz()
 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   AliTracker::PropagateTrackToBxByBz(track, radius+step, mass, step, kTRUE,0.99);
 360   AliTracker::PropagateTrackToBxByBz(track, radius+0.5, mass, step*0.1, kTRUE,0.99);
 361   Double_t xyz[3] = {particle->Vx(),particle->Vy(),particle->Vz()};
 362   Double_t sxyz[3]={0.0,0.0,0.0};
 363   AliESDVertex vertex(xyz,sxyz);
 364   Double_t dca[2], cov[3];
 365   //Bool_t isOK = track->PropagateToDCA(&vertex,AliTracker::GetBz(),10,dca,cov);
 366   Double_t field[3];  track->GetBxByBz(field);
 367   Bool_t isOK = track->PropagateToDCABxByBz(&vertex,field,10,dca,cov);
 368   if(!isOK) return;
 369
 370   // Fill histogram
 371   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 372   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 373
 374   Float_t mceta =  particle->Eta();
 375   Float_t mcphi =  particle->Phi();
 376   if(mcphi<0) mcphi += 2.*TMath::Pi();
 377   Float_t mcpt = particle->Pt();
 378   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
 379   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
 380
 381   //if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 382   //{
 383     if(mcpt == 0) return;
 384
 385     //deltaYTPC= track->GetY()-particle->Vy();
 386     deltaYTPC= track->GetY(); // it is already rotated
 387     deltaZTPC = track->GetZ()-particle->Vz();
 388     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 389     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
 390     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
 391
 392     /*
 393     printf("------------------------------ \n");
 394     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() );
 395     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 );
 396     */
 397
 398
 399     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
 400     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 401
 402     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
 403     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
 404
 405     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 406     else pull1PtTPC = 0.;
 407
 408     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
 409     fResolHisto->Fill(vResolHisto);
 410
 411     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
 412     fPullHisto->Fill(vPullHisto);
 413   //}
 414   if(track) delete track;
 415 }
 416
 417 //_____________________________________________________________________________
 418 void AliPerformanceMC::ProcessInnerTPC(AliTrackReference* const refIn,  AliTrackReference *const  refOut, TParticle *const particle) {
 419   //
 420   // Test propagation from Out to In
 421   //
 422   if(!refIn) return;
 423   if(!refOut) return;
 424   if(!particle) return;
 425
 426   AliExternalTrackParam *track=MakeTrack(refOut,particle);
 427   if (!track) return;
 428
 429   Double_t mass = particle->GetMass();
 430   Double_t step=1;
 431
 432   Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
 433   //Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
 434   //track->Rotate(alphaIn-alphaOut);
 435   track->Rotate(alphaIn);
 436   //printf("alphaIn %f, alphaOut %f \n",alphaIn,alphaOut);
 437
 438   //
 439   //Bool_t isOK = AliTracker::PropagateTrackTo(track, refIn->R(), mass, step, kTRUE,0.99);
 440   //Bool_t isOK = AliTracker::PropagateTrackTo(track, refIn->R(), mass, step, kFALSE,0.99);
 441   Bool_t isOK = AliTracker::PropagateTrackToBxByBz(track, refIn->R(), mass, step, kFALSE,0.99);
 442   if(!isOK) return;
 443
 444   // calculate alpha angle
 445   //Double_t xyz[3] = {refIn->X(),refIn->Y(),refIn0->Z()};
 446   //Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 447   //if (alpha<0) alpha += TMath::TwoPi();
 448
 449   // rotate outer track to inner
 450   // and propagate to the radius of the first track referenco of TPC
 451   //Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
 452   //Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
 453   //if(!isOK) return;
 454
 455   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * refIn->Theta()));
 456   Float_t mcphi =  refIn->Phi();
 457   if(mcphi<0) mcphi += 2.*TMath::Pi();
 458   Float_t mcpt = refIn->Pt();
 459   Float_t mcsnp = TMath::Sin(TMath::ATan2(refIn->Py(),refIn->Px()));
 460   Float_t mctgl = TMath::Tan(TMath::ATan2(refIn->Pz(),refIn->Pt()));
 461
 462   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 463   Float_t pull1PtTPC=0., pullYTPC=0., pullZTPC=0., pullPhiTPC=0., pullLambdaTPC=0.;
 464
 465   if(mcpt == 0) return;
 466
 467   //deltaYTPC = track->GetY()-refIn->Y();
 468   deltaYTPC = track->GetY();
 469   deltaZTPC = track->GetZ()-refIn->Z();
 470   deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(refIn->Pz(),refIn->Pt());
 471   deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(refIn->Py(),refIn->Px());
 472   deltaPtTPC = (track->Pt()-mcpt) / mcpt;
 473
 474   /*
 475   printf("------------------------------ \n");
 476   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() );
 477   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 );
 478   */
 479
 480   if(TMath::Sqrt(track->GetSigmaY2())) pullYTPC = track->GetY() / TMath::Sqrt(track->GetSigmaY2());
 481   if(TMath::Sqrt(track->GetSigmaZ2())) pullZTPC = (track->GetZ()-refIn->Z()) / TMath::Sqrt(track->GetSigmaZ2());
 482   if(TMath::Sqrt(track->GetSigmaSnp2())) pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
 483   if(TMath::Sqrt(track->GetSigmaTgl2())) pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
 484   if(TMath::Sqrt(track->GetSigma1Pt2())) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 485
 486   //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);
 487
 488   Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,refIn->Y(),refIn->Z(),mcphi,mceta,mcpt};
 489   fResolHisto->Fill(vResolHisto);
 490
 491   Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt};
 492   fPullHisto->Fill(vPullHisto);
 493
 494   if(track) delete track;
 495 }
 496
 497 //_____________________________________________________________________________
 498 void AliPerformanceMC::ProcessOuterTPCExt(TParticle *const part, TClonesArray * const trefs)
 499 {
 500   const Int_t kMinRefs=5;
 501   Int_t nrefs = trefs->GetEntries();
 502   if (nrefs<kMinRefs) return; // we should have enough references
 503   Int_t iref0 =-1;
 504   Int_t iref1 =-1;
 505
 506   for (Int_t iref=0; iref<nrefs; iref++){
 507     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 508     if (!ref) continue;
 509     Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
 510     if (dir<0) break;
 511     if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
 512     if (iref0<0) iref0 = iref;
 513     iref1 = iref;
 514   }
 515   if (iref1-iref0<kMinRefs) return;
 516   Double_t covar[15];
 517   for (Int_t icov=0; icov<15; icov++) covar[icov]=0;
 518   covar[0]=1;
 519   covar[2]=1;
 520   covar[5]=1;
 521   covar[9]=1;
 522   covar[14]=1;
 523
 524   AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
 525   AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
 526   AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
 527   AliExternalTrackParam *paramUpdate   = MakeTrack(refIn,part);
 528   paramUpdate->AddCovariance(covar);
 529
 530   Bool_t isOKP=kTRUE;
 531   Bool_t isOKU=kTRUE;
 532   AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
 533   if(!field) return;
 534   for (Int_t iref = iref0; iref<=iref1; iref++){
 535     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 536     if(!ref) continue;
 537     Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
 538     Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
 539     Double_t mag[3];
 540     field->Field(pos,mag);
 541     isOKP&=paramPropagate->Rotate(alphaC);
 542     isOKU&=paramUpdate->Rotate(alphaC);
 543     for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
 544       //isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
 545       //isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
 546       isOKP&=paramPropagate->PropagateToBxByBz(xref, mag);
 547       isOKU&=paramUpdate->PropagateToBxByBz(xref, mag);
 548     }
 549     //isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
 550     //isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
 551     isOKP&=paramPropagate->PropagateToBxByBz(ref->R(), mag);
 552     isOKU&=paramUpdate->PropagateToBxByBz(ref->R(), mag);
 553     Double_t clpos[2] = {0, ref->Z()};
 554     Double_t clcov[3] = { 0.005,0,0.005};
 555     isOKU&= paramUpdate->Update(clpos, clcov);
 556   }
 557
 558   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * refOut->Theta()));
 559   Float_t mcphi =  refOut->Phi();
 560   if(mcphi<0) mcphi += 2.*TMath::Pi();
 561   Float_t mcpt = refOut->Pt();
 562   Float_t mcsnp = TMath::Sin(TMath::ATan2(refOut->Py(),refOut->Px()));
 563   Float_t mctgl = TMath::Tan(TMath::ATan2(refOut->Pz(),refOut->Pt()));
 564
 565   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 566   Float_t pull1PtTPC=0., pullYTPC=0., pullZTPC=0., pullPhiTPC=0., pullLambdaTPC=0.;
 567
 568   if(mcpt == 0) return;
 569
 570   //deltaYTPC = track->GetY()-refIn->Y();
 571   deltaYTPC = paramPropagate->GetY();
 572   deltaZTPC = paramPropagate->GetZ()-refOut->Z();
 573   deltaLambdaTPC = TMath::ATan2(paramPropagate->Pz(),paramPropagate->Pt())-TMath::ATan2(refOut->Pz(),refOut->Pt());
 574   deltaPhiTPC = TMath::ATan2(paramPropagate->Py(),paramPropagate->Px())-TMath::ATan2(refOut->Py(),refOut->Px());
 575   deltaPtTPC = (paramPropagate->Pt()-mcpt) / mcpt;
 576
 577   if(TMath::Sqrt(paramPropagate->GetSigmaY2())) pullYTPC = paramPropagate->GetY() / TMath::Sqrt(paramPropagate->GetSigmaY2());
 578   if(TMath::Sqrt(paramPropagate->GetSigmaZ2())) pullZTPC = (paramPropagate->GetZ()-refOut->Z()) / TMath::Sqrt(paramPropagate->GetSigmaZ2());
 579   if(TMath::Sqrt(paramPropagate->GetSigmaSnp2())) pullPhiTPC = (paramPropagate->GetSnp() - mcsnp) / TMath::Sqrt(paramPropagate->GetSigmaSnp2());
 580   if(TMath::Sqrt(paramPropagate->GetSigmaTgl2())) pullLambdaTPC = (paramPropagate->GetTgl() - mctgl) / TMath::Sqrt(paramPropagate->GetSigmaTgl2());
 581   if(TMath::Sqrt(paramPropagate->GetSigma1Pt2())) pull1PtTPC = (paramPropagate->OneOverPt()-1./mcpt) / TMath::Sqrt(paramPropagate->GetSigma1Pt2());
 582
 583   //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);
 584
 585   Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,refIn->Y(),refIn->Z(),mcphi,mceta,mcpt};
 586   fResolHisto->Fill(vResolHisto);
 587
 588   Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt};
 589   fPullHisto->Fill(vPullHisto);
 590 }
 591
 592 //_____________________________________________________________________________
 593 AliExternalTrackParam * AliPerformanceMC::MakeTrack(const AliTrackReference* const ref, TParticle* const part)
 594 {
 595   //
 596   // Make track out of the track ref
 597   // part - TParticle used to determine chargr
 598   // the covariance matrix - equal 0 - starting from ideal MC position
 599   Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
 600   Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
 601   Double_t cv[21];
 602   for (Int_t i=0; i<21;i++) cv[i]=0;
 603   if (!part->GetPDG()) return 0;
 604   AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
 605   return param;
 606 }
 607
 608 //_____________________________________________________________________________
 609 TH1F* AliPerformanceMC::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
 610   // Create resolution histograms
 611
 612   TH1F *hisr, *hism;
 613   if (!gPad) new TCanvas;
 614   hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
 615   if (type) return hism;
 616   else
 617     return hisr;
 618 }
 619
 620 //_____________________________________________________________________________
 621 void AliPerformanceMC::Analyse() {
 622   // Analyse comparison information and store output histograms
 623   // in the folder "folderRes"
 624   //
 625   TH1::AddDirectory(kFALSE);
 626   TH1F *h=0;
 627   TH2F *h2D=0;
 628   TObjArray *aFolderObj = new TObjArray;
 629   if(!aFolderObj) return;
 630
 631   // write results in the folder
 632   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
 633   c->cd();
 634
 635   char name[256];
 636   char title[256];
 637
 638   for(Int_t i=0; i<5; i++)
 639   {
 640     for(Int_t j=5; j<10; j++)
 641     {
 642       //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
 643       fResolHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
 644       if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.,0.9); // eta window
 645       //fResolHisto->GetAxis(9)->SetRangeUser(0.16,3.); // pt window
 646       if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
 647
 648       h2D = (TH2F*)fResolHisto->Projection(i,j);
 649
 650       h = AliPerformanceMC::MakeResol(h2D,1,0,20);
 651       snprintf(name,256,"h_res_%d_vs_%d",i,j);
 652       h->SetName(name);
 653
 654       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 655       snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
 656       h->GetYaxis()->SetTitle(title);
 657       snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 658       h->SetTitle(title);
 659
 660       if(j==9) h->SetBit(TH1::kLogX);
 661       aFolderObj->Add(h);
 662
 663       h = AliPerformanceMC::MakeResol(h2D,1,1,20);
 664       //h = (TH1F*)arr->At(1);
 665       snprintf(name,256,"h_mean_res_%d_vs_%d",i,j);
 666       h->SetName(name);
 667
 668       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 669       snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
 670       h->GetYaxis()->SetTitle(title);
 671
 672       snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 673       h->SetTitle(title);
 674
 675       if(j==9) h->SetBit(TH1::kLogX);
 676       aFolderObj->Add(h);
 677
 678       fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
 679       fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
 680
 681       //
 682       if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(-1.0,0.99); // theta window
 683       else  fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.5); // theta window
 684       fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.);  // pt threshold
 685       if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
 686
 687       h2D = (TH2F*)fPullHisto->Projection(i,j);
 688
 689       h = AliPerformanceMC::MakeResol(h2D,1,0,20);
 690       snprintf(name,256,"h_pull_%d_vs_%d",i,j);
 691       h->SetName(name);
 692
 693       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 694       snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
 695       h->GetYaxis()->SetTitle(title);
 696       snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 697       h->SetTitle(title);
 698
 699       //if(j==9) h->SetBit(TH1::kLogX);
 700       aFolderObj->Add(h);
 701
 702       h = AliPerformanceMC::MakeResol(h2D,1,1,20);
 703       snprintf(name,256,"h_mean_pull_%d_vs_%d",i,j);
 704       h->SetName(name);
 705
 706       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 707       snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
 708       h->GetYaxis()->SetTitle(title);
 709       snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 710       h->SetTitle(title);
 711
 712       //if(j==9) h->SetBit(TH1::kLogX);
 713       aFolderObj->Add(h);
 714     }
 715   }
 716
 717   // export objects to analysis folder
 718   fAnalysisFolder = ExportToFolder(aFolderObj);
 719
 720   // delete only TObjArray
 721   if(aFolderObj) delete aFolderObj;
 722 }
 723
 724 //_____________________________________________________________________________
 725 TFolder* AliPerformanceMC::ExportToFolder(TObjArray * array)
 726 {
 727   // recreate folder avery time and export objects to new one
 728   //
 729   AliPerformanceMC * comp=this;
 730   TFolder *folder = comp->GetAnalysisFolder();
 731
 732   TString name, title;
 733   TFolder *newFolder = 0;
 734   Int_t i = 0;
 735   Int_t size = array->GetSize();
 736
 737   if(folder) {
 738      // get name and title from old folder
 739      name = folder->GetName();
 740      title = folder->GetTitle();
 741
 742          // delete old one
 743      delete folder;
 744
 745          // create new one
 746      newFolder = CreateFolder(name.Data(),title.Data());
 747      newFolder->SetOwner();
 748
 749          // add objects to folder
 750      while(i < size) {
 751            newFolder->Add(array->At(i));
 752            i++;
 753          }
 754   }
 755
 756 return newFolder;
 757 }
 758
 759 //_____________________________________________________________________________
 760 Long64_t AliPerformanceMC::Merge(TCollection* const list)
 761 {
 762   // Merge list of objects (needed by PROOF)
 763
 764   if (!list)
 765   return 0;
 766
 767   if (list->IsEmpty())
 768   return 1;
 769
 770   TIterator* iter = list->MakeIterator();
 771   TObject* obj = 0;
 772
 773   // collection of generated histograms
 774   Int_t count=0;
 775   while((obj = iter->Next()) != 0)
 776   {
 777   AliPerformanceMC* entry = dynamic_cast<AliPerformanceMC*>(obj);
 778   if (entry == 0) continue;
 779
 780   fResolHisto->Add(entry->fResolHisto);
 781   fPullHisto->Add(entry->fPullHisto);
 782
 783   count++;
 784   }
 785
 786 return count;
 787 }
 788
 789 //_____________________________________________________________________________
 790 TFolder* AliPerformanceMC::CreateFolder(TString name,TString title) {
 791 // create folder for analysed histograms
 792 //
 793 TFolder *folder = 0;
 794   folder = new TFolder(name.Data(),title.Data());
 795
 796   return folder;
 797 }