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