TPC/CalibMacros/CalibLaserVscan.C

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