TPC/CalibMacros/CalibLaserVscan.C

   1
   2 /*
   3    //
   4    // 0. Make a calibration
   5    // 1. Make a laser scan list
   6    //    e.g in TPC workscape
   7    find `pwd`/*/laserMean.root >laserScan.txt
   8    // 2. Define a reference data
   9    rrunA=84469/; for a in `cat laserScan.txt`; do echo `pwd`/$rrunA/laserMean.root; done >laserScanRefA.txt
  10    rrunC=84469; for a in `cat laserScan.txt`; do echo `pwd`/$rrunC/laserMean.root; done >laserScanRefC.txt
  11    rrun=84469; for a in `cat laserScan.txt`; do echo `pwd`/$rrun/laserMean.root; done >laserScanRef.txt
  12                                                                                   //
  13    //
  14    .x ~/rootlogon.C
  15    gSystem->Load("libANALYSIS");
  16    gSystem->Load("libTPCcalib");
  17    gSystem->Load("libSTAT");
  18
  19    gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  20    gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  21    .L $ALICE_ROOT/TPC/CalibMacros/CalibLaserVscan.C+
  22
  23    AliXRDPROOFtoolkit tool;
  24    chain = tool.MakeChainRandom("laserScan.txt","Mean",0,10200);
  25    chain->Lookup();
  26    chainRef = tool.MakeChain("laserScanRef.txt","Mean",0,10200);
  27    chain->AddFriend(chainRef,"R")
  28    chainRefA = tool.MakeChain("laserScanRefA.txt","Mean",0,10200);
  29    chain->AddFriend(chainRefA,"RA")
  30    chainRefC = tool.MakeChain("laserScanRefC.txt","Mean",0,10200);
  31    chain->AddFriend(chainRefC,"RC")
  32    //
  33    // MakeMeanBundle();
  34    // SaveResult();
  35    //
  36    ReadRunSetup();
  37    ReadResult();
  38
  39    MakeAnalysisBeam();
  40    TFile fbundle("scanDeltaBeam.root");
  41    chain->Draw("mphi:GetValueBundle(id,1)","isOK&&GetValueBundle(id,0)>3&&LTr.fSide==0","")
  42 */
  43
  44 #include <fstream>
  45 #include "TFile.h"
  46 #include "TMatrixD.h"
  47 #include "TVectorD.h"
  48
  49 #include "TChain.h"
  50 #include "TCut.h"
  51 #include "TH1F.h"
  52 #include "TObjArray.h"
  53 #include "TProfile.h"
  54 #include "TCanvas.h"
  55 #include "TH1F.h"
  56 #include "TGraph.h"
  57 #include "TF1.h"
  58 #include "TTreeStream.h"
  59 #include "TLegend.h"
  60 #include "AliTPCLaserTrack.h"
  61 #include "AliTPCcalibDB.h"
  62 #include "TStatToolkit.h"
  63
  64
  65 TMatrixD matrixP4Corr(7,2);       // P4 correction for P0 - y position and P2 -snp(phi)
  66 TMatrixD matrixP2RMSIROC(7,5);    // rms of second derivative IROC
  67 TMatrixD matrixP2RMSOROC(7,5);    // rms of second derivative OROC
  68 // indexes - beam rod
  69 //           rod 5 means all rods
  70 //
  71 TMatrixD matrixMeanBundle(336,6); // mean position/delta per bundle
  72 TMatrixD matrixRMSBundle(336,6);  // rms of distribution
  73 // indexes:
  74 // 0 - number of entries
  75 // 1 - dy
  76 // 2 - dphi
  77 // 3 - dp4
  78 // 4 - dy0
  79 // 5 - p0:p4 slope
  80 //
  81 map<int,TVectorD*> mapRunVoltage;             // run to voltage
  82 // run#  -0 ggA 1- ggC 2- coA 3-coC 4-skA 5-skC
  83 map<int,int> mapRunVgg;           // run to gating grid voltage
  84 map<int,int> mapRunVskirt;        //
  85 map<int,int> mapRunVcover;        //
  86 TArrayI runlist(10000);
  87 //
  88 TChain * chain=0;
  89 TObjArray apic;
  90 //
  91 //
  92 TCut tA="isOK&&RA.isOK";
  93 TCut tC="isOK&&RC.isOK";
  94 TCut cA="eY.fElements<0.01&&RA.eY.fElements<0.01&&X.fElements>10&&RA.X.fElements>10";
  95 TCut cC="eY.fElements<0.01&&RA.eY.fElements<0.01&&X.fElements>10&&RA.X.fElements>10";
  96
  97
  98
  99
 100 Double_t GetValueBundle(Int_t id, Int_t type){
 101   //
 102   //
 103   return matrixMeanBundle(id,type);
 104 }
 105
 106 Double_t GetP4Corr(Int_t id, Int_t value){
 107   AliTPCLaserTrack *ltrack =(AliTPCLaserTrack *) AliTPCLaserTrack::GetTracks()->At(id);
 108   Int_t beam = ltrack->GetBeam();
 109   return matrixP4Corr(beam,value);
 110 }
 111
 112 Double_t GetDyBundle(Int_t id){
 113   Double_t dy = matrixMeanBundle(id,4);
 114   Double_t dx = matrixMeanBundle(id,5);
 115   AliTPCLaserTrack *ltrack =(AliTPCLaserTrack *) AliTPCLaserTrack::GetTracks()->At(id);
 116   Double_t p2 = ltrack->GetParameter()[2];
 117   Double_t delta = dy+dx*p2;
 118   return delta;
 119 }
 120
 121 Double_t GetRMSBundle(Int_t id, Int_t type){
 122   //
 123   //
 124   return matrixRMSBundle(id,type);
 125 }
 126 Int_t GetVoltage(Int_t run, Int_t type){
 127   //
 128   // Get the voltage
 129   // run#  -0 ggA 1- ggC 2- coA 3-coC 4-skA 5-skC
 130   TVectorD *runVoltage = mapRunVoltage[run];
 131   if (!runVoltage) return -1;
 132   return (*runVoltage)[type];
 133 }
 134
 135
 136
 137 void SaveResult(){
 138   TFile f("laserData.root","recreate");
 139   matrixMeanBundle.Write("matrixMeanBundle");
 140   matrixRMSBundle.Write("matrixRMSBundle");
 141   matrixP4Corr.Write("matrixP4Corr");
 142   matrixP2RMSIROC.Write("matrixP2RMSIROC");
 143   matrixP2RMSOROC.Write("matrixP2RMSOROC");
 144   f.Close();
 145 }
 146 void ReadResult(){
 147   AliTPCLaserTrack::LoadTracks();
 148   TFile f("laserData.root");
 149   matrixMeanBundle.Read("matrixMeanBundle");
 150   matrixRMSBundle.Read("matrixRMSBundle");
 151   matrixP4Corr.Read("matrixP4Corr");
 152   matrixP2RMSIROC.Read("matrixP2RMSIROC");
 153   matrixP2RMSOROC.Read("matrixP2RMSOROC");
 154   f.Close();
 155 }
 156
 157 void ReadRunSetup(){
 158   ifstream in;
 159   in.open("runSetup.txt");
 160   TString objfile;
 161   string line;
 162   TObjArray *arr = 0;
 163   Int_t counter=0;
 164
 165   while(in.good()) {
 166     //in >> objfile;
 167     getline(in,line);
 168     objfile=line;
 169     printf("%s\n",objfile.Data());
 170     arr = objfile.Tokenize(" ");
 171     if (!arr) continue;
 172     if (arr->GetEntries()>=7){
 173       Int_t run = atoi(arr->At(0)->GetName());
 174       Int_t vcover = atoi(arr->At(0)->GetName());
 175       Int_t vskirt = atoi(arr->At(1)->GetName());
 176       Int_t vgg = atoi(arr->At(3)->GetName());
 177       TVectorD * vsetup = new TVectorD(6);
 178       for (Int_t iv=0; iv<6; iv++){
 179         (*vsetup)[iv] = atoi(arr->At(iv+1)->GetName());
 180       }
 181       mapRunVoltage[run]=vsetup;
 182       mapRunVgg[run]=vgg;
 183       mapRunVskirt[run]=vskirt;
 184       mapRunVcover[run]=vcover;
 185       runlist[counter]=run;
 186       counter++;
 187     }
 188     delete arr;
 189   }
 190   runlist[counter]=-1;
 191 }
 192
 193
 194
 195
 196
 197 void MakeAnalysisBeam(){
 198   //
 199   //
 200   //
 201   TTreeSRedirector *pcstream = new TTreeSRedirector("scanDeltaBeam.root");
 202   TH1* phisP0=new TH1F("hhisP0","hhisP0",100,-5.,5.);
 203   TH1* phisP0X=new TH1F("hhisP0X","hhisP0X",100,-5.,5.);
 204   TH1* phisP2=new TH1F("hhisP2","hhisP2",100,-6,6);
 205   TH1* phisP4=new TH1F("hhisP4","hhisP4",100,-0.05,0.05);
 206   // second derivative
 207   TH1* phisP2IROC=new TH1F("hhisP2IROC","hhisPIROC",100,-2.,2.);
 208   TH1* phisP2OROC=new TH1F("hhisP2OROC","hhisPOROC",100,-2.,2.);
 209   //
 210   //
 211   //
 212   for (Int_t iside=0;iside<2;iside++){
 213     for (Int_t ibeam=0; ibeam<8; ibeam++){
 214       for (Int_t irun=0; irun<100; irun++){
 215         Int_t run =runlist[irun];
 216         if (run==-1) break;
 217         if (run==0) continue;
 218         Int_t vskirt = GetVoltage(run,1);
 219         Int_t vcover = GetVoltage(run,2);
 220         Int_t vgg    = GetVoltage(run,0);
 221         //
 222         TCut cut = Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&LTr.fBeam==%d&&run==%d", iside, ibeam,run);
 223         if (ibeam==7){
 224           cut = Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&run==%d", iside,run);
 225         }
 226         Int_t entries = chain->Draw("gp41-GetValueBundle(id,3)>>hhisP4",cut,"goff");
 227         if (entries==0) continue;
 228         chain->Draw("10*(mphi-GetValueBundle(id,1)-GetP4Corr(id,0)*gp41)>>hhisP0",cut,"goff");
 229         chain->Draw("10*(mphi-GetDyBundle(id)-GetP4Corr(id,0)*gp41)>>hhisP0X",cut,"goff");
 230         chain->Draw("1000*(mphiP-GetValueBundle(id,2)-GetP4Corr(id,1)*gp41)>>hhisP2",cut,"goff");
 231         //
 232         printf("%d\t%d\t%d\t%f\t%f\t%d\n", iside, ibeam,run, phisP0X->GetMean(), phisP0X->GetRMS(), Int_t(phisP0X->GetEntries()));
 233         //
 234         chain->Draw("10*(mphi-GetValueBundle(id,1)-GetP4Corr(id,0)*gp41)>>hhisP0",cut,"goff");
 235         //
 236         chain->Draw("mPy2vP2In*100^2>>hhisP2IROC",cut,"goff");
 237         chain->Draw("mPy2vP2Out*100^2>>hhisP2OROC",cut,"goff");
 238         //
 239         //
 240         Double_t mp0  = phisP0->GetMean();
 241         Double_t mp0X = phisP0X->GetMean();
 242         Double_t mp2  = phisP2->GetMean();
 243         Double_t mp4  = phisP4->GetMean();
 244         Double_t mp2I = phisP2IROC->GetMean();
 245         Double_t mp2O = phisP2OROC->GetMean();
 246         Double_t sp0  = phisP0->GetRMS();
 247         Double_t sp0X = phisP0X->GetRMS();
 248         Double_t sp2  = phisP2->GetRMS();
 249         Double_t sp4  = phisP4->GetRMS();
 250         Double_t sp2I = phisP2IROC->GetRMS();
 251         Double_t sp2O = phisP2OROC->GetRMS();
 252         (*pcstream)<<"vScanBeam"<<
 253           "side="<<iside<<
 254           "run="<<run<<
 255           "ibeam="<<ibeam<<
 256           "vgg="<<vgg<<
 257           "vskirt="<<vskirt<<
 258           "vcover="<<vcover<<
 259           "entries="<<entries<<
 260           "mp0="<<mp0<<
 261           "mp0X="<<mp0X<<
 262           "mp2="<<mp2<<
 263           "mp4="<<mp4<<
 264           "mp2I="<<mp2I<<
 265           "mp2O="<<mp2O<<
 266           "sp0="<<sp0<<
 267           "sp0X="<<sp0X<<
 268           "sp2="<<sp2<<
 269           "sp4="<<sp4<<
 270           "sp2I="<<sp2I<<
 271           "sp2O="<<sp2O<<
 272           "\n";
 273       }
 274     }
 275   }
 276   delete pcstream;
 277 }
 278
 279
 280
 281
 282 void MakeMeanBundle(){
 283   //
 284   //
 285   //
 286   AliTPCLaserTrack::LoadTracks();
 287   AliTPCLaserTrack *ltrack;
 288   TF1 * fp1 = 0;
 289   TCut ccut;
 290   TH1F * phisP0 = 0;
 291   TH1F * phisP2 = 0;
 292   TH1F * phisP4 = 0;
 293   //
 294   // get p4 corr
 295   //
 296   for (Int_t ibeam=0; ibeam<7; ibeam++){
 297     Int_t entries0=chain->Draw("mphi:gp41",Form("isOK&&LTr.fBeam==%d",ibeam));
 298     TGraph gr0(entries0,chain->GetV2(),chain->GetV1());
 299     gr0.Draw();
 300     gr0.Fit("pol1","Q","Q");
 301     fp1 = gr0.GetFunction("pol1");
 302     matrixP4Corr(ibeam,0)=fp1->GetParameter(1);
 303     Int_t entries2=chain->Draw("mphiP:gp41",Form("isOK&&LTr.fBeam==%d",ibeam));
 304     TGraph gr2(entries2,chain->GetV2(),chain->GetV1());
 305     gr2.Draw();
 306     gr2.Fit("pol1","Q","Q");
 307     fp1 = gr2.GetFunction("pol1");
 308     matrixP4Corr(ibeam,1)=fp1->GetParameter(1);
 309   }
 310   //
 311   // get RMS of second derivative
 312   //
 313   phisP0= new TH1F("hisP0","hisP0",500,-1.,1.);
 314   for (Int_t irod=0; irod<5;irod++)
 315     for (Int_t ibeam=0; ibeam<7; ibeam++){
 316       //
 317       TCut cut=Form("isOK&&LTr.fBeam==%d&&LTr.fRod==%d",ibeam,irod);
 318       if (irod==5) cut=Form("isOK&&LTr.fBeam==%d",ibeam);
 319       chain->Draw("mPy2vP2In*100^2>>hisP0",cut);
 320       matrixP2RMSIROC(ibeam,irod) = phisP0->GetRMS();
 321       chain->Draw("mPy2vP2Out*100^2>>hisP0",cut);
 322       matrixP2RMSOROC(ibeam,irod) = phisP0->GetRMS();
 323   }
 324   delete phisP0;
 325
 326   //
 327   // Get mean bundle position
 328   //
 329   for (Int_t id=0; id<336; id+=7){
 330     ltrack =(AliTPCLaserTrack *) AliTPCLaserTrack::GetTracks()->At(id);
 331     Int_t side   = ltrack->GetSide();
 332     Int_t rod    = ltrack->GetRod();
 333     Int_t bundle = ltrack->GetBundle();
 334     //    Int_t beam   = ltrack->GetBeam();
 335     ccut = Form("isOK&&LTr.fSide==%d&&LTr.fBundle==%d&&LTr.fRod==%d",side,bundle,rod);
 336     phisP0= new TH1F("hisP0","hisP0",500,-0.5,0.5);
 337     phisP2= new TH1F("hisP2","hisP2",500,-0.02,0.02);
 338     phisP4= new TH1F("hisP4","hisP4",500,-0.2,0.2);
 339     Int_t entries =chain->Draw("mphi",Form("isOK&&id==%d",id),"goff");
 340     chain->Draw("mphi-GetP4Corr(id,0)*gp41>>hisP0",ccut,"");
 341     chain->Draw("mphiP-GetP4Corr(id,1)*gp41>>hisP2",ccut,"");
 342     chain->Draw("gp41>>hisP4",ccut,"");
 343     Int_t entriesG =chain->Draw("mphi-GetP4Corr(id,0)*gp41:tan(asin(LTr.fP[2]))",ccut,"");
 344     if (entriesG<3) continue;
 345     TGraph gr(entriesG,chain->GetV2(),chain->GetV1());
 346     gr.Draw();
 347     gr.Fit("pol1","Q","Q");
 348     fp1 = gr.GetFunction("pol1");
 349     for (Int_t id2=0;id2<6; id2++){
 350       Int_t jd = id+id2;
 351       matrixMeanBundle(jd,0) = entries;
 352       matrixMeanBundle(jd,1) =phisP0->GetMean();
 353       matrixRMSBundle(jd,1) =phisP0->GetRMS();
 354       matrixMeanBundle(jd,2) =phisP2->GetMean();
 355       matrixRMSBundle(jd,2) =phisP2->GetRMS();
 356       matrixMeanBundle(jd,3) =phisP4->GetMean();
 357       matrixRMSBundle(jd,3) =phisP4->GetRMS();
 358       matrixMeanBundle(jd,4) = fp1->GetParameter(0);  // delta y
 359       matrixMeanBundle(jd,5) = fp1->GetParameter(1);  // delta x
 360     }
 361     //
 362     printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\n",id,entries,matrixMeanBundle(id,1),matrixMeanBundle(id,2),matrixMeanBundle(id,3),matrixMeanBundle(id,4),matrixMeanBundle(id,5));
 363     printf("%d\t%d\t%f\t%f\t%f\n",id,entries,matrixRMSBundle(id,1),matrixRMSBundle(id,2),matrixRMSBundle(id,3));
 364     delete phisP0;
 365     delete phisP2;
 366     delete phisP4;
 367   }
 368 }
 369
 370
 371
 372
 373
 374
 375
 376
 377 void MakeGraphsdY(){
 378   //
 379   // Make delta Y pictures from voltage scan
 380   //
 381   TObjArray *aprofY = new TObjArray(14);
 382   for (Int_t ib=0;ib<14;ib++){
 383     TProfile *profY = new TProfile("py","py",100,0,150);
 384     chain->Draw("10*(mphi-GetValueBundle(id,1)):bz>>py",Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&LTr.fBeam==%d",ib/7,ib%7),"prof");
 385     profY->SetName(Form("#DeltaY side %d beam %d",ib/7,ib%7));
 386     profY->SetTitle(Form("#DeltaY side %d beam %d",ib/7,ib%7));
 387     aprofY->AddAt(profY,ib);
 388     profY->SetMaximum(2.5);
 389     profY->SetMinimum(-2.5);
 390     profY->SetMarkerColor((ib%7)+1);
 391     profY->SetMarkerStyle((ib%7)+22);
 392     profY->SetMarkerSize(1.2);
 393     profY->SetXTitle("U_{gg} (V)");
 394     profY->SetYTitle("#Delta_{y} (mm)");
 395     profY->Fit("pol1");
 396   }
 397   TCanvas *cY = new TCanvas("deltaY","deltaY",900,600);
 398   cY->Divide(5,3);
 399   cY->Draw();
 400   {
 401     for (Int_t ib=0;ib<14;ib++){
 402       cY->cd(ib+1);
 403       aprofY->At(ib)->Draw("p");
 404     }
 405   }
 406   cY->SaveAs("pic/deltaYlaserVscna.eps");
 407   cY->SaveAs("pic/deltaYlaserVscna.gif");
 408   apic.AddLast(cY);
 409 }
 410
 411 void MakeGraphsdP2(){
 412   //
 413   // Make delta Y pictures from voltage scan
 414   //
 415   TObjArray *aprofP2 = new TObjArray(14);
 416   for (Int_t ib=0;ib<14;ib++){
 417     TProfile *profP2 = new TProfile("pyP","pyP",100,0,150);
 418     chain->Draw("(mphiP-GetValueBundle(id,2)):bz>>pyP",Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&LTr.fBeam==%d",ib/7,ib%7),"prof");
 419     profP2->SetName(Form("#Delta_{#phi} side %d beam %d",ib/7,ib%7));
 420     profP2->SetTitle(Form("#Delta_{#phi} side %d beam %d",ib/7,ib%7));
 421     aprofP2->AddAt(profP2,ib);
 422     profP2->SetMaximum(0.005);
 423     profP2->SetMinimum(-0.005);
 424     profP2->SetMarkerColor((ib%7)+1);
 425     profP2->SetMarkerStyle((ib%7)+22);
 426     profP2->SetMarkerSize(1.2);
 427     profP2->SetXTitle("U_{gg} (V)");
 428     profP2->SetYTitle("#Delta_{#phi}");
 429     profP2->Fit("pol1");
 430   }
 431   TCanvas *cP2 = new TCanvas("deltaP2","deltaP2",900,600);
 432   cP2->Divide(5,3);
 433   cP2->Draw();
 434   {
 435     for (Int_t ib=0;ib<14;ib++){
 436       cP2->cd(ib+1);
 437       aprofP2->At(ib)->Draw("p");
 438     }
 439   }
 440   cP2->SaveAs("pic/deltaP2laserVscna.eps");
 441   cP2->SaveAs("pic/deltaP2laserVscna.gif");
 442   apic.AddLast(cP2);
 443 }
 444
 445
 446
 447
 448 void MakeGraphsP4(){
 449   //
 450   // Make delta Y pictures from voltage scan
 451   //
 452   TObjArray *aprofP4 = new TObjArray(14);
 453   for (Int_t ib=0;ib<14;ib++){
 454     TProfile *profP4 = new TProfile("pp4","pp4",100,0,150);
 455     chain->Draw("(gp41-GetValueBundle(id,3)):bz>>pp4",Form("isOK&&GetValueBundle(id,0)>3&&LTr.fSide==%d&&LTr.fBeam==%d",ib/7,ib%7),"prof");
 456     profP4->SetName(Form("1/p_{t} side %d beam %d",ib/7,ib%7));
 457     profP4->SetTitle(Form("1/p_{t} side %d beam %d",ib/7,ib%7));
 458     aprofP4->AddAt(profP4,ib);
 459     profP4->SetMaximum(0.025);
 460     profP4->SetMinimum(-0.025);
 461     profP4->SetMarkerColor((ib%7)+1);
 462     profP4->SetMarkerStyle((ib%7)+22);
 463     profP4->SetMarkerSize(1.2);
 464     profP4->SetXTitle("U_{gg} (V)");
 465     profP4->SetYTitle("1/p_{t} (GeV/c)");
 466     profP4->Fit("pol1");
 467   }
 468   TCanvas *cY = new TCanvas("P4","P4",900,600);
 469   cY->Divide(5,3);
 470   cY->Draw();
 471   {
 472     for (Int_t ib=0;ib<14;ib++){
 473       cY->cd(ib+1);
 474       aprofP4->At(ib)->Draw("p");
 475     }
 476   }
 477   cY->SaveAs("pic/m1ptlaserVscna.eps");
 478   cY->SaveAs("pic/m1ptlaserVscna.gif");
 479   apic.AddLast(cY);
 480 }
 481
 482
 483 void MakePlotsP2GG(TCut ucut){
 484   //
 485   //
 486   //
 487   TFile fbundle("scanDeltaBeam.root");
 488   TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
 489   TGraph *graph[4];
 490   Int_t  mstyle[4]={22,23,24,25};
 491   Int_t  mcolor[4]={2,3,4,6};
 492   Int_t entries=0;
 493   entries = treeScan->Draw("10*sp2I/4:vgg","ibeam==7&&side==0"+ucut,"");
 494   graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 495   entries = treeScan->Draw("10*sp2O/4:vgg","ibeam==7&&side==0"+ucut,"");
 496   graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 497   entries = treeScan->Draw("10*sp2I/4:vgg","ibeam==7&&side==1"+ucut,"");
 498   graph[2]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 499   entries = treeScan->Draw("10*sp2O/4:vgg","ibeam==7&&side==1"+ucut,"");
 500   graph[3]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 501
 502   for (Int_t i=0; i<4; i++){
 503     graph[i]->GetXaxis()->SetTitle("U_{gg} (V)");
 504     graph[i]->GetYaxis()->SetTitle("Sagita (mm)");
 505     graph[i]->SetMinimum(0);
 506     graph[i]->SetMaximum(2.5);
 507     graph[i]->SetMarkerStyle(mstyle[i]);
 508     graph[i]->SetMarkerSize(1);
 509     graph[i]->SetMarkerColor(mcolor[i]);
 510     if (i==0) graph[i]->Draw("ap");
 511     graph[i]->Draw("p");
 512   }
 513   TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan - Sagita on 50 cm");
 514   legend->AddEntry(graph[0],"IROC A side");
 515   legend->AddEntry(graph[1],"OROC A side");
 516   legend->AddEntry(graph[2],"IROC C side");
 517   legend->AddEntry(graph[3],"OROC C side");
 518   legend->Draw();
 519   gPad->SaveAs("pic/scansagita_Vgg.eps");
 520   gPad->SaveAs("pic/scansagita_Vgg.gif");
 521 }
 522
 523
 524
 525
 526
 527 void MakePlotsP2Cover(TCut ucut){
 528   //
 529   //
 530   //
 531   TFile fbundle("scanDeltaBeam.root");
 532   TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
 533   TGraph *graph[4];
 534   Int_t  mstyle[4]={22,23,24,25};
 535   Int_t  mcolor[4]={2,3,4,6};
 536   Int_t entries=0;
 537   entries = treeScan->Draw("10*sp2I/4:vcover","ibeam==7&&side==0"+ucut,"");
 538   graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 539   entries = treeScan->Draw("10*sp2O/4:vcover","ibeam==7&&side==0"+ucut,"");
 540   graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 541   entries = treeScan->Draw("10*sp2I/4:vcover","ibeam==7&&side==1"+ucut,"");
 542   graph[2]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 543   entries = treeScan->Draw("10*sp2O/4:vcover","ibeam==7&&side==1"+ucut,"");
 544   graph[3]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 545
 546   for (Int_t i=0; i<4; i++){
 547     graph[i]->GetXaxis()->SetTitle("U_{cover} (V)");
 548     graph[i]->GetYaxis()->SetTitle("Sagita (mm)");
 549     graph[i]->SetMinimum(0.25);
 550     graph[i]->SetMaximum(0.7);
 551     graph[i]->SetMarkerStyle(mstyle[i]);
 552     graph[i]->SetMarkerSize(1);
 553     graph[i]->SetMarkerColor(mcolor[i]);
 554     if (i==0) graph[i]->Draw("ap");
 555     graph[i]->Draw("p");
 556   }
 557   TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan - Sagita on 50 cm");
 558   legend->AddEntry(graph[0],"IROC A side");
 559   legend->AddEntry(graph[1],"OROC A side");
 560   legend->AddEntry(graph[2],"IROC C side");
 561   legend->AddEntry(graph[3],"OROC C side");
 562   legend->Draw();
 563   gPad->SaveAs("pic/scansagita_Cover.eps");
 564   gPad->SaveAs("pic/scansagita_Cover.gif");
 565 }
 566
 567 void MakePlotsP2Skirt(TCut ucut){
 568   //
 569   //
 570   //
 571   TFile fbundle("scanDeltaBeam.root");
 572   TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
 573   TGraph *graph[4];
 574   Int_t  mstyle[4]={22,23,24,25};
 575   Int_t  mcolor[4]={2,3,4,6};
 576   Int_t entries=0;
 577   entries = treeScan->Draw("10*sp2I/4:vskirt","ibeam==7&&side==0"+ucut,"");
 578   graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 579   entries = treeScan->Draw("10*sp2O/4:vskirt","ibeam==7&&side==0"+ucut,"");
 580   graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 581   entries = treeScan->Draw("10*sp2I/4:vskirt","ibeam==7&&side==1"+ucut,"");
 582   graph[2]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 583   entries = treeScan->Draw("10*sp2O/4:vskirt","ibeam==7&&side==1"+ucut,"");
 584   graph[3]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 585
 586   for (Int_t i=0; i<4; i++){
 587     graph[i]->GetXaxis()->SetTitle("U_{skirt} (V)");
 588     graph[i]->GetYaxis()->SetTitle("Sagita (mm)");
 589     graph[i]->SetMinimum(0.25);
 590     graph[i]->SetMaximum(0.7);
 591     graph[i]->SetMarkerStyle(mstyle[i]);
 592     graph[i]->SetMarkerSize(1);
 593     graph[i]->SetMarkerColor(mcolor[i]);
 594     if (i==0) graph[i]->Draw("ap");
 595     graph[i]->Draw("p");
 596   }
 597   TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan - Sagita on 50 cm");
 598   legend->AddEntry(graph[0],"IROC A side");
 599   legend->AddEntry(graph[1],"OROC A side");
 600   legend->AddEntry(graph[2],"IROC C side");
 601   legend->AddEntry(graph[3],"OROC C side");
 602   legend->Draw();
 603   gPad->SaveAs("pic/scansagita_Skirt.eps");
 604   gPad->SaveAs("pic/scansagita_Skirt.gif");
 605 }
 606
 607
 608
 609
 610 void MakePlotsdYGG(){
 611   //
 612   //
 613   //
 614   TFile fbundle("scanDeltaBeam.root");
 615   TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
 616   TGraph *graph[4];
 617   Int_t  mstyle[4]={22,23,24,25};
 618   Int_t  mcolor[4]={2,3,4,6};
 619   Int_t entries=0;
 620   entries = treeScan->Draw("sp0X:vgg","ibeam==7&&side==0","");
 621   graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 622   entries = treeScan->Draw("sp0X:vgg","ibeam==7&&side==1","");
 623   graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 624
 625   for (Int_t i=0; i<2; i++){
 626     graph[i]->GetXaxis()->SetTitle("U_{gg} (V)");
 627     graph[i]->GetYaxis()->SetTitle("#sigma_{r#phi} (mm)");
 628     graph[i]->SetMinimum(0);
 629     graph[i]->SetMaximum(1);
 630     graph[i]->SetMarkerStyle(mstyle[i]);
 631     graph[i]->SetMarkerSize(1);
 632     graph[i]->SetMarkerColor(mcolor[i]);
 633     if (i==0) graph[i]->Draw("ap");
 634     graph[i]->Draw("p");
 635   }
 636   TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan-#sigma_{r#phi}");
 637   legend->AddEntry(graph[0],"A side");
 638   legend->AddEntry(graph[1],"C side");
 639   legend->Draw();
 640   gPad->SaveAs("pic/scandy_Vgg.eps");
 641   gPad->SaveAs("pic/scandy_Vgg.gif");
 642 }
 643
 644
 645 void MakePlotsdYCover(TCut ucut){
 646   //
 647   //
 648   //
 649   TFile fbundle("scanDeltaBeam.root");
 650   TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
 651   TGraph *graph[4];
 652   Int_t  mstyle[4]={22,23,24,25};
 653   Int_t  mcolor[4]={2,3,4,6};
 654   Int_t entries=0;
 655   entries = treeScan->Draw("sp0X:vcover","ibeam==7&&side==0"+ucut,"");
 656   graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 657   entries = treeScan->Draw("sp0X:vcover","ibeam==7&&side==1"+ucut,"");
 658   graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 659
 660   for (Int_t i=0; i<2; i++){
 661     graph[i]->GetXaxis()->SetTitle("U_{cover} (V)");
 662     graph[i]->GetYaxis()->SetTitle("#sigma_{r#phi} (mm)");
 663     graph[i]->SetMinimum(0);
 664     graph[i]->SetMaximum(0.75);
 665     graph[i]->SetMarkerStyle(mstyle[i]);
 666     graph[i]->SetMarkerSize(1);
 667     graph[i]->SetMarkerColor(mcolor[i]);
 668     if (i==0) graph[i]->Draw("ap");
 669     graph[i]->Draw("p");
 670   }
 671   TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan-#sigma_{r#phi}");
 672   legend->AddEntry(graph[0],"A side");
 673   legend->AddEntry(graph[1],"C side");
 674   legend->Draw();
 675   gPad->SaveAs("pic/scandy_Vcover.eps");
 676   gPad->SaveAs("pic/scandy_Vcover.gif");
 677 }
 678
 679 void MakePlotsdYSkirt(TCut ucut){
 680   //
 681   //
 682   //
 683   TFile fbundle("scanDeltaBeam.root");
 684   TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
 685   TGraph *graph[4];
 686   Int_t  mstyle[4]={22,23,24,25};
 687   Int_t  mcolor[4]={2,3,4,6};
 688   Int_t entries=0;
 689   entries = treeScan->Draw("sp0X:vskirt","ibeam==7&&side==0"+ucut,"");
 690   graph[0]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 691   entries = treeScan->Draw("sp0X:vskirt","ibeam==7&&side==1"+ucut,"");
 692   graph[1]= new TGraph(entries, treeScan->GetV2(), treeScan->GetV1());
 693
 694   for (Int_t i=0; i<2; i++){
 695     graph[i]->GetXaxis()->SetTitle("U_{skirt} (V)");
 696     graph[i]->GetYaxis()->SetTitle("#sigma_{r#phi} (mm)");
 697     graph[i]->SetMinimum(0);
 698     graph[i]->SetMaximum(0.75);
 699     graph[i]->SetMarkerStyle(mstyle[i]);
 700     graph[i]->SetMarkerSize(1);
 701     graph[i]->SetMarkerColor(mcolor[i]);
 702     if (i==0) graph[i]->Draw("ap");
 703     graph[i]->Draw("p");
 704   }
 705   TLegend *legend = new TLegend(0.45,0.70,0.9,0.9, "Voltage scan-#sigma_{r#phi}");
 706   legend->AddEntry(graph[0],"A side");
 707   legend->AddEntry(graph[1],"C side");
 708   legend->Draw();
 709   gPad->SaveAs("pic/scandy_Vskirt.eps");
 710   gPad->SaveAs("pic/scandy_Vskirt.gif");
 711 }
 712
 713
 714
 715 void GetOptimalSetting(){
 716   //
 717   //
 718   //
 719   TFile fbundle("scanDeltaBeam.root");
 720   TTree * treeScan = (TTree*)fbundle.Get("vScanBeam");
 721   //
 722   Float_t meansp0X=0;
 723   Float_t meansp2I=0;
 724   Float_t meansp2O=0;
 725   TH1F * his = new TH1F("hisDelta","hisDelta",300,0,1);
 726   treeScan->Draw("sp0X>>hisDelta","ibeam==7","");
 727   meansp0X = his->GetMean();
 728   treeScan->Draw("sp2I>>hisDelta","ibeam==7","");
 729   meansp2I = his->GetMean();
 730   treeScan->Draw("sp2O>>hisDelta","ibeam==7","");
 731   meansp2O = his->GetMean();
 732   //
 733   //
 734   //
 735   Int_t indexes[10000];
 736   for (Int_t iside=0;iside<2;iside++){
 737     printf("*************************************\n");
 738     printf("Side%d\tRun\tchi2\t\tVgg\tVskirt\tVcover\n",iside);
 739
 740     Int_t entries  =  treeScan->Draw(Form("sp0X/%f+sp2I/%f+sp2O/%f:run",meansp0X,meansp2I,meansp2O),Form("ibeam==7&&side==%d",iside));
 741     TMath::Sort(entries,treeScan->GetV1(),indexes,kFALSE);
 742     for (Int_t i=0; i<entries; i++){
 743       Int_t index = indexes[i];
 744       Int_t run   = treeScan->GetV2()[index];
 745       Float_t chi2=treeScan->GetV1()[index];
 746       Float_t vgg =
 747         printf("%d\t%d\t%f\t%d\t%d\t%d\n",iside,run, chi2,GetVoltage(run,0), GetVoltage(run,1),GetVoltage(run,2));
 748     }
 749   }
 750
 751    for (Int_t iside=0;iside<2;iside++){
 752     printf("***********************\n");
 753     printf("Side%d\tRun\tchi2\t\tVgg\tVskirt\tVcover\n",iside);
 754
 755     Int_t entries  =  treeScan->Draw(Form("sp2I/%f+sp2O/%f:run",meansp0X,meansp2I,meansp2O),Form("ibeam==7&&side==%d",iside));
 756     TMath::Sort(entries,treeScan->GetV1(),indexes,kFALSE);
 757     for (Int_t i=0; i<entries; i++){
 758       Int_t index = indexes[i];
 759       Int_t run   = treeScan->GetV2()[index];
 760       Float_t chi2=treeScan->GetV1()[index];
 761       Float_t vgg =
 762         printf("%d\t%d\t%f\t%d\t%d\t%d\n",iside,run, chi2,GetVoltage(run,0), GetVoltage(run,1),GetVoltage(run,2));
 763     }
 764   }
 765 }
 766
 767
 768
 769
 770 void MakeAliases(){
 771   //
 772   // use table
 773   //
 774   chain->SetAlias("VggA","GetVoltage(run,0)");
 775   chain->SetAlias("VggC","GetVoltage(run,1)");
 776   chain->SetAlias("VcoA","GetVoltage(run,2)");
 777   chain->SetAlias("VcoC","GetVoltage(run,3)");
 778   chain->SetAlias("VskA","GetVoltage(run,4)");
 779   chain->SetAlias("VskC","GetVoltage(run,5)");
 780   //
 781   // cuts
 782   //
 783   chain->SetAlias("TisOK","mdEdx>5&&entries>400");
 784   chain->SetAlias("CisOK","nCl.fElements>entries*0.5&&eY.fElements<0.01");
 785   chain->SetAlias("ATisOK","RA.mdEdx>5&&RA.entries>400");
 786   chain->SetAlias("ACisOK","RA.nCl.fElements>RA.entries*0.5&&abs(RA.dY.fElements-dY.fElements)<0.3&&RA.eY.fElements<0.01");
 787   //
 788   // shortcuts
 789   //
 790   chain->SetAlias("lX","X.fElements");   //
 791   chain->SetAlias("tY","kY.fElements");   //
 792   chain->SetAlias("dE","sign(Y.fElements)*(tan(10*pi/180.)*X.fElements-abs(Y.fElements))");
 793   chain->SetAlias("RdY","(dY.fElements-RA.dY.fElements)");
 794
 795   chain->SetAlias("dIF","lX-85.2");
 796   chain->SetAlias("dOF","245.8-lX");
 797   chain->SetAlias("dUtY","((GetVoltage(run,0)-GetVoltage(RA.run,0))/400.*tY)");
 798   //
 799   chain->SetAlias("cIF","dUtY/(1+dIF/13.0)");
 800   chain->SetAlias("cOF","dUtY/(1+dOF/22.0)");
 801   chain->SetAlias("cIF2","(dUtY*dIF/13.0)/((1.+dIF/13.0)^2)");
 802   chain->SetAlias("cOF2","(dUtY*dOF/22.0)/((1.+dOF/22.0)^2)");
 803
 804   chain->SetAlias("cE","(GetVoltage(run,0)-GetVoltage(RA.run,0))/400.*sign(dE)/(1+(abs(dE)-1.55)/1.4)");
 805   chain->SetAlias("cE2","(GetVoltage(run,0)-GetVoltage(RA.run,0))/400.*sign(dE)*((abs(dE)-1.55)/1.4)/((1+(abs(dE)-1.55)/1.4)^2)");
 806 }
 807
 808
 809
 810 void MakeAliasesBoth(){
 811   //
 812   // cuts - slect good tracks
 813   //
 814   chain->SetAlias("TisOK","mdEdx>5&&entries>400");
 815   chain->SetAlias("ATisOK","(LTr.fSide==0)*(RA.mdEdx>5&&RA.entries>500)");
 816   chain->SetAlias("CTisOK","(LTr.fSide==1)*(RC.mdEdx>5&&RC.entries>500)");
 817   //
 818   chain->Draw(">>run","TisOK&&(ATisOK||CTisOK)","entryList");
 819   TEntryList *elist = (TEntryList*)gDirectory->Get("run");
 820   chain->SetEntryList(elist);
 821   //
 822   //
 823   chain->SetAlias("CisOK","(nCl.fElements>entries*0.5&&eY.fElements<0.01)");
 824   chain->SetAlias("ACisOK","RA.nCl.fElements>RA.entries*0.5&&abs(RA.dY.fElements-dY.fElements)<0.3&&RA.eY.fElements<0.01");
 825   chain->SetAlias("CCisOK","RC.nCl.fElements>RC.entries*0.5&&abs(RC.dY.fElements-dY.fElements)<0.3&&abs(RC.dZ.fElements-dZ.fElements)<0.3&&RC.eY.fElements<0.01");
 826
 827
 828   //
 829   // voltage table
 830   //
 831   chain->SetAlias("Vgg","((LTr.fSide==0)*GetVoltage(run,0)+(LTr.fSide==1)*GetVoltage(run,1))");
 832   chain->SetAlias("VcoA","((LTr.fSide==0)*GetVoltage(run,2)+(LTr.fSide==1)*GetVoltage(run,3))");
 833   chain->SetAlias("VskA","((LTr.fSide==0)*GetVoltage(run,4)+(LTr.fSide==1)*GetVoltage(run,5))");
 834
 835   chain->SetAlias("dVgg","(((LTr.fSide==0)*(GetVoltage(run,0)-GetVoltage(RA.run,0))+((LTr.fSide==0)*(GetVoltage(run,0)-GetVoltage(RA.run,0))))/400.)");
 836   //
 837   // shortcuts
 838   //
 839   chain->SetAlias("RdY","((LTr.fSide==0)*(dY.fElements-RA.dY.fElements)+(LTr.fSide==1)*(dY.fElements-RC.dY.fElements))");
 840
 841   chain->SetAlias("drift","(1.-abs(LTr.fP[1]+0.)/250)");
 842   chain->SetAlias("lX","X.fElements");   //
 843   chain->SetAlias("tY","kY.fElements");   //
 844   chain->SetAlias("dE","sign(Y.fElements)*(tan(10*pi/180.)*X.fElements-abs(Y.fElements))");
 845
 846   chain->SetAlias("dIF","lX-85.2");
 847   chain->SetAlias("dOF","245.8-lX");
 848   chain->SetAlias("dUtY","dVgg*tY");
 849   //
 850   //
 851   //
 852   chain->SetAlias("cIF","dUtY/(1+dIF/13.0)");
 853   chain->SetAlias("cOF","dUtY/(1+dOF/22.0)");
 854   chain->SetAlias("cIF2","(dUtY*dIF/13.0)/((1.+dIF/13.0)^2)");
 855   chain->SetAlias("cOF2","(dUtY*dOF/22.0)/((1.+dOF/22.0)^2)");
 856
 857   chain->SetAlias("cE","dVgg*sign(dE)/(1+(abs(dE)-1.5)/1.3)");
 858   chain->SetAlias("cE2","dVgg*sign(dE)*((abs(dE)-1.5)/1.3)/((1+(abs(dE)-1.5)/1.3)^2)");
 859
 860 }
 861
 862
 863
 864
 865 void MakeFit(){
 866
 867   Int_t  ntracks=3000000;
 868   TStatToolkit toolkit;
 869   Double_t chi2=0;
 870   Int_t    npoints=0;
 871   TVectorD fitParam;
 872   TMatrixD covMatrix;
 873   TString fstring="";
 874   fstring+="cIF++";
 875   fstring+="cIF2++";
 876   fstring+="cOF++";
 877   fstring+="cOF2++";
 878   fstring+="cE++";
 879   fstring+="cE2++";
 880   TCut cutA="LTr.fBeam>-1&&CisOK&&(ACisOK||CCisOK)";
 881
 882   TString * strA = TStatToolkit::FitPlane(chain,"RdY", fstring.Data(),cutA, chi2,npoints,fitParam,covMatrix,-1.,0, ntracks);
 883   chain->SetAlias("fdY",strA->Data());
 884
 885   printf("sqrt(Chi2/npoints)=%f\n",TMath::Sqrt(chi2/npoints));
 886   chain->Draw("RdY-fdY",cutA);
 887
 888
 889   TF1 fe("fe","[0]/(1+(x-[2])/[1])",2,10);
 890   fe.SetParameters(0.4,1,1.5);
 891
 892
 893 }
 894
 895
 896
 897 //Examples:
 898 // chain->Draw("RdY:(GetVoltage(run,0)-GetVoltage(RA.run,0))*tY/(1+dIF/10.)","TisOK&&ATisOK&&CisOK&&ACisOK&&LTr.fBundle>0&&LTr.fBeam!=3&&dIF<30","",10000)
 899
 900
 901
 902