TPC/CalibMacros/CalibTimeVdrift.C

   1 /* small macro to generate and update OCDB entries for a given run:
   2
   3 // How to use it...
   4 //
   5 // 1. go to the directory with the "CalibObjectsTrain1.root" file and start your aliroot
   6 //
   7 // 2. copy and paste the following commands
   8 //
   9 gSystem->Load("libSTEER");
  10 gSystem->Load("libANALYSIS");
  11 gSystem->Load("libSTAT");
  12 gSystem->Load("libTPCcalib");
  13 gSystem->AddIncludePath("-I$ALICE_ROOT/STEER");
  14 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
  15 .L $ALICE_ROOT/TPC/CalibMacros/CalibTimeVdrift.C+
  16
  17 //
  18 // 3. setup your OCDB output path e.g:
  19 //
  20 ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
  21 //
  22 // 4. simple execution
  23 //
  24 CalibTimeVdriftGlobal()
  25 //
  26 // 5. try to visualize new entry
  27 //
  28 Int_t run=90000; //
  29 //.x ConfigOCDB.C(run);
  30
  31 ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
  32 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB")
  33 AliCDBManager::Instance()->SetSpecificStorage("TPC/Calib/TimeDrift",ocdbStorage.Data());
  34
  35 AliCDBEntry* entry = AliCDBManager::Instance()->Get("TPC/Calib/TimeDrift",run)
  36 TObjArray * arr = (TObjArray*)entry->GetObject();
  37
  38 TObjArray *picArray = new TObjArray;
  39 MakeDefaultPlots(arr,picArray);
  40
  41 //
  42 //  7. Tricky part - hopefully not neccessary
  43 //     change the default time 0
  44 //     BUG in OCDB - we can not change the SpecificStorage once we read given
  45 //                   entry
  46
  47    .x ConfigOCDB.C(
  48
  49 */
  50
  51 #include "TMap.h"
  52 #include "TGraphErrors.h"
  53 #include "AliExternalTrackParam.h"
  54 #include "TFile.h"
  55 #include "TGraph.h"
  56 #include "TMultiGraph.h"
  57 #include "TCanvas.h"
  58 #include "THnSparse.h"
  59 #include "TLegend.h"
  60 #include "TPad.h"
  61
  62
  63 #include "AliTPCcalibTime.h"
  64 #include "AliSplineFit.h"
  65 #include "AliCDBMetaData.h"
  66 #include "AliCDBId.h"
  67 #include "AliCDBManager.h"
  68 #include "AliCDBStorage.h"
  69 #include "AliTPCcalibBase.h"
  70 #include "AliTPCcalibDB.h"
  71 #include "AliTPCcalibDButil.h"
  72 #include "AliRelAlignerKalman.h"
  73 #include "AliTPCParamSR.h"
  74
  75 TString ocdbStorage="dummy";
  76 //
  77 //
  78 //
  79 void AddHistoGraphs(  TObjArray * vdriftArray, AliTPCcalibTime *timeDrift, Int_t minEntries);
  80 void AddAlignmentGraphs(  TObjArray * vdriftArray, AliTPCcalibTime *timeDrift);
  81 void AddLaserGraphs(  TObjArray * vdriftArray, AliTPCcalibTime *timeDrift);
  82 //
  83
  84 TGraphErrors* FilterGraphMedianAbs(TGraphErrors * graph, Float_t cut,Double_t &medianY);
  85 TGraphErrors* FilterGraphDrift(TGraphErrors * graph, Float_t errSigmaCut, Float_t medianCutAbs);
  86
  87
  88
  89
  90 void CalibTimeVdriftGlobal(Char_t* file="CalibObjectsTrain1.root"){
  91
  92   const Int_t    kMinEntries=500;     // minimal number of entries
  93   //
  94   // 1. Initialization and run range setting
  95   TFile fcalib(file);
  96   AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
  97   Int_t startRun=0, endRun=0;
  98   Int_t startTime=0, endTime=0;
  99   //Int_t startPT=-1, endPT=-1;
 100   //
 101   // find the fist and last run
 102   //
 103   TObjArray *hisArray =timeDrift->GetHistoDrift();
 104   {for (Int_t i=0; i<hisArray->GetEntriesFast(); i++){
 105     THnSparse* addHist=(THnSparse*)hisArray->UncheckedAt(i);
 106     if (addHist->GetEntries()<kMinEntries) continue;
 107     if (!addHist) continue;
 108     TH1D* histo    =addHist->Projection(3);
 109     TH1D* histoTime=addHist->Projection(0);
 110     if (startRun==0){
 111       startRun=histo->FindFirstBinAbove(0);
 112       endRun  =histo->FindLastBinAbove(0);
 113     }else{
 114       startRun=TMath::Min(histo->FindFirstBinAbove(0),startRun);
 115       endRun  =TMath::Max(histo->FindLastBinAbove(0),endRun);
 116     }
 117     if (startTime==0){
 118       startTime=histoTime->FindFirstBinAbove(0);
 119       endTime  =histoTime->FindLastBinAbove(0);
 120     }else{
 121       startTime=TMath::Min(histoTime->FindFirstBinAbove(0),startTime);
 122       endTime  =TMath::Max(histoTime->FindLastBinAbove(0),endTime);
 123     }
 124     delete histo;
 125     delete histoTime;
 126   }}
 127   {for (Int_t i=0; i<hisArray->GetEntriesFast(); i++){
 128       THnSparse* addHist=(THnSparse*)hisArray->At(i);
 129       if (!addHist) continue;
 130       addHist->GetAxis(0)->SetRange(startTime-1,endTime+1);
 131       addHist->GetAxis(3)->SetRange(startRun-1,endRun+1);
 132     }
 133   }
 134   printf("Run range  :\t%d-%d\n", startRun, endRun);
 135   printf("Time range :\t%d-%d\n", startTime, endTime);
 136   //
 137   // 2. extraction of the information
 138   //
 139   TObjArray * vdriftArray = new TObjArray();
 140   AddAlignmentGraphs(vdriftArray,timeDrift);
 141   AddHistoGraphs(vdriftArray,timeDrift,kMinEntries);
 142   AddLaserGraphs(vdriftArray,timeDrift);
 143   //
 144   //
 145   // 3. update of OCDB
 146   //
 147   //
 148   AliCDBMetaData *metaData= new AliCDBMetaData();
 149   metaData->SetObjectClassName("TObjArray");
 150   metaData->SetResponsible("Dag Toppe Larsen");
 151   metaData->SetBeamPeriod(1);
 152   metaData->SetAliRootVersion("05-25-01"); //root version
 153   metaData->SetComment("Calibration of the time dependence of the drift velocity due to pressure and temperature changes");
 154   AliCDBId* id1=NULL;
 155   id1=new AliCDBId("TPC/Calib/TimeDrift", startRun, AliCDBRunRange::Infinity());
 156   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
 157   gStorage->Put(vdriftArray, (*id1), metaData);
 158
 159 }
 160
 161
 162 void UpdateDriftParam(AliTPCParam *param, TObjArray *arr, Int_t startRun){
 163   //
 164   //  update the OCDB entry for the nominal time0
 165   //
 166   //
 167   //  AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
 168   AliTPCParam *paramNew = (AliTPCParam *)param->Clone();
 169   TGraphErrors *grT =  (TGraphErrors *)arr->FindObject("ALIGN_ITSM_TPC_T0");
 170   Double_t deltaTcm = TMath::Median(grT->GetN(),grT->GetY());
 171   Double_t deltaT   = deltaTcm/param->GetDriftV();
 172   paramNew->SetL1Delay(param->GetL1Delay()-deltaT);
 173   paramNew->Update();
 174
 175   AliCDBMetaData *metaData= new AliCDBMetaData();
 176   metaData->SetObjectClassName("TObjArray");
 177   metaData->SetResponsible("Marian Ivanov");
 178   metaData->SetBeamPeriod(1);
 179   metaData->SetAliRootVersion("05-25-02"); //root version
 180   metaData->SetComment("Updated calibration of nominal time 0");
 181   AliCDBId* id1=NULL;
 182   id1=new AliCDBId("TPC/Calib/Parameters", startRun, AliCDBRunRange::Infinity());
 183   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
 184   gStorage->Put(param, (*id1), metaData);
 185
 186 }
 187
 188 void CalibTimeVdrift(Int_t runNumber, AliTPCcalibTime* vdrift, Int_t end=false){
 189   TObjArray* vdriftArray = new TObjArray();
 190   TObjArray * array=vdrift->GetHistoDrift();
 191   if(!array) return;
 192   TIterator* iterator = array->MakeIterator();
 193   iterator->Reset();
 194   THnSparse* hist=NULL;
 195   while((hist=(THnSparseF*)iterator->Next())){
 196     if(!hist) continue;
 197     hist->Print();
 198     if(end) hist->GetAxis(3)->SetRangeUser(runNumber-0.5,       end+0.5);
 199     else    hist->GetAxis(3)->SetRangeUser(runNumber-0.5, runNumber+0.5);
 200     TString name=hist->GetName();
 201     Int_t dim[4]={0,1,2,3};
 202     THnSparse* newHist=hist->Projection(4,dim);
 203     newHist->SetName(name);
 204     vdriftArray->Add(newHist);
 205     TGraphErrors* graph=AliTPCcalibBase::FitSlices(newHist,2,0,400,100,0.05,0.95, kTRUE);
 206     printf("name=%s graph=%i\n", name.Data(), graph==0);
 207     if(!graph || !graph->GetN()) continue;
 208     printf("name=%s graph=%i, N=%i\n", name.Data(), graph==0, graph->GetN());
 209     Int_t pos=name.Index("_");
 210     name=name(pos,name.Capacity()-pos);
 211     TString graphName=graph->ClassName();
 212     graphName+=name;
 213     graphName.ToUpper();
 214     graph->SetName(graphName);
 215     printf("name=%s\n", graphName.Data());
 216     vdriftArray->Add(graph);
 217 //      AliSplineFit* fit=new AliSplineFit();
 218 //      fit->SetGraph(graph);
 219 //      fit->SetMinPoints(graph->GetN()+1);
 220 //      fit->InitKnots(graph,2,0,0.001);
 221 //      fit->SplineFit(0);
 222 //      TString fiName=fit->ClassName();
 223 //      fiName+=type;
 224 //      fiName+=trigger;
 225 //      fiName.ToUpper();
 226 //      fit->SetName(fiName.Data());
 227 //      printf("name=%s\n", fiName.Data());
 228 //      vdriftArray->Add(fit);
 229   }
 230   THnSparse* laserHist=NULL;
 231   TGraphErrors* laserGraph=NULL;
 232   TString laserName="";
 233
 234   laserHist=vdrift->GetHistVdriftLaserA(1);
 235   laserName=laserHist->ClassName();
 236   laserName+="_MEAN_DRIFT_LASER_ALL_A";
 237   laserName.ToUpper();
 238   laserHist->SetName(laserName);
 239   vdriftArray->Add(laserHist);
 240   laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,400,100,0.05,0.95, kTRUE);
 241   if(laserGraph && laserGraph->GetN()){
 242     laserName=laserGraph->GetName();
 243     laserName+="_MEAN_DRIFT_LASER_ALL_A";
 244     laserName.ToUpper();
 245     laserGraph->SetName(laserName);
 246     vdriftArray->Add(laserGraph);
 247   }
 248
 249   laserHist=vdrift->GetHistVdriftLaserC(1);
 250   laserName=laserHist->ClassName();
 251   laserName+="_MEAN_DRIFT_LASER_ALL_C";
 252   laserName.ToUpper();
 253   laserHist->SetName(laserName);
 254   vdriftArray->Add(laserHist);
 255   laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,400,100,0.05,0.95, kTRUE);
 256   if(laserGraph && laserGraph->GetN()){
 257     laserName=laserGraph->GetName();
 258     laserName+="_MEAN_DRIFT_LASER_ALL_C";
 259     laserName.ToUpper();
 260     laserGraph->SetName(laserName);
 261     vdriftArray->Add(laserGraph);
 262   }
 263
 264   AliCDBMetaData *metaData= new AliCDBMetaData();
 265   metaData->SetObjectClassName("TObjArray");
 266   metaData->SetResponsible("Dag Toppe Larsen");
 267   metaData->SetBeamPeriod(1);
 268   metaData->SetAliRootVersion("05-25-01"); //root version
 269   metaData->SetComment("Calibration of the time dependence of the drift velocity due to pressure and temperature changes");
 270   AliCDBId* id1=NULL;
 271   if(end) id1=new AliCDBId("TPC/Calib/TimeDrift", runNumber, end);
 272   else    id1=new AliCDBId("TPC/Calib/TimeDrift", runNumber, runNumber);
 273   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
 274   gStorage->Put(vdriftArray, (*id1), metaData);
 275   printf("done runNumber=%i, end=%i\n", runNumber, end);
 276 }
 277
 278 void CalibrateAll(Char_t* file="CalibObjectsTrain1.root"){
 279   TFile fcalib(file);
 280   AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
 281   THnSparse* addHist=timeDrift->GetHistoDrift("all");
 282   TH1D* histo=addHist->Projection(3);
 283   Int_t start=histo->FindFirstBinAbove(0);
 284   Int_t end  =histo->FindLastBinAbove(0);
 285   CalibTimeVdrift(start-1, timeDrift, end-1);
 286   for(Int_t i=start;i<end;i++){
 287     if(histo->GetBinContent(i)){
 288       printf("start=%i i=%i end=%i i=%i center=%f", start, i, end, i, histo->GetBinCenter(i));
 289       CalibTimeVdrift(i, timeDrift);
 290     }
 291   }
 292 }
 293
 294 void MakePlot(Char_t* file="CalibObjects.root", Char_t* trigger="", /*Int_t useError=true,*/ Int_t separateFits=false, Int_t pointLimit=0, Int_t showHisto=true, Int_t otherColour=false, Int_t spline=true, Double_t yMin=0.0, Double_t yMax=0.0){
 295   TFile fcalib(file);
 296   AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
 297
 298 //Get graphs
 299   TObjArray* grArray=new TObjArray();
 300   TObjArray* hiArray=new TObjArray();
 301   Int_t counter=0;
 302   Int_t pos=counter+1;
 303   Int_t colour=2;
 304   TObjArray * addArray=timeDrift->GetHistoDrift();
 305   if(!addArray) return;
 306   TIterator* iterator = addArray->MakeIterator();
 307   iterator->Reset();
 308   THnSparse* addHist=NULL;
 309
 310   THnSparse* laserHist=NULL;
 311   TString laserName="";
 312
 313   laserHist=timeDrift->GetHistVdriftLaserA(1);
 314   laserName=laserHist->ClassName();
 315   laserName+="_MEAN_VDRIFT_LASER_ALL_A";
 316   laserName.ToUpper();
 317   laserHist->SetName(laserName);
 318   addArray->Add(laserHist);
 319
 320   laserHist=timeDrift->GetHistVdriftLaserC(1);
 321   laserName=laserHist->ClassName();
 322   laserName+="_MEAN_VDRIFT_LASER_ALL_C";
 323   laserName.ToUpper();
 324   laserHist->SetName(laserName);
 325   addArray->Add(laserHist);
 326
 327   while((addHist=(THnSparseF*)iterator->Next())){
 328     if (!addHist) continue;
 329     TString histName=addHist->GetName();
 330     if(!histName.Contains(trigger) && !histName.Contains("COSMICS_ALL")) continue;
 331     addHist->Print();
 332     TGraph* graph=AliTPCcalibBase::FitSlices(addHist,2,0,400,100,0.05,0.95, kTRUE);
 333     TH2D* histo=addHist->Projection(2,0);
 334     if(!graph || !histo || !graph->GetN()) continue;
 335     if(spline){
 336       AliSplineFit fit;
 337       fit.SetGraph(graph);
 338       fit.SetMinPoints(graph->GetN()+1);
 339       fit.InitKnots(graph,2,0,0.001);
 340       fit.SplineFit(0);
 341       graph=fit.MakeGraph(graph->GetX()[0],graph->GetX()[graph->GetN()-1],50000,0);
 342     }
 343     graph->SetNameTitle(addHist->GetName(),"Relative drift velocity change");
 344     graph->GetXaxis()->SetTitle("time");
 345     graph->GetXaxis()->SetTimeFormat();
 346     graph->GetXaxis()->SetTimeDisplay(1);
 347     graph->GetYaxis()->SetTitleOffset(1.2);
 348     graph->GetXaxis()->SetLabelSize(0.035);
 349     graph->GetYaxis()->SetLabelSize(0.04);
 350     graph->GetYaxis()->SetTitle("relative drift velocity change");
 351     if(yMin) graph->SetMinimum(yMin);
 352     if(yMax) graph->SetMaximum(yMax);
 353     printf("bins=%i name=%s couter=%i pos=%i colour=%i Ymax=%f Ymin=%f Xmax=%f Xmin=%f\n", graph->GetXaxis()->GetNbins(), addHist->GetName(), counter, pos, colour, graph->GetYaxis()->GetXmax(), graph->GetYaxis()->GetXmin(), graph->GetXaxis()->GetXmax(), graph->GetXaxis()->GetXmin());
 354     if(graph->GetXaxis()->GetNbins()<=pointLimit) continue;
 355 //    if(!useError&&!spline) for(Int_t j=0; j<graph->GetXaxis()->GetNbins(); j++) graph->SetPointError(j,0,0);
 356 //    if( *((Int_t*)(graph->GetName()))==*((Int_t*)("all")) ){
 357     TString name=graph->GetName();
 358     if(name.Contains("COSMICS_ALL")){
 359       graph->SetLineColor(0+1);
 360       graph->SetMarkerColor(0+1);
 361       grArray->AddAtAndExpand(graph,0);
 362       histo->SetLineColor(1);
 363       histo->SetMarkerColor(1);
 364       hiArray->AddAtAndExpand(histo,0);
 365       pos=counter;
 366     }
 367     else{
 368       graph->SetLineColor(colour);
 369       graph->SetMarkerColor(colour);
 370       grArray->AddAtAndExpand(graph,pos);
 371       histo->SetLineColor(colour);
 372       histo->SetMarkerColor(colour);
 373       hiArray->AddAtAndExpand(histo,pos);
 374       colour++;
 375       if(colour==10) colour++;
 376     }
 377     pos++;
 378     counter++;
 379     graph=0;
 380     addHist=0;
 381   }
 382
 383   if(!separateFits){
 384     TCanvas* plot = new TCanvas("plot", "Relative drift velocity change", 1500, 900);
 385     //Draw graphs
 386     TGraph* grTemp=0;
 387     TH2D* hiTemp=0;
 388     for(Int_t i=0; i<counter; i++){
 389       if(!(*grArray)[i] || !(*hiArray)[i]) continue;
 390       grTemp=(TGraph*)(*grArray)[i];
 391       hiTemp=(TH2D*)(*hiArray)[i];
 392       if(i==0) grTemp->Draw("alp");
 393       else grTemp->Draw("lp");
 394       if(showHisto) hiTemp->Draw("same");
 395     }
 396     //Add legends
 397      TLegend* leg=NULL;
 398      if(pointLimit){
 399        leg = new TLegend(0.4,0.75,1,0.95,NULL,"brNDC");
 400        leg->SetTextSize(0.03);
 401      }
 402      else{
 403        leg = new TLegend(0.4,0.65,1,1,NULL,"brNDC");
 404        leg->SetTextSize(0.02);
 405      }
 406     leg->SetLineColor(0);
 407     leg->SetLineWidth(0);
 408     leg->SetFillColor(0);
 409     TLegendEntry* entry=0;
 410     for(Int_t i=0; i<counter; i++){
 411       if(!(*grArray)[i]) continue;
 412       TString name=(*grArray)[i]->GetName();
 413       Int_t pos=name.Index("_")+1+12;
 414       name=name(pos,name.Capacity()-pos);
 415       entry = leg->AddEntry((*grArray)[i],name,"l");
 416       leg->Draw();
 417     }
 418     plot->SaveAs("relative_drift_velocity_change.pdf");
 419     plot->SaveAs("relative_drift_velocity_change.png");
 420   }
 421   else{
 422     //Draw graphs
 423     TGraph* grTemp=0;
 424     TH2D* hiTemp=0;
 425     for(Int_t i=0; i<counter; i++){
 426       if(!(*grArray)[i] || !(*hiArray)[i]) continue;
 427       grTemp=(TGraph*)(*grArray)[i];
 428       hiTemp=(TH2D*)(*hiArray)[i];
 429       TString name=grTemp->GetName();
 430       Int_t pos=name.Index("_")+1+12;
 431       name=name(pos,name.Capacity()-pos);
 432       grTemp->SetTitle(name);
 433       if(otherColour){
 434         grTemp->SetLineColor(2);
 435         grTemp->SetMarkerColor(2);
 436         hiTemp->SetLineColor(1);
 437         hiTemp->SetMarkerColor(1);
 438       }
 439       TCanvas* plot = new TCanvas(name, name, 1500, 900);
 440       grTemp->Draw("alp");
 441       if(showHisto) hiTemp->Draw("same");
 442       TString name2=name;
 443       name+=".pdf";
 444       name2+=".png";
 445       plot->SaveAs(name);
 446       plot->SaveAs(name2);
 447     }
 448   }
 449 }
 450
 451 void MakeQaPlot(Char_t* file="CalibObjects.root"){
 452   TFile fcalib(file);
 453   AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
 454   TCanvas* plotQa = new TCanvas("plotQa", "Relative drift velocity change QA", 1600, 600);
 455   plotQa->Divide(5,2);
 456   for(Int_t i=0; i<10; i++){
 457     plotQa->cd(i+1);
 458     timeDrift->GetCosmiMatchingHisto(i)->Draw();
 459   }
 460   plotQa->cd(0);
 461   plotQa->SaveAs("relative_drift_velocity_change_qa.pdf");
 462 }
 463
 464 void MakeHistPlot(Char_t* file="CalibObjects.root", Char_t* trigger="COSMICS_ALL"){
 465   TFile fcalib(file);
 466   AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
 467   THnSparse* addHist=timeDrift->GetHistoDrift(trigger);
 468   TH2D* histo=addHist->Projection(2,0);
 469
 470   histo->SetTitle("Relative drift velocity change");
 471   histo->GetXaxis()->SetTitle("time");
 472   histo->GetXaxis()->SetTimeFormat();
 473   histo->GetXaxis()->SetTimeDisplay(1);
 474   histo->GetYaxis()->SetTitleOffset(1.3);
 475   histo->GetYaxis()->SetTitle("relative drift velocity change");
 476   histo->SetLineColor(1);
 477   histo->SetMarkerColor(1);
 478   TGraphErrors* graph=AliTPCcalibBase::FitSlices(addHist,2,0,400,100,0.05,0.95, kTRUE);
 479   if(!graph) return;
 480   graph->SetNameTitle(addHist->GetName(),"Relative drift velocity change");
 481   graph->GetXaxis()->SetTitle("time");
 482   graph->GetXaxis()->SetTimeFormat();
 483   graph->GetXaxis()->SetTimeDisplay(1);
 484   graph->GetYaxis()->SetTitleOffset(1.3);
 485   graph->GetYaxis()->SetTitle("relative drift velocity change");
 486   graph->SetMaximum(0.003);
 487   graph->SetMinimum(-0.003);
 488   graph->SetLineColor(2);
 489   graph->SetMarkerColor(2);
 490   TCanvas* histo_fit = new TCanvas("histo_fit", "Relative drift velocity change histogram/fit", 800, 600);
 491   histo->Draw();
 492   graph->Draw("lp");
 493   histo_fit->SaveAs("relative_drift_velocity_change_histo_fit.pdf");
 494   histo_fit->SaveAs("relative_drift_velocity_change_histo_fit.png");
 495 }
 496
 497
 498
 499 void PrintArray(TObjArray *array){
 500   //
 501   //
 502   //
 503   Int_t entries = array->GetEntries();
 504   for (Int_t i=0; i<entries; i++){
 505     printf("%s\n", array->At(i)->GetName());
 506   }
 507 }
 508
 509
 510 TMultiGraph * MakeMultiGraph(TObjArray *array, const char  *mask, Int_t color, Int_t style){
 511   //
 512   // Make a multi graph with graphs selected
 513   //
 514   if (!array) return 0;
 515   TMultiGraph *mgraph=new TMultiGraph("MultiGraph", "MultiGraph");
 516   TObjArray * maskArray=0;
 517   if (mask){
 518     TString grmask(mask);
 519     maskArray = grmask.Tokenize("*");
 520   }
 521   Double_t ymin=0,ymax=-1;
 522   for (Int_t i=0; i<array->GetEntries();i++){
 523     TGraphErrors *graph= (TGraphErrors*)array->At(i);
 524     if (!graph) continue;
 525     if (maskArray){
 526       Bool_t isOK=kTRUE;
 527       TString str(graph->GetName());
 528       for (Int_t imask=0; imask<maskArray->GetEntries();imask++)
 529         if (!str.Contains(maskArray->At(imask)->GetName())==0) isOK=kFALSE;
 530       if (!isOK) continue;
 531     }
 532     if (ymax<ymin){
 533       ymin=graph->GetMinimum();
 534       ymax=graph->GetMaximum();
 535     }
 536     ymin = TMath::Min(ymin,graph->GetMinimum());
 537     ymax = TMath::Max(ymax,graph->GetMaximum());
 538   }
 539   ymin=ymin-(ymax-ymin)*0.2;
 540   ymax=ymax+(ymax-ymin)*0.5;
 541   Bool_t first=kTRUE;
 542   for (Int_t i=0; i<array->GetEntries();i++){
 543     TGraphErrors *graph= (TGraphErrors*)array->At(i);
 544     if (!graph) continue;
 545     if (maskArray){
 546       Bool_t isOK=kTRUE;
 547       TString str(graph->GetName());
 548       for (Int_t imask=0; imask<maskArray->GetEntries();imask++)
 549         if (str.Contains(maskArray->At(imask)->GetName())==0) isOK=kFALSE;
 550       if (!isOK) continue;
 551     }
 552     graph->SetMarkerColor(color);
 553     graph->SetMarkerStyle(style);
 554     graph->SetLineColor(color);
 555     graph->SetMaximum(ymax);
 556     graph->SetMinimum(ymin);
 557     if (first)  {first=kFALSE;}
 558     mgraph->Add(graph);
 559   }
 560   return mgraph;
 561 }
 562
 563
 564 TGraphErrors* FilterGraphDrift(TGraphErrors * graph, Float_t errSigmaCut, Float_t medianCutAbs){
 565   // 2 filters:
 566   //    1. filter graph - error cut errSigmaCut
 567   //    2. filter graph - medianCutAbs around median
 568   //
 569   // errSigmaCut   - cut on error
 570   // medianCutAbs  - cut on value around median
 571   Double_t dummy=0;               //
 572   //
 573   // 1. filter graph - error cut errSigmaCut
 574   //
 575   TGraphErrors *graphF;
 576   graphF = AliTPCcalibDButil::FilterGraphMedianErr(graph,errSigmaCut,dummy);
 577   delete graph;
 578   if (!graphF) return 0;
 579   graph = AliTPCcalibDButil::FilterGraphMedianErr(graphF,errSigmaCut,dummy);
 580   delete graphF;
 581   if (!graph) return 0;
 582   //
 583   // filter graph - kMedianCutAbs around median
 584   //
 585   graphF=FilterGraphMedianAbs(graph, medianCutAbs,dummy);
 586   delete graph;
 587   if (!graphF) return 0;
 588   graph=FilterGraphMedianAbs(graphF, medianCutAbs,dummy);
 589   delete graphF;
 590   if (!graph) return 0;
 591   return graph;
 592 }
 593
 594
 595
 596 TGraphErrors* FilterGraphMedianAbs(TGraphErrors * graph, Float_t cut,Double_t &medianY){
 597   //
 598   // filter outlyer measurement
 599   // Only points around median +- cut filtered
 600   //
 601   if (!graph) return  0;
 602   Int_t kMinPoints=2;
 603   Int_t npoints0 = graph->GetN();
 604   Int_t npoints=0;
 605   Float_t  rmsY=0;
 606   Double_t *outx=new Double_t[npoints0];
 607   Double_t *outy=new Double_t[npoints0];
 608   Double_t *errx=new Double_t[npoints0];
 609   Double_t *erry=new Double_t[npoints0];
 610   //
 611   //
 612   if (npoints0<kMinPoints) return 0;
 613   for (Int_t iter=0; iter<3; iter++){
 614     npoints=0;
 615     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
 616       if (graph->GetY()[ipoint]==0) continue;
 617       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
 618       outx[npoints]  = graph->GetX()[ipoint];
 619       outy[npoints]  = graph->GetY()[ipoint];
 620       errx[npoints]  = graph->GetErrorX(ipoint);
 621       erry[npoints]  = graph->GetErrorY(ipoint);
 622       npoints++;
 623     }
 624     if (npoints<=1) break;
 625     medianY  =TMath::Median(npoints,outy);
 626     rmsY   =TMath::RMS(npoints,outy);
 627   }
 628   TGraphErrors *graphOut=0;
 629   if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
 630   return graphOut;
 631 }
 632
 633
 634 void AddHistoGraphs(  TObjArray * vdriftArray, AliTPCcalibTime *timeDrift, Int_t minEntries){
 635   //
 636   // Add graphs corresponding to the alignment
 637   //
 638   const Double_t kErrSigmaCut=5;      // error sigma cut - for filtering
 639   const Double_t kMedianCutAbs=0.03;  // error sigma cut - for filtering
 640   //
 641   TObjArray * array=timeDrift->GetHistoDrift();
 642   if (array){
 643     THnSparse* hist=NULL;
 644     // 2.a) cosmics with different triggers
 645     for (Int_t i=0; i<array->GetEntriesFast();i++){
 646       hist=(THnSparseF*)array->UncheckedAt(i);
 647       if(!hist) continue;
 648       if (hist->GetEntries()<minEntries) continue;
 649       //hist->Print();
 650       TString name=hist->GetName();
 651       Int_t dim[4]={0,1,2,3};
 652       THnSparse* newHist=hist->Projection(4,dim);
 653       newHist->SetName(name);
 654       TGraphErrors* graph=AliTPCcalibBase::FitSlices(newHist,2,0,400,100,0.05,0.95, kTRUE);
 655       printf("name=%s graph=%i, N=%i\n", name.Data(), graph==0, graph->GetN());
 656       Int_t pos=name.Index("_");
 657       name=name(pos,name.Capacity()-pos);
 658       TString graphName=graph->ClassName();
 659       graphName+=name;
 660       graphName.ToUpper();
 661       //
 662       graph = FilterGraphDrift(graph, kErrSigmaCut, kMedianCutAbs);
 663       if (!graph) {
 664         printf("Graph =%s filtered out\n", name.Data());
 665         continue;
 666       }
 667       //
 668       graph->SetMarkerStyle(i%8+20);
 669       graph->SetMarkerColor(i%7);
 670       graph->GetXaxis()->SetTitle("Time");
 671       graph->GetYaxis()->SetTitle("v_{dcor}");
 672       graph->SetName(graphName);
 673       graph->SetTitle(graphName);
 674       printf("Graph %d\t=\t%s\n", i, graphName.Data());
 675       vdriftArray->Add(graph);
 676     }
 677   }
 678 }
 679
 680
 681
 682
 683 void AddAlignmentGraphs(  TObjArray * vdriftArray, AliTPCcalibTime *timeDrift){
 684   //
 685   // Add graphs corresponding to alignment to the object array
 686   //
 687   TObjArray *arrayITS=0;
 688   TObjArray *arrayTOF=0;
 689   TObjArray *arrayTRD=0;
 690   TMatrixD *mstatITS=0;
 691   TMatrixD *mstatTOF=0;
 692   TMatrixD *mstatTRD=0;
 693   //
 694   arrayITS=timeDrift->GetAlignITSTPC();
 695   arrayTRD=timeDrift->GetAlignTRDTPC();
 696   arrayTOF=timeDrift->GetAlignTOFTPC();
 697   mstatITS= AliTPCcalibDButil::MakeStatRelKalman(arrayITS,0.9,50,0.025);
 698   mstatTOF= AliTPCcalibDButil::MakeStatRelKalman(arrayTOF,0.9,1000,0.025);
 699   mstatTRD= AliTPCcalibDButil::MakeStatRelKalman(arrayTRD,0.9,50,0.025);
 700   //
 701   TObjArray * arrayITSP= AliTPCcalibDButil::SmoothRelKalman(arrayITS,*mstatITS, 0, 5.);
 702   TObjArray * arrayITSM= AliTPCcalibDButil::SmoothRelKalman(arrayITS,*mstatITS, 1, 5.);
 703   TObjArray * arrayITSB= AliTPCcalibDButil::SmoothRelKalman(arrayITSP,arrayITSM);
 704   TObjArray * arrayTOFP= AliTPCcalibDButil::SmoothRelKalman(arrayTOF,*mstatTOF, 0, 5.);
 705   TObjArray * arrayTOFM= AliTPCcalibDButil::SmoothRelKalman(arrayTOF,*mstatTOF, 1, 5.);
 706   TObjArray * arrayTOFB= AliTPCcalibDButil::SmoothRelKalman(arrayTOFP,arrayTOFM);
 707
 708   TObjArray * arrayTRDP= AliTPCcalibDButil::SmoothRelKalman(arrayTRD,*mstatTRD, 0, 5.);
 709   TObjArray * arrayTRDM= AliTPCcalibDButil::SmoothRelKalman(arrayTRD,*mstatTRD, 1, 5.);
 710   TObjArray * arrayTRDB= AliTPCcalibDButil::SmoothRelKalman(arrayTRDP,arrayTRDM);
 711
 712   //
 713   //
 714   Int_t entries=TMath::Max(arrayITS->GetEntriesFast(),arrayTOF->GetEntriesFast());
 715   TObjArray *arrays[12]={arrayITS, arrayITSP, arrayITSM, arrayITSB,
 716                          arrayTRD, arrayTRDP, arrayTRDM, arrayTRDB,
 717                          arrayTOF, arrayTOFP, arrayTOFM, arrayTOFB};
 718   TString   grnames[12]={"ALIGN_ITS", "ALIGN_ITSP", "ALIGN_ITSM", "ALIGN_ITSB",
 719                          "ALIGN_TRD", "ALIGN_TRDP", "ALIGN_TRDM","ALIGN_TRDB",
 720                          "ALIGN_TOF", "ALIGN_TOFP", "ALIGN_TOFM","ALIGN_TOFB"};
 721   TString   grpar[9]={"DELTAPSI", "DELTATHETA", "DELTAPHI",
 722                       "DELTAX", "DELTAY", "DELTAZ",
 723                       "DRIFTVD", "T0", "VDGY"};
 724
 725
 726   TVectorD vX(entries);
 727   TVectorD vY(entries);
 728   TVectorD vEx(entries);
 729   TVectorD vEy(entries);
 730   TObjArray *arr=0;
 731   for (Int_t iarray=0; iarray<12; iarray++){
 732     arr = arrays[iarray];
 733     if (arr==0) continue;
 734     for (Int_t ipar=0; ipar<9; ipar++){
 735       Int_t counter=0;
 736       for (Int_t itime=0; itime<arr->GetEntriesFast(); itime++){
 737         AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) arr->UncheckedAt(itime);
 738         if (!kalman) continue;
 739         vX[counter]=kalman->GetTimeStamp();
 740         vY[counter]=(*(kalman->GetState()))[ipar];
 741         if (ipar==6) vY[counter]=1./(*(kalman->GetState()))[ipar]-1;
 742         vEx[counter]=0;
 743         vEy[counter]=TMath::Sqrt((*(kalman->GetStateCov()))(ipar,ipar));
 744         counter++;
 745       }
 746
 747       TGraphErrors * graph=new TGraphErrors(counter, vX.GetMatrixArray(),
 748                                           vY.GetMatrixArray(),
 749                                           vEx.GetMatrixArray(),
 750                                           vEy.GetMatrixArray());
 751       TString grName=grnames[iarray];
 752       grName+="_TPC_";
 753       grName+=grpar[ipar];
 754       graph->SetName(grName.Data());
 755       vdriftArray->AddLast(graph);
 756     }
 757   }
 758 }
 759
 760
 761
 762 void AddLaserGraphs(  TObjArray * vdriftArray, AliTPCcalibTime *timeDrift){
 763   //
 764   // Add laser graphs
 765   //
 766   THnSparse* laserHist=NULL;
 767   TGraphErrors* laserGraph=NULL;
 768   TString laserName="";
 769   laserHist=timeDrift->GetHistVdriftLaserA(1);
 770   laserName=laserHist->ClassName();
 771   laserName+="_MEAN_DRIFT_LASER_ALL_A";
 772   laserName.ToUpper();
 773   laserHist->SetName(laserName);
 774   laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,300,10,0.05,0.95, kTRUE);
 775   laserGraph=FilterGraphDrift(laserGraph, 5,0.04);
 776   //
 777   if(laserGraph && laserGraph->GetN()){
 778     laserName=laserGraph->GetName();
 779     laserName+="_MEAN_DRIFT_LASER_ALL_A";
 780     laserName.ToUpper();
 781     laserGraph->SetName(laserName);
 782     vdriftArray->Add(laserGraph);
 783   }
 784
 785   laserHist=NULL;
 786   laserGraph=NULL;
 787   laserName="";
 788
 789   laserHist=timeDrift->GetHistVdriftLaserC(1);
 790   laserName=laserHist->ClassName();
 791   laserName+="_MEAN_DRIFT_LASER_ALL_C";
 792   laserName.ToUpper();
 793   laserHist->SetName(laserName);
 794   //vdriftArray->Add(laserHist);// really add the whole histogram ??
 795   laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,300,10,0.05,0.95, kTRUE);
 796   laserGraph = FilterGraphDrift(laserGraph, 5,0.04);
 797   if(laserGraph && laserGraph->GetN()){
 798     laserName=laserGraph->GetName();
 799     laserName+="_MEAN_DRIFT_LASER_ALL_C";
 800     laserName.ToUpper();
 801     laserGraph->SetName(laserName);
 802     vdriftArray->Add(laserGraph);
 803   }
 804
 805   // 2.c) additional informaion from laser for A and C side : DELAY
 806   laserHist=NULL;
 807   laserGraph=NULL;
 808   laserName="";
 809
 810   laserHist=timeDrift->GetHistVdriftLaserA(0);
 811   laserName=laserHist->ClassName();
 812   laserName+="_MEAN_DELAY_LASER_ALL_A";
 813   laserName.ToUpper();
 814   laserHist->SetName(laserName);
 815   laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,100,10,0.05,0.95, kTRUE);
 816
 817   if(laserGraph && laserGraph->GetN()){
 818     laserName=laserGraph->GetName();
 819     laserName+="_MEAN_DELAY_LASER_ALL_A";
 820     laserName.ToUpper();
 821     laserGraph->SetName(laserName);
 822     vdriftArray->Add(laserGraph);
 823   }
 824
 825   laserHist=NULL;
 826   laserGraph=NULL;
 827   laserName="";
 828
 829   laserHist=timeDrift->GetHistVdriftLaserC(0);
 830   laserName=laserHist->ClassName();
 831   laserName+="_MEAN_DELAY_LASER_ALL_C";
 832   laserName.ToUpper();
 833   laserHist->SetName(laserName);
 834   //vdriftArray->Add(laserHist);// really add the whole histogram ??
 835   laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,100,10,0.05,0.95, kTRUE);
 836   if(laserGraph && laserGraph->GetN()){
 837     laserName=laserGraph->GetName();
 838     laserName+="_MEAN_DELAY_LASER_ALL_C";
 839     laserName.ToUpper();
 840     laserGraph->SetName(laserName);
 841     vdriftArray->Add(laserGraph);
 842   }
 843
 844   // 2.d) additional informaion from laser for A and C side : GlobalY Gradient
 845   laserHist=NULL;
 846   laserGraph=NULL;
 847   laserName="";
 848
 849   laserHist=timeDrift->GetHistVdriftLaserA(2);
 850   laserName=laserHist->ClassName();
 851   laserName+="_MEAN_GLOBALYGRADIENT_LASER_ALL_A";
 852   laserName.ToUpper();
 853   laserHist->SetName(laserName);
 854   //vdriftArray->Add(laserHist);  // really add the whole histogram ??
 855   laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,100,10,0.05,0.95, kTRUE);
 856   laserGraph=FilterGraphDrift(laserGraph, 5,0.04);
 857   if(laserGraph && laserGraph->GetN()){
 858     laserName=laserGraph->GetName();
 859     laserName+="_MEAN_GLOBALYGRADIENT_LASER_ALL_A";
 860     laserName.ToUpper();
 861     laserGraph->SetName(laserName);
 862     vdriftArray->Add(laserGraph);
 863   }
 864
 865   laserHist=NULL;
 866   laserGraph=NULL;
 867   laserName="";
 868
 869   laserHist=timeDrift->GetHistVdriftLaserC(2);
 870   laserName=laserHist->ClassName();
 871   laserName+="_MEAN_GLOBALYGRADIENT_LASER_ALL_C";
 872   laserName.ToUpper();
 873   laserHist->SetName(laserName);
 874   laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,100,10,0.05,0.95, kTRUE);
 875   laserGraph=FilterGraphDrift(laserGraph, 5,0.04);
 876   if(laserGraph && laserGraph->GetN()){
 877     laserName=laserGraph->GetName();
 878     laserName+="_MEAN_GLOBALYGRADIENT_LASER_ALL_C";
 879     laserName.ToUpper();
 880     laserGraph->SetName(laserName);
 881     vdriftArray->Add(laserGraph);
 882   }
 883 }
 884
 885 void SetDefaultGraphDrift(TGraph *graph, Int_t color, Int_t style){
 886   //
 887   //
 888   //
 889   graph->GetXaxis()->SetTimeDisplay(kTRUE);
 890   graph->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
 891   graph->SetMaximum( 0.025);
 892   graph->SetMinimum(-0.025);
 893   graph->GetXaxis()->SetTitle("Time");
 894   graph->GetYaxis()->SetTitle("v_{dcorr}");
 895   //
 896   graph->GetYaxis()->SetLabelSize(0.03);
 897   graph->GetXaxis()->SetLabelSize(0.03);
 898   //
 899   graph->GetXaxis()->SetNdivisions(10,5,0);
 900   graph->GetYaxis()->SetNdivisions(10,5,0);
 901   //
 902   graph->GetXaxis()->SetLabelOffset(0.02);
 903   graph->GetYaxis()->SetLabelOffset(0.005);
 904   //
 905   graph->GetXaxis()->SetTitleOffset(1.3);
 906   graph->GetYaxis()->SetTitleOffset(1.2);
 907   //
 908   graph->SetMarkerColor(color);
 909   graph->SetLineColor(color);
 910   graph->SetMarkerStyle(style);
 911 }
 912
 913 void SetPadStyle(TPad *pad, Float_t mx0, Float_t mx1, Float_t my0, Float_t my1){
 914   //
 915   //
 916   pad->SetTicks(1,1);
 917   pad->SetMargin(mx0,mx1,my0,my1);
 918 }
 919
 920
 921 void MakeDefaultPlots(TObjArray * arr, TObjArray *picArray){
 922   //
 923   //
 924   //
 925   // margins
 926   Float_t mx0=0.12, mx1=0.1, my0=0.15, my1=0.1;
 927   //
 928   TGraphErrors* laserA       =(TGraphErrors*)arr->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
 929   TGraphErrors* laserC       =(TGraphErrors*)arr->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
 930   TGraphErrors* cosmic       =(TGraphErrors*)arr->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
 931   TGraphErrors* cross        =(TGraphErrors*)arr->FindObject("TGRAPHERRORS_VDRIFT_CROSS_ALL");
 932   TGraphErrors* itstpcP       =(TGraphErrors*)arr->FindObject("ALIGN_ITSP_TPC_DRIFTVD");
 933   TGraphErrors* itstpcM       =(TGraphErrors*)arr->FindObject("ALIGN_ITSM_TPC_DRIFTVD");
 934   //
 935   if (laserA)  SetDefaultGraphDrift(laserA,2,25);
 936   if (laserC)  SetDefaultGraphDrift(laserC,4,26);
 937   if (cosmic)  SetDefaultGraphDrift(cosmic,3,27);
 938   if (cross)   SetDefaultGraphDrift(cross,4,28);
 939   if (itstpcP) SetDefaultGraphDrift(itstpcP,2,29);
 940   if (itstpcM) SetDefaultGraphDrift(itstpcM,4,30);
 941   //
 942   //
 943   TPad *pad=0;
 944   //
 945   // Laser-Laser
 946   //
 947   if (laserA&&laserC){
 948     pad = new TCanvas("TPCLaserVDrift","TPCLaserVDrift");
 949     laserA->Draw("alp");
 950     SetPadStyle(pad,mx0,mx1,my0,my1);
 951     laserA->Draw("apl");
 952     laserC->Draw("p");
 953     TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
 954     legend->AddEntry(laserA,"Laser A side");
 955     legend->AddEntry(laserC,"Laser C side");
 956     legend->Draw();
 957     picArray->AddLast(pad);
 958   }
 959
 960   if (itstpcP&&itstpcM){
 961     pad = new TCanvas("ITSTPC","ITSTPC");
 962     itstpcP->Draw("alp");
 963     SetPadStyle(pad,mx0,mx1,my0,my1);
 964     itstpcP->Draw("alp");
 965     gPad->Clear();
 966     itstpcM->Draw("apl");
 967     itstpcP->Draw("p");
 968     TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
 969     legend->AddEntry(itstpcP,"ITS-TPC smooth forward");
 970     legend->AddEntry(itstpcM,"ITS-TPC smooth back");
 971     legend->Draw();
 972     picArray->AddLast(pad);
 973   }
 974
 975   if (itstpcP&&laserA){
 976     pad = new TCanvas("ITSTPC_LASER","ITSTPC_LASER");
 977     SetPadStyle(pad,mx0,mx1,my0,my1);
 978     laserA->Draw("alp");
 979     itstpcM->Draw("p");
 980     TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
 981     legend->AddEntry(laserA,"TPC laser");
 982     legend->AddEntry(itstpcM,"ITS-TPC smooth back");
 983     legend->Draw();
 984     picArray->AddLast(pad);
 985   }
 986
 987   if (itstpcP&&cross){
 988     pad = new TCanvas("ITSTPC_CROSS","ITSTPC_CROSS");
 989     SetPadStyle(pad,mx0,mx1,my0,my1);
 990     itstpcP->Draw("alp");
 991     pad->Clear();
 992     cross->Draw("ap");
 993     itstpcP->Draw("p");
 994     //
 995     TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
 996
 997     legend->AddEntry(cross,"TPC cross tracks");
 998     legend->AddEntry(itstpcP,"ITS-TPC smooth back");
 999     legend->Draw();
1000     picArray->AddLast(pad);
1001   }
1002   if (itstpcP&&cosmic){
1003     pad = new TCanvas("ITSTPC_COSMIC","ITSTPC_COSMIC");
1004     SetPadStyle(pad,mx0,mx1,my0,my1);
1005     itstpcP->Draw("alp");
1006     pad->Clear();
1007     cosmic->Draw("ap");
1008     itstpcP->Draw("p");
1009     //
1010     TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
1011
1012     legend->AddEntry(cosmic,"TPC cross tracks0 up-down");
1013     legend->AddEntry(itstpcP,"ITS-TPC smooth back");
1014     legend->Draw();
1015     picArray->AddLast(pad);
1016   }
1017 }
1018
1019
1020 TObjArray * SelectEntries(TObjArray* array, char * mask){
1021   //
1022   //
1023   //
1024   TObjArray * arraySelect = new TObjArray;
1025   TObjArray * maskArray=0;
1026   if (mask){
1027     TString grmask(mask);
1028     maskArray = grmask.Tokenize("*");
1029   }
1030   for (Int_t i=0; i<array->GetEntries();i++){
1031     TObject *graph= array->At(i);
1032     if (!graph) continue;
1033     if (maskArray){
1034       Bool_t isOK=kTRUE;
1035       TString str(graph->GetName());
1036       for (Int_t imask=0; imask<maskArray->GetEntries();imask++)
1037         if (str.Contains(maskArray->At(imask)->GetName())==0) isOK=kFALSE;
1038       if (!isOK) continue;
1039       if (isOK) {
1040         printf("%s\n",graph->GetName());
1041         arraySelect->AddLast(graph->Clone());
1042       }
1043     }
1044   }
1045   return arraySelect;
1046 }
1047
1048
1049 void Substract(TGraphErrors *refgr, TObjArray *arraySelect){
1050   //
1051   //
1052   /*
1053   TObjArray *arraySelect = SelectEntries(arr,"COSMIC");
1054   TGraphErrors *refgr= (TGraphErrors*)arr->FindObject("ALIGN_ITSP_TPC_DRIFTVD");
1055   Substract(refgr,arraySelect);
1056   */
1057   Float_t mx0=0.12, mx1=0.1, my0=0.15, my1=0.5;
1058   {
1059   for (Int_t igr=0; igr<arraySelect->GetEntries(); igr++){
1060     TGraphErrors *gr = (TGraphErrors *)arraySelect->At(igr);
1061     Int_t ngr = gr->GetN();
1062      for (Int_t i=0; i<ngr;i++){
1063       gr->GetY()[i]-=refgr->Eval(gr->GetX()[i]);
1064       gr->GetY()[i]*=250/0.264;
1065     }
1066     gr->GetXaxis()->SetRangeUser(refgr->GetXaxis()->GetXmin(),refgr->GetXaxis()->GetXmax());
1067     SetDefaultGraphDrift(gr,((igr)),20+igr);
1068     gr->SetMaximum(8.);
1069     gr->SetMinimum(4.);
1070     gr->GetYaxis()->SetTitle("#Delta_{t}(Time Bin)");
1071   }
1072   //
1073   TPad * pad = new TCanvas("Delays","Delays",1000,800);
1074   SetPadStyle(pad,mx0,mx1,my0,my1);
1075   TLegend *legend = new TLegend(mx0+0.01,1-my1, 1-mx1, 1-0.01, "Time Offset");
1076   {
1077     arraySelect->At(0)->Draw("ap");
1078     for (Int_t igr=0; igr<arraySelect->GetEntries(); igr++){
1079       TGraphErrors *gr = (TGraphErrors *)arraySelect->At(igr);
1080       if (gr->GetN()<5) continue;
1081       if (gr) gr->Draw("p");
1082       legend->AddEntry(gr,gr->GetName());
1083     }
1084   }
1085   legend->Draw();
1086   }
1087 }