| | 320 | Let's assume the output file is ''crab-data-only.root'' ({{{rootifysql --out=crab-data-only.root}}}). Requesting the data and writing the file took me about 60s. |
| | 321 | |
| | 322 | To run an analysis on the data you can use the following root macro "''ganymed.C''". It produces a theta-square plot. Its execution took about 5s ({{{root ganymed.C++}}}) |
| | 323 | |
| | 324 | {{{ |
| | 325 | #include <iostream> |
| | 326 | |
| | 327 | #include <TMath.h> |
| | 328 | #include <TH1.h> |
| | 329 | #include <TChain.h> |
| | 330 | #include <TStopwatch.h> |
| | 331 | |
| | 332 | void ganymed() |
| | 333 | { |
| | 334 | // Create chain for the tree Result |
| | 335 | // This is just easier than using TFile/TTree |
| | 336 | TChain c("Result"); |
| | 337 | |
| | 338 | // Add the input file to the |
| | 339 | c.AddFile("crab-data-only.root"); |
| | 340 | |
| | 341 | // Define variables for all leaves to be accessed |
| | 342 | // By definition rootifysql writes only doubles |
| | 343 | double X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta, |
| | 344 | M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size, |
| | 345 | ConcCore, ConcCOG, NumIslands, NumUsedPixels; |
| | 346 | |
| | 347 | // Connect the variables to the cordesponding leaves |
| | 348 | c.SetBranchAddress("X", &X); |
| | 349 | c.SetBranchAddress("Y", &Y); |
| | 350 | c.SetBranchAddress("MeanX", &MeanX); |
| | 351 | c.SetBranchAddress("MeanY", &MeanY); |
| | 352 | c.SetBranchAddress("Width", &Width); |
| | 353 | c.SetBranchAddress("Length", &Length); |
| | 354 | c.SetBranchAddress("CosDelta", &CosDelta); |
| | 355 | c.SetBranchAddress("SinDelta", &SinDelta); |
| | 356 | c.SetBranchAddress("M3Long", &M3Long); |
| | 357 | c.SetBranchAddress("SlopeLong", &SlopeLong); |
| | 358 | c.SetBranchAddress("Leakage1", &Leakage1); |
| | 359 | c.SetBranchAddress("NumIslands", &NumIslands); |
| | 360 | c.SetBranchAddress("NumUsedPixels", &NumUsedPixels); |
| | 361 | c.SetBranchAddress("SlopeSpreadWeighted", &SlopeSpreadWeighted); |
| | 362 | c.SetBranchAddress("Size", &Size); |
| | 363 | c.SetBranchAddress("ConcCOG", &ConcCOG); |
| | 364 | c.SetBranchAddress("ConcCore", &ConcCore); |
| | 365 | |
| | 366 | // Set some constants (they could be included in the database |
| | 367 | // in the future) |
| | 368 | double mm2deg = +0.0117193246260285378; |
| | 369 | double abberation = 1.02; |
| | 370 | |
| | 371 | // -------------------- Source dependent parameter calculation ------------------- |
| | 372 | |
| | 373 | // Create a histogram for on- and off-data |
| | 374 | TH1F hon("on", "", 55, 0, 1); |
| | 375 | TH1F hoff("off", "", 55, 0, 1); |
| | 376 | |
| | 377 | // Loop over all events |
| | 378 | TStopwatch clock; |
| | 379 | for (int i=0; i<c.GetEntries(); i++) |
| | 380 | { |
| | 381 | // read the i-th event from the file |
| | 382 | c.GetEntry(i); |
| | 383 | |
| | 384 | // First calculate all cuts to speedup the analysis |
| | 385 | double area = TMath::Pi()*Width*Length; |
| | 386 | |
| | 387 | bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1<0.1; |
| | 388 | |
| | 389 | bool cut0 = |
| | 390 | SlopeSpreadWeighted<0.68 && |
| | 391 | SlopeSpreadWeighted>0.18 && |
| | 392 | log10(area*mm2deg*mm2deg)<(log10(Size)-2)*1.1-1.55 && |
| | 393 | ConcCore>0.13 && |
| | 394 | ConcCOG >0.15; |
| | 395 | |
| | 396 | if (!cutq || !cut0) |
| | 397 | continue; |
| | 398 | |
| | 399 | // Loop over all wobble positions in the camera |
| | 400 | for (int angle=0; angle<360; angle+=60) |
| | 401 | { |
| | 402 | // -------------------- Source dependent parameter calculation ------------------- |
| | 403 | |
| | 404 | double cr = cos(angle*TMath::DegToRad()); |
| | 405 | double sr = sin(angle*TMath::DegToRad()); |
| | 406 | |
| | 407 | double px = cr*X-sr*Y; |
| | 408 | double py = cr*Y+sr*X; |
| | 409 | |
| | 410 | double dx = MeanX - px/abberation; |
| | 411 | double dy = MeanY - py/abberation; |
| | 412 | |
| | 413 | double dist = sqrt(dx*dx + dy*dy); |
| | 414 | |
| | 415 | double alpha = asin((CosDelta*dy - SinDelta*dx)/dist); |
| | 416 | double sgn = TMath::Sign(1., (CosDelta*dx + SinDelta*dy)/dist); |
| | 417 | |
| | 418 | // ------------------------------- Application ---------------------------------- |
| | 419 | |
| | 420 | double m3l = M3Long*sgn*mm2deg; |
| | 421 | double slope = SlopeLong*sgn/mm2deg; |
| | 422 | |
| | 423 | // --------------------------------- Aanlysis ----------------------------------- |
| | 424 | |
| | 425 | double xi = 1.39 + 0.154*slope + 1.679*(1-1/(1+4.86*Leakage1)); |
| | 426 | |
| | 427 | double sign1 = m3l+0.07; |
| | 428 | double sign2 = (dist*mm2deg-0.5)*7.2-slope; |
| | 429 | |
| | 430 | double disp = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length)/mm2deg; |
| | 431 | |
| | 432 | double thetasq = (disp*disp + dist*dist - 2*disp*dist*cos(alpha))*mm2deg*mm2deg; |
| | 433 | |
| | 434 | // Fill the on- and off-histograms |
| | 435 | if (angle==0) |
| | 436 | hon.Fill(thetasq); |
| | 437 | else |
| | 438 | hoff.Fill(thetasq, 1./5); |
| | 439 | } |
| | 440 | } |
| | 441 | clock.Print(); |
| | 442 | |
| | 443 | // Plot the result |
| | 444 | hon.SetMinimum(0); |
| | 445 | hon.DrawCopy(); |
| | 446 | hoff.DrawCopy("same"); |
| | 447 | } |
| | 448 | }}} |
| | 449 | |