PWG0/dNdEta/AlidNdEtaTask.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaTask.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TVector3.h>
   9 #include <TChain.h>
  10 #include <TFile.h>
  11 #include <TH1F.h>
  12 #include <TH2F.h>
  13 #include <TH3F.h>
  14 #include <TParticle.h>
  15 #include <TRandom.h>
  16 #include <TNtuple.h>
  17 #include <TObjString.h>
  18 #include <TF1.h>
  19
  20 #include <AliLog.h>
  21 #include <AliESDVertex.h>
  22 #include <AliESDEvent.h>
  23 #include <AliStack.h>
  24 #include <AliHeader.h>
  25 #include <AliGenEventHeader.h>
  26 #include <AliMultiplicity.h>
  27 #include <AliAnalysisManager.h>
  28 #include <AliMCEventHandler.h>
  29 #include <AliMCEvent.h>
  30 #include <AliESDInputHandler.h>
  31
  32 #include "esdTrackCuts/AliESDtrackCuts.h"
  33 #include "AliPWG0Helper.h"
  34 #include "AliCorrection.h"
  35 #include "AliCorrectionMatrix3D.h"
  36 #include "dNdEta/dNdEtaAnalysis.h"
  37
  38 ClassImp(AlidNdEtaTask)
  39
  40 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
  41   AliAnalysisTask("AlidNdEtaTask", ""),
  42   fESD(0),
  43   fOutput(0),
  44   fOption(opt),
  45   fAnalysisMode(AliPWG0Helper::kTPC),
  46   fReadMC(kFALSE),
  47   fUseMCVertex(kFALSE),
  48   fEsdTrackCuts(0),
  49   fdNdEtaAnalysisESD(0),
  50   fMult(0),
  51   fMultVtx(0),
  52   fEvents(0),
  53   fVertexResolution(0),
  54   fdNdEtaAnalysis(0),
  55   fdNdEtaAnalysisTr(0),
  56   fdNdEtaAnalysisTrVtx(0),
  57   fdNdEtaAnalysisTracks(0),
  58   fVertex(0),
  59   fPartPt(0),
  60   fDeltaPhi(0)
  61 {
  62   //
  63   // Constructor. Initialization of pointers
  64   //
  65
  66   // Define input and output slots here
  67   DefineInput(0, TChain::Class());
  68   DefineOutput(0, TList::Class());
  69 }
  70
  71 AlidNdEtaTask::~AlidNdEtaTask()
  72 {
  73   //
  74   // Destructor
  75   //
  76
  77   // histograms are in the output list and deleted when the output
  78   // list is deleted by the TSelector dtor
  79
  80   if (fOutput) {
  81     delete fOutput;
  82     fOutput = 0;
  83   }
  84 }
  85
  86 //________________________________________________________________________
  87 void AlidNdEtaTask::ConnectInputData(Option_t *)
  88 {
  89   // Connect ESD
  90   // Called once
  91
  92   Printf("AlidNdEtaTask::ConnectInputData called");
  93
  94   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  95   if (!tree) {
  96     Printf("ERROR: Could not read tree from input slot 0");
  97   } else {
  98     // Disable all branches and enable only the needed ones
  99     //tree->SetBranchStatus("*", 0);
 100
 101     tree->SetBranchStatus("fTriggerMask", 1);
 102     tree->SetBranchStatus("fSPDVertex*", 1);
 103     // PrimaryVertex also needed
 104
 105     if (fAnalysisMode == AliPWG0Helper::kSPD) {
 106       tree->SetBranchStatus("fSPDMult*", 1);
 107       //AliPWG0Helper::SetBranchStatusRecursive(tree, "AliMultiplicity", kTRUE);
 108     }
 109
 110     if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
 111       AliESDtrackCuts::EnableNeededBranches(tree);
 112       tree->SetBranchStatus("fTracks.fLabel", 1);
 113       //AliPWG0Helper::SetBranchStatusRecursive(tree, "Tracks", kTRUE);
 114     }
 115
 116     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 117
 118     if (!esdH) {
 119       Printf("ERROR: Could not get ESDInputHandler");
 120     } else
 121       fESD = esdH->GetEvent();
 122   }
 123
 124   // disable info messages of AliMCEvent (per event)
 125   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 126 }
 127
 128 void AlidNdEtaTask::CreateOutputObjects()
 129 {
 130   // create result objects and add to output list
 131
 132   Printf("AlidNdEtaTask::CreateOutputObjects");
 133
 134   fOutput = new TList;
 135   fOutput->SetOwner();
 136
 137   fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
 138   fOutput->Add(fdNdEtaAnalysisESD);
 139
 140   fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
 141   fOutput->Add(fMult);
 142
 143   fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
 144   fOutput->Add(fMultVtx);
 145
 146   for (Int_t i=0; i<3; ++i)
 147   {
 148     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
 149     fPartEta[i]->Sumw2();
 150     fOutput->Add(fPartEta[i]);
 151   }
 152
 153   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 40, -20, 20);
 154   fOutput->Add(fEvents);
 155
 156   fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
 157   fOutput->Add(fVertexResolution);
 158
 159   if (fReadMC)
 160   {
 161     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
 162     fOutput->Add(fdNdEtaAnalysis);
 163
 164     fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
 165     fOutput->Add(fdNdEtaAnalysisTr);
 166
 167     fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
 168     fOutput->Add(fdNdEtaAnalysisTrVtx);
 169
 170     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
 171     fOutput->Add(fdNdEtaAnalysisTracks);
 172
 173     fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
 174     fOutput->Add(fVertex);
 175
 176     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
 177     fPartPt->Sumw2();
 178     fOutput->Add(fPartPt);
 179   }
 180
 181   if (fAnalysisMode == AliPWG0Helper::kSPD)
 182     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 1000, -3.14, 3.14);
 183 }
 184
 185 void AlidNdEtaTask::Exec(Option_t*)
 186 {
 187   // process the event
 188
 189   // Check prerequisites
 190   if (!fESD)
 191   {
 192     AliDebug(AliLog::kError, "ESD branch not available");
 193     return;
 194   }
 195
 196   // trigger definition
 197   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2);
 198
 199   // get the ESD vertex
 200   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 201
 202   Double_t vtx[3];
 203   if (vtxESD) {
 204     vtxESD->GetXYZ(vtx);
 205   }
 206   else
 207     Printf("No vertex");
 208
 209   // fill z vertex resolution
 210   if (vtxESD)
 211     fVertexResolution->Fill(vtxESD->GetZRes());
 212
 213   // post the data already here
 214   PostData(0, fOutput);
 215
 216   // create list of (label, eta, pt) tuples
 217   Int_t inputCount = 0;
 218   Int_t* labelArr = 0;
 219   Float_t* etaArr = 0;
 220   Float_t* ptArr = 0;
 221   if (fAnalysisMode == AliPWG0Helper::kSPD)
 222   {
 223     // get tracklets
 224     const AliMultiplicity* mult = fESD->GetMultiplicity();
 225     if (!mult)
 226     {
 227       AliDebug(AliLog::kError, "AliMultiplicity not available");
 228       return;
 229     }
 230
 231     labelArr = new Int_t[mult->GetNumberOfTracklets()];
 232     etaArr = new Float_t[mult->GetNumberOfTracklets()];
 233     ptArr = new Float_t[mult->GetNumberOfTracklets()];
 234
 235     // get multiplicity from ITS tracklets
 236     for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 237     {
 238       //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 239
 240       // This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
 241       if (mult->GetDeltaPhi(i) < -1000)
 242         continue;
 243
 244       fDeltaPhi->Fill(mult->GetDeltaPhi(i));
 245
 246       etaArr[inputCount] = mult->GetEta(i);
 247       labelArr[inputCount] = mult->GetLabel(i);
 248       ptArr[inputCount] = 0; // no pt for tracklets
 249       ++inputCount;
 250     }
 251
 252     Printf("Accepted %d tracklets", inputCount);
 253   }
 254   else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
 255   {
 256     if (!fEsdTrackCuts)
 257     {
 258       AliDebug(AliLog::kError, "fESDTrackCuts not available");
 259       return;
 260     }
 261
 262     // get multiplicity from ESD tracks
 263     TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
 264     Int_t nGoodTracks = list->GetEntries();
 265
 266     labelArr = new Int_t[nGoodTracks];
 267     etaArr = new Float_t[nGoodTracks];
 268     ptArr = new Float_t[nGoodTracks];
 269
 270     // loop over esd tracks
 271     for (Int_t i=0; i<nGoodTracks; i++)
 272     {
 273       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 274       if (!esdTrack)
 275       {
 276         AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 277         continue;
 278       }
 279
 280       etaArr[inputCount] = esdTrack->Eta();
 281       labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 282       ptArr[inputCount] = esdTrack->Pt();
 283       ++inputCount;
 284     }
 285
 286     Printf("Accepted %d tracks", nGoodTracks);
 287
 288     delete list;
 289   }
 290   else
 291     return;
 292
 293   if (fUseMCVertex) {
 294     Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
 295     if (!fReadMC) {
 296       Printf("ERROR: fUseMCVertex set without fReadMC set!");
 297       return;
 298     }
 299
 300     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 301     if (!eventHandler) {
 302       Printf("ERROR: Could not retrieve MC event handler");
 303       return;
 304     }
 305
 306     AliMCEvent* mcEvent = eventHandler->MCEvent();
 307     if (!mcEvent) {
 308       Printf("ERROR: Could not retrieve MC event");
 309       return;
 310     }
 311
 312     AliHeader* header = mcEvent->Header();
 313     if (!header)
 314     {
 315       AliDebug(AliLog::kError, "Header not available");
 316       return;
 317     }
 318
 319     // get the MC vertex
 320     AliGenEventHeader* genHeader = header->GenEventHeader();
 321     TArrayF vtxMC(3);
 322     genHeader->PrimaryVertex(vtxMC);
 323
 324     vtx[2] = vtxMC[2];
 325   }
 326
 327   // Processing of ESD information (always)
 328   if (eventTriggered)
 329   {
 330     // control hist
 331     fMult->Fill(inputCount);
 332     if (vtxESD)
 333     {
 334       // control hist
 335       fMultVtx->Fill(inputCount);
 336
 337       for (Int_t i=0; i<inputCount; ++i)
 338       {
 339         Float_t eta = etaArr[i];
 340         Float_t pt  = ptArr[i];
 341
 342         fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
 343
 344         if (TMath::Abs(vtx[2]) < 20)
 345         {
 346           fPartEta[0]->Fill(eta);
 347
 348           if (vtx[2] < 0)
 349             fPartEta[1]->Fill(eta);
 350           else
 351             fPartEta[2]->Fill(eta);
 352         }
 353       }
 354
 355       // for event count per vertex
 356       fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
 357
 358       // control hist
 359       fEvents->Fill(vtx[2]);
 360     }
 361   }
 362
 363   if (fReadMC)   // Processing of MC information (optional)
 364   {
 365     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 366     if (!eventHandler) {
 367       Printf("ERROR: Could not retrieve MC event handler");
 368       return;
 369     }
 370
 371     AliMCEvent* mcEvent = eventHandler->MCEvent();
 372     if (!mcEvent) {
 373       Printf("ERROR: Could not retrieve MC event");
 374       return;
 375     }
 376
 377     AliStack* stack = mcEvent->Stack();
 378     if (!stack)
 379     {
 380       AliDebug(AliLog::kError, "Stack not available");
 381       return;
 382     }
 383
 384     AliHeader* header = mcEvent->Header();
 385     if (!header)
 386     {
 387       AliDebug(AliLog::kError, "Header not available");
 388       return;
 389     }
 390
 391     // get the MC vertex
 392     AliGenEventHeader* genHeader = header->GenEventHeader();
 393     TArrayF vtxMC(3);
 394     genHeader->PrimaryVertex(vtxMC);
 395
 396     // loop over mc particles
 397     Int_t nPrim  = stack->GetNprimary();
 398
 399     Int_t nAcceptedParticles = 0;
 400
 401     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 402     {
 403       //Printf("Getting particle %d", iMc);
 404       TParticle* particle = stack->Particle(iMc);
 405
 406       if (!particle)
 407         continue;
 408
 409       if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 410         continue;
 411
 412       AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
 413       Float_t eta = particle->Eta();
 414       Float_t pt = particle->Pt();
 415
 416        // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
 417       if (TMath::Abs(eta) < 1.5 && pt > 0.3)
 418         nAcceptedParticles++;
 419
 420       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
 421       fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
 422
 423       if (eventTriggered)
 424       {
 425         fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
 426         if (vtxESD)
 427           fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
 428       }
 429
 430       if (TMath::Abs(eta) < 0.8)
 431         fPartPt->Fill(pt);
 432     }
 433
 434     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
 435     if (eventTriggered)
 436     {
 437       fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
 438       if (vtxESD)
 439         fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
 440     }
 441
 442     if (eventTriggered && vtxESD)
 443     {
 444       // from tracks is only done for triggered and vertex reconstructed events
 445
 446       for (Int_t i=0; i<inputCount; ++i)
 447       {
 448         Int_t label = labelArr[i];
 449
 450         if (label < 0)
 451           continue;
 452
 453         //Printf("Getting particle of track %d", label);
 454         TParticle* particle = stack->Particle(label);
 455
 456         if (!particle)
 457         {
 458           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
 459           continue;
 460         }
 461
 462         fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
 463       } // end of track loop
 464
 465       // for event count per vertex
 466       fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
 467     }
 468   }
 469
 470   delete[] etaArr;
 471   delete[] labelArr;
 472   delete[] ptArr;
 473 }
 474
 475 void AlidNdEtaTask::Terminate(Option_t *)
 476 {
 477   // The Terminate() function is the last function to be called during
 478   // a query. It always runs on the client, it can be used to present
 479   // the results graphically or save the results to file.
 480
 481   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 482   if (!fOutput) {
 483     Printf("ERROR: fOutput not available");
 484     return;
 485   }
 486
 487   fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
 488   fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
 489   fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
 490   fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
 491   fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
 492   for (Int_t i=0; i<3; ++i)
 493     fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
 494
 495   if (!fdNdEtaAnalysisESD)
 496   {
 497     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
 498     return;
 499   }
 500
 501   if (fMult && fMultVtx)
 502   {
 503     new TCanvas;
 504     fMult->Draw();
 505     fMultVtx->SetLineColor(2);
 506     fMultVtx->Draw("SAME");
 507   }
 508
 509   if (fPartEta[0])
 510   {
 511     Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
 512     Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
 513
 514     Printf("%d events with vertex used", events1 + events2);
 515
 516     if (events1 > 0 && events2 > 0)
 517     {
 518     fPartEta[0]->Scale(1.0 / (events1 + events2));
 519     fPartEta[1]->Scale(1.0 / events1);
 520     fPartEta[2]->Scale(1.0 / events2);
 521
 522     for (Int_t i=0; i<3; ++i)
 523       fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
 524
 525     new TCanvas("control", "control", 500, 500);
 526     for (Int_t i=0; i<3; ++i)
 527     {
 528       fPartEta[i]->SetLineColor(i+1);
 529       fPartEta[i]->Draw((i==0) ? "" : "SAME");
 530     }
 531     }
 532   }
 533
 534   if (fEvents)
 535   {
 536     new TCanvas("control3", "control3", 500, 500);
 537     fEvents->Draw();
 538   }
 539
 540   TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
 541
 542   if (fdNdEtaAnalysisESD)
 543     fdNdEtaAnalysisESD->SaveHistograms();
 544
 545   if (fEsdTrackCuts)
 546     fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
 547
 548   if (fMult)
 549     fMult->Write();
 550
 551   if (fMultVtx)
 552     fMultVtx->Write();
 553
 554   for (Int_t i=0; i<3; ++i)
 555     if (fPartEta[i])
 556       fPartEta[i]->Write();
 557
 558   if (fEvents)
 559     fEvents->Write();
 560
 561   if (fVertexResolution)
 562     fVertexResolution->Write();
 563
 564   if (fDeltaPhi)
 565     fDeltaPhi->Write();
 566
 567   fout->Write();
 568   fout->Close();
 569
 570   Printf("Writting result to analysis_esd_raw.root");
 571
 572   if (fReadMC)
 573   {
 574     fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
 575     fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
 576     fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
 577     fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
 578     fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
 579
 580     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
 581     {
 582       AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
 583       return;
 584     }
 585
 586     fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
 587     fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
 588     fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
 589     fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
 590
 591     Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
 592     fPartPt->Scale(1.0/events);
 593     fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
 594
 595     TFile* fout = new TFile("analysis_mc.root","RECREATE");
 596
 597     fdNdEtaAnalysis->SaveHistograms();
 598     fdNdEtaAnalysisTr->SaveHistograms();
 599     fdNdEtaAnalysisTrVtx->SaveHistograms();
 600     fdNdEtaAnalysisTracks->SaveHistograms();
 601     fPartPt->Write();
 602
 603     fout->Write();
 604     fout->Close();
 605
 606     Printf("Writting result to analysis_mc.root");
 607
 608     if (fPartPt)
 609     {
 610       new TCanvas("control2", "control2", 500, 500);
 611       fPartPt->Draw();
 612     }
 613   }
 614 }