PWGPP/QA/Tracking/ExpertQA/AnalyzedEdx.C

   1 /*
   2   Int_t year=2012;
   3   Int_t period=5; // 0=a, 1=b, 2=c
   4   .x $HOME/rootlogon.C
   5   .L ~/alice/reconstruction/trunk/QAplots/source/highPt/AnalyzedEdx.C+
   6   //SelectV0sforPID();
   7   MakeSelectedArrays(4,4000);
   8   FitSelected(10000,100000,year,period);
   9
  10  */
  11
  12 /*
  13   To implement:
  14   1. Selection of runs:
  15      a.)  Redo run wise gain calibration. Important for the bad periods
  16      b.)  Remove runs with outlier peak position.
  17      c.)  Make sector vise calibration. (To eliminate trigger bias)
  18
  19
  20
  21 */
  22
  23 #include "TFile.h"
  24 #include "TTree.h"
  25 #include "TVectorD.h"
  26 #include "TMatrixD.h"
  27 #include "TH2.h"
  28 #include "TF1.h"
  29 #include "TTreeStream.h"
  30 #include "AliMathBase.h"
  31 #include "TSystem.h"
  32 #include "TChain.h"
  33 #include "TDatabasePDG.h"
  34 #include "TRandom.h"
  35 #include "AliTPCcalibBase.h"
  36 //
  37 #include "AliESDv0.h"
  38 #include "AliESDtrack.h"
  39 #include "TMath.h"
  40 #include "AliXRDPROOFtoolkit.h"
  41
  42
  43 TTree * tree  = 0;
  44 Int_t run=0;
  45 Int_t period=0;
  46 TTreeSRedirector *pcstream = 0; //new TTreeSRedirector("trend.root");
  47 TObjArray * arrayPions=0;
  48 TObjArray * arrayProtons=0;
  49 TObjArray * arrayElectrons=0;
  50 TMatrixD *ppullPionsV0=0;
  51 TMatrixD *ppullProtonsV0=0;
  52 TMatrixD *ppullElectronsV0=0;
  53
  54
  55
  56 //void SetPionAliases();
  57 //void FitSlopeROC();
  58 //void AnalyzedEdxFile(const char * fname="");
  59
  60
  61
  62 void SelectV0sforPID( const char * finput="highptAll.list"){
  63   //
  64   // Code to select identified V0 for the PID
  65   // As an input chain of filter trees is used
  66   // Parameter:
  67   //   finput - name of the list file
  68   // Oputput:
  69   //   file - V0Selected.root
  70   TChain * chain  = AliXRDPROOFtoolkit::MakeChainRandom(finput,"V0s",0,1000);
  71   chain->SetCacheSize(1000000000);
  72   //
  73   TDatabasePDG pdg;
  74   Double_t massLambda = pdg.GetParticle("Lambda0")->Mass();
  75   Double_t massK0 = pdg.GetParticle("K0")->Mass();
  76   Double_t massPion = pdg.GetParticle("pi+")->Mass();
  77   Double_t massProton = pdg.GetParticle("proton")->Mass();
  78   //
  79   //
  80   chain->SetAlias("massPion",Form("(%f+0)",massPion));
  81   chain->SetAlias("massProton",Form("(%f+0)",massProton));
  82   chain->SetAlias("massK0",Form("(%f+0)",massK0));
  83   chain->SetAlias("massLambda",Form("(%f+0)",massLambda));
  84   // delta of mass
  85   chain->SetAlias("K0Delta","(v0.GetEffMass(2,2)-massK0)");
  86   chain->SetAlias("LDelta","(v0.GetEffMass(4,2)-massLambda)");
  87   chain->SetAlias("ALDelta","(v0.GetEffMass(2,4)-massLambda)");
  88   chain->SetAlias("EDelta","(v0.GetEffMass(0,0))");
  89   // pull of the mass
  90   chain->SetAlias("K0Pull","(v0.GetEffMass(2,2)-massK0)/v0.GetKFInfo(2,2,1)");
  91   chain->SetAlias("LPull","(v0.GetEffMass(4,2)-massLambda)/v0.GetKFInfo(4,2,1)");
  92   chain->SetAlias("ALPull","(v0.GetEffMass(2,4)-massLambda)/v0.GetKFInfo(2,4,1)");
  93   chain->SetAlias("EPull","EDelta/v0.GetKFInfo(0,0,1)");
  94   // effective pull of the mass - (empirical values form fits)
  95   chain->SetAlias("K0PullEff","K0Delta/sqrt((3.63321e-03)**2+(5.68795e-04*v0.Pt())**2)");
  96   chain->SetAlias("LPullEff","LDelta/sqrt((1.5e-03)**2+(1.8e-04*v0.Pt())**2)");
  97   chain->SetAlias("ALPullEff","ALDelta/sqrt((1.5e-03)**2+(1.8e-04*v0.Pt())**2)");
  98   chain->SetAlias("EPullEff","v0.GetEffMass(0,0)/sqrt((5e-03)**2+(1.e-04*v0.Pt())**2)");
  99   //
 100   //
 101   chain->SetAlias("dEdx0DProton","AliMathBase::BetheBlochAleph(track0.fIp.P()/massProton)");
 102   chain->SetAlias("dEdx1DProton","AliMathBase::BetheBlochAleph(track1.fIp.P()/massProton)");
 103   chain->SetAlias("dEdx0DPion","AliMathBase::BetheBlochAleph(track0.fIp.P()/massPion)");
 104   chain->SetAlias("dEdx1DPion","AliMathBase::BetheBlochAleph(track1.fIp.P()/massPion)");
 105   //
 106   // V0 - cuts -PID,
 107   //
 108   chain->SetAlias("cutDist","sqrt((track0.fIp.fP[0]-track1.fIp.fP[0])**2+(track0.fIp.fP[1]-track1.fIp.fP[1])**2)>3");
 109   chain->SetAlias("cutLong","track0.GetTPCClusterInfo(3,1,0)-5*abs(track0.fP[4])>130&&track1.GetTPCClusterInfo(3,1,0)>130-5*abs(track0.fP[4])");
 110   chain->SetAlias("cutPID","track0.fTPCsignal>0&&track1.fTPCsignal>0");
 111   chain->SetAlias("cutResol","sqrt(track0.fC[14]/track0.fP[4])<0.15&&sqrt(track1.fC[14]/track1.fP[4])<0.15");
 112   chain->SetAlias("cutV0","cutPID&&cutDist&&cutLong&&cutResol");
 113   //
 114   //
 115   chain->SetAlias("K0Selected",      "abs(K0Pull)<3. &&abs(K0PullEff)<3.  && abs(LPull)>3  && abs(ALPull)>3  &&v0.PtArmV0()>0.11");
 116   chain->SetAlias("LambdaSelected",  "abs(LPull)<3.  &&abs(LPullEff)<3.   && abs(K0Pull)>3 && abs(EPull)>3 && abs(EDelta)>0.05");
 117   chain->SetAlias("ALambdaSelected", "abs(ALPull)<3. &&abs(ALPullEff)<3   && abs(K0Pull)>3 && abs(EPull)>3 &&abs(EDelta)>0.05");
 118   //
 119   chain->SetAlias("GammaSelected", "abs(EPull)<3     && abs(K0Pull)>3 && abs(LPull)>3 && abs(ALPull)>3");
 120   //
 121   // Gamma PID selection
 122   //
 123
 124
 125   //
 126   //
 127   TFile *fselected = TFile::Open("V0Selected.root","recreate");
 128   TTree * treeK0     =   chain->CopyTree("type==8&&cutV0&&K0Selected");
 129   TTree * treeLambda =   chain->CopyTree("type==4&&cutV0&&LambdaSelected");
 130   TTree * treeALambda =   chain->CopyTree("type==2&&cutV0&&ALambdaSelected");
 131   TTree * treeGamma =   chain->CopyTree("type==1&&cutV0&&GammaSelected");
 132   //
 133   TTree * trees[4]={treeK0,treeLambda, treeGamma,treeALambda};
 134   TList * aliases = chain->GetListOfAliases();
 135   Int_t nalias= aliases->GetEntries();
 136
 137   for (Int_t i=0; i<4; i++){
 138     for (Int_t ialias=0; ialias<nalias; ialias++){
 139       TNamed *alias = (TNamed*)aliases->At(ialias);
 140       trees[i]->SetAlias(alias->GetName(),alias->GetTitle());
 141     }
 142   }
 143   treeK0->Write("treeK0");
 144   treeLambda->Write("treeLambda");
 145   treeALambda->Write("treeALambda");
 146   treeGamma->Write("treeGamma");
 147   fselected->Close();
 148   //
 149 }
 150
 151
 152
 153 Double_t GetPullMass(AliESDv0 * v0, Int_t p0, Int_t p1, Int_t pdgCode){
 154   //
 155   // reeturn mass pull
 156   //  Test values: p0=2; p1=2; pdgCode=321
 157   //
 158   TDatabasePDG *pdg= TDatabasePDG::Instance();
 159   Double_t pdgMass = pdg->GetParticle(pdgCode)->Mass();
 160   Double_t recMass = v0->GetKFInfo(p0,p1, 0);
 161   Double_t rmsMass = v0->GetKFInfo(p0,p1, 1);
 162   if (rmsMass<=0) return -1;
 163   return   (recMass-pdgMass)/rmsMass;
 164 }
 165
 166
 167
 168 void MakeSelectedArrays(Int_t scaling=20,  const Int_t maxEntries=10000){
 169   //
 170   //  Select pions, proton and electron based on the V0 information
 171   //  Selected tracks are written to the array in parallel with the V0 pull information
 172   //
 173   //
 174
 175
 176   TFile *fselected = TFile::Open("V0Selected.root");
 177   TTree * treeK0= (TTree*)fselected->Get("treeK0");
 178   TTree * treeLambda= (TTree*)fselected->Get("treeLambda");
 179   TTree * treeALambda= (TTree*)fselected->Get("treeALambda");
 180   TTree * treeGamma= (TTree*)fselected->Get("treeGamma");
 181   Int_t entriesK0= treeK0->GetEntries();
 182   Int_t entriesLambda= treeLambda->GetEntries();
 183   Int_t entriesALambda= treeALambda->GetEntries();
 184   Int_t entriesGamma= treeGamma->GetEntries();
 185
 186   //overestimated size
 187   arrayPions=new TObjArray(entriesK0*2+entriesLambda+entriesALambda);
 188   arrayProtons=new TObjArray(entriesLambda+entriesALambda);
 189   arrayElectrons=new TObjArray(entriesGamma*2);
 190   //
 191   ppullElectronsV0=new TMatrixD(entriesGamma*2,2);
 192   ppullPionsV0=new TMatrixD(entriesK0*2+entriesLambda+entriesALambda,2);
 193   ppullProtonsV0=new TMatrixD(entriesLambda+entriesALambda,2);
 194   TMatrixD & pullElectronsV0=*ppullElectronsV0;
 195   TMatrixD & pullPionsV0=*ppullPionsV0;
 196   TMatrixD & pullProtonsV0=*ppullProtonsV0;
 197   //
 198   // Use K0s, Lambda and ALambda to select pions
 199   //
 200   AliESDtrack * track0=0, *track1=0;
 201   AliESDv0 * v0=0;
 202   treeK0->SetBranchAddress("track0.",&track0);
 203   treeK0->SetBranchAddress("track1.",&track1);
 204   treeK0->SetBranchAddress("v0.",&v0);
 205   //
 206   treeLambda->SetBranchAddress("track0.",&track0);
 207   treeLambda->SetBranchAddress("track1.",&track1);
 208   treeLambda->SetBranchAddress("v0.",&v0);
 209   //
 210   treeALambda->SetBranchAddress("track0.",&track1);
 211   treeALambda->SetBranchAddress("track1.",&track0);
 212   treeALambda->SetBranchAddress("v0.",&v0);
 213   //
 214   treeGamma->SetBranchAddress("track0.",&track1);
 215   treeGamma->SetBranchAddress("track1.",&track0);
 216   treeGamma->SetBranchAddress("v0.",&v0);
 217   //
 218   Int_t counterPions=0;
 219   Int_t counterProtons=0;
 220   Int_t counterElectrons=0;
 221
 222   TH2F *hisPionPt = new TH2F("hisPionPt","hisPionPt", 20,0.2,20,20,-1,1);
 223   TH2F *hisProtonPt = new TH2F("hisProtonPt","hisProtonPt", 20,0.2,20,20,-1,1);
 224   TH2F *hisElectronPt = new TH2F("hisElectronPt","hisElectronPt", 20,0.2,20,20,-1,1);
 225   //
 226   AliTPCcalibBase::BinLogX(hisPionPt->GetXaxis());
 227   AliTPCcalibBase::BinLogX(hisProtonPt->GetXaxis());
 228   AliTPCcalibBase::BinLogX(hisElectronPt->GetXaxis());
 229   //
 230   //
 231   // Filter K0s
 232   //
 233   //
 234   for (Int_t iv0=0; iv0<entriesK0/scaling; iv0++){
 235     treeK0->GetEntry(iv0);
 236     Double_t pullK0= GetPullMass(v0,2,2,kK0);
 237     //
 238     hisPionPt->Fill(track0->Pt(),track0->GetTgl());
 239     Int_t bin =  hisPionPt->FindBin(track0->Pt(),track0->GetTgl());
 240     Int_t entriesPt=(bin>0)? hisPionPt->GetBinContent(bin):0;
 241     if (entriesPt< maxEntries) {
 242       pullPionsV0(counterPions,0) =pullK0;
 243       arrayPions->AddAt(track0->Clone(),counterPions); counterPions++;
 244     }
 245
 246     hisPionPt->Fill(track1->Pt(), track1->GetTgl());
 247     bin =  hisPionPt->FindBin(track1->Pt(),track1->GetTgl());
 248     entriesPt=(bin>0)? hisPionPt->GetBinContent(bin):0;
 249     if (entriesPt< maxEntries) {
 250       pullPionsV0(counterPions,0) =pullK0;
 251       arrayPions->AddAt(track1->Clone(),counterPions);
 252       counterPions++;
 253     }
 254     if (iv0%200==0) printf("K0s %d\t%d\n",iv0, counterPions);
 255   }
 256   //
 257   // Filter lambdas
 258   //
 259   for (Int_t ilambda=0; ilambda<=1; ilambda++){
 260     TTree * treeL = (ilambda==0)? treeLambda: treeALambda;
 261     Int_t entriesTree=treeL->GetEntries();
 262     for (Int_t iv0=0; iv0<entriesTree/scaling; iv0++){
 263       treeL->GetEntry(iv0);
 264       Double_t pullLambda= 0;
 265       if (ilambda==0) pullLambda=GetPullMass(v0,4,2,kLambda0);
 266       if (ilambda==0) pullLambda=GetPullMass(v0,2,4,kLambda0);
 267       // 1 Additional cut on Gammas
 268       if (TMath::Abs(v0->GetEffMass(0,0))<0.05) continue;
 269       //
 270       hisPionPt->Fill(track1->Pt(), track1->GetTgl());
 271       Int_t bin =  hisPionPt->FindBin(track1->Pt(),track1->GetTgl());
 272       Int_t entriesPt=(bin>0)? hisPionPt->GetBinContent(bin):0;
 273       if (entriesPt< maxEntries) {
 274         arrayPions->AddAt(track1->Clone(),counterPions);
 275         pullPionsV0(counterPions,0) =pullLambda;
 276         counterPions++;
 277       }
 278       //
 279       hisProtonPt->Fill(track0->Pt(), track0->GetTgl());
 280       bin =  hisProtonPt->FindBin(track0->Pt(),track0->GetTgl());
 281       entriesPt=(bin>0)? hisProtonPt->GetBinContent(bin):0;
 282       if (entriesPt< maxEntries) {
 283         arrayProtons->AddAt(track0->Clone(),counterProtons);
 284         pullProtonsV0(counterProtons,0) =pullLambda;
 285         counterProtons++;
 286       }
 287       if (iv0%200==0) printf("Lambda Protons %d\t%d\t%d\n",ilambda, iv0,counterProtons);
 288     }
 289   }
 290   //
 291   // Electrons selection
 292   //
 293
 294   Int_t nel= treeGamma->Draw("track0.fTPCsignal","abs(track0.fTPCsignal/track1.fTPCsignal-1)<0.05","goff",10000);
 295   Double_t meanElectron0=0, sigmaElectron0=0, meanElectron1 =0 , sigmaElectron1=0;
 296   AliMathBase::EvaluateUni(nel, treeGamma->GetV1(),meanElectron0,sigmaElectron0,0.55*nel);
 297   nel= treeGamma->Draw("track1.fTPCsignal",Form("abs(track0.fTPCsignal-%f)<%f&&abs(track0.fTPCsignal/track1.fTPCsignal-1)<0.15",meanElectron0,sigmaElectron0),"goff",10000);
 298   AliMathBase::EvaluateUni(nel, treeGamma->GetV1(),meanElectron1,sigmaElectron1,0.55*nel);
 299   //
 300   for (Int_t iv0=0; iv0<entriesGamma/scaling; iv0++){
 301     treeGamma->GetEntry(iv0);
 302     Double_t pullGamma= GetPullMass(v0,0,0,kGamma);
 303     //
 304     //
 305     hisElectronPt->Fill(track0->Pt(), track0->GetTgl());
 306     Int_t bin =  hisElectronPt->FindBin(track0->Pt(),track0->GetTgl());
 307     Int_t entriesPt=(bin>0)? hisElectronPt->GetBinContent(bin):0;
 308     if (entriesPt< maxEntries && TMath::Abs(track1->GetTPCsignal()-meanElectron1)<sigmaElectron1) {
 309       arrayElectrons->AddAt(track0->Clone(),counterElectrons);
 310       pullElectronsV0(counterElectrons,0) =pullGamma;
 311       counterElectrons++;
 312     }
 313     hisElectronPt->Fill(track1->Pt(), track1->GetTgl());
 314     bin =  hisElectronPt->FindBin(track1->Pt(),track1->GetTgl());
 315     entriesPt=(bin>0)? hisElectronPt->GetBinContent(bin):0;
 316     if (entriesPt< maxEntries  && TMath::Abs(track0->GetTPCsignal()-meanElectron1)<sigmaElectron1) {
 317       arrayElectrons->AddAt(track1->Clone(),counterElectrons);
 318       pullElectronsV0(counterElectrons,0) =pullGamma;
 319       counterElectrons++;
 320     }
 321     if (iv0%100==0) printf("Electrons %d\t%d\n",iv0, counterElectrons);
 322   }
 323
 324
 325
 326 }
 327
 328 void SelectHPT(){
 329   //
 330   Int_t counter=0;
 331   AliESDtrack * track0=0;
 332   TChain * chainHPT  = AliXRDPROOFtoolkit::MakeChainRandom("highptAll.list","highPt",0,1000);
 333   Int_t entriesHPT = chainHPT->GetEntries();
 334   chainHPT->SetCacheSize(1000000000);
 335   chainHPT->SetBranchAddress("esdTrack.",&track0);
 336   TBranch * branchFlags = chainHPT->GetBranch("esdTrack.fFlags");
 337   //  TBranch * branchTOFSignal = chainHPT->GetBranch("esdTrack.fITSsignal");
 338   TBranch * branchTOFSignalDX = chainHPT->GetBranch("esdTrack.fTOFsignalDx");
 339   TBranch * branchParam = chainHPT->GetBranch("esdTrack.AliExternalTrackParam.fP[5]");
 340   TBranch * branchTrack = chainHPT->GetBranch("esdTrack.");
 341   counter=0;
 342
 343   for (Int_t i=0; i<entriesHPT; i++){
 344     //
 345     if (i%100000==0) printf("Entrye\t %d\t%f\n",i,counter);
 346     branchFlags->GetEntry(i);
 347     if (track0->IsOn(0x4)==0) continue;
 348     //
 349     branchTOFSignalDX->GetEntry(i);
 350     if (track0->GetTOFsignalDx()>3) continue;
 351     //
 352     branchParam->GetEntry(i);
 353     if (TMath::Abs(track0->GetParameter()[4])<0.25) continue;
 354     branchTrack->GetEntry(i);
 355     // Select pions
 356     //esdTrack.fITSncls>4&&abs(esdTrack.fTOFsignalDx**2+esdTrack.fTOFsignalDz**2)<3&&esdTrack.fTOFsignal<esdTrack.fTrackTime[3]-400&&abs(esdTrack.fTOFsignal-esdTrack.fTrackTime[2])<300"
 357     //
 358     // Select Protons
 359     // chainHPT->Draw("esdTrack.fTPCsignal:esdTrack.P()>>his(60,0.3,6,50,30,200)","esdTrack.fITSncls>4&&abs(esdTrack.fTOFsignalDx**2+esdTrack.fTOFsignalDz**2)<4*(1+abs(esdTrack.fP[4]))&&esdTrack.fTOFsignal>esdTrack.fTrackTime[3]+500&&abs(esdTrack.fTOFsignal-esdTrack.fTrackTime[4])<4*(100+abs(esdTrack.fP[4])*100)","colz",10000000);
 360     counter++;
 361     if (counter%100==0) printf("%d\t%d\t%f\n",i, counter,track0->GetTOFsignal());
 362   }
 363 }
 364
 365
 366
 367
 368
 369
 370
 371 void FitSelected( const Int_t ngener=1000, const Int_t maxEntries=10000, Int_t year=0, Int_t period=0){
 372   //
 373   // 1. fit the roubust means and the
 374   //
 375   Int_t counter=0;
 376   TVectorD vecPt0(maxEntries);
 377   TVectorD vecP(maxEntries);
 378   TVectorD vecTheta(maxEntries);
 379   TVectorD vecdEdxTheor(maxEntries);
 380   //
 381   TVectorD vecdEdx0(maxEntries);
 382   TVectorD vecdEdx1(maxEntries);
 383   TVectorD vecdEdx2(maxEntries);
 384   TVectorD vecdEdxTRD(maxEntries);
 385   TVectorD vecdEdxOROC(maxEntries);
 386   TVectorD vecdEdxAll(maxEntries);
 387   TVectorD * vecArray[6]={&vecdEdx0, &vecdEdx1, &vecdEdx2, &vecdEdxTRD, &vecdEdxOROC, &vecdEdxAll};
 388   TVectorD vecMeanTrunc60(6);
 389   TVectorD vecRMSTrunc60(6);
 390   TVectorD vecMeanTrunc80(6);
 391   TVectorD vecRMSTrunc80(6);
 392   TVectorD vecMeanTrunc(6);
 393   TVectorD vecRMSTrunc(6);
 394   TVectorD vecMeanTruncAll(6);
 395   TVectorD vecRMSTruncAll(6);
 396   TVectorD vecSelected(6);
 397   //
 398   //
 399   //
 400   TDatabasePDG pdg;
 401   Double_t massPion = pdg.GetParticle("pi+")->Mass();
 402   Double_t massProton = pdg.GetParticle("proton")->Mass();
 403   Double_t massElectron = pdg.GetParticle("e+")->Mass();
 404
 405   TTreeSRedirector * pcstream = new TTreeSRedirector("pidDEDX.root","recreate");
 406   //
 407   //
 408   //
 409   for (Int_t isel=0; isel<ngener; isel++){
 410     Double_t mcp=gRandom->Rndm()*20;
 411     Double_t mctheta=2*(gRandom->Rndm()-0.5);
 412     Int_t ipid= TMath::Nint(gRandom->Rndm()*3.);
 413     Double_t mass=massPion;
 414     TObjArray *arrayTrack = arrayPions;
 415     if (ipid==1) {arrayTrack=arrayProtons; mass=massProton;}
 416     if (ipid==2) {arrayTrack=arrayElectrons; mass=massElectron;}
 417     Int_t entriesArray= arrayTrack->GetEntries();
 418     //
 419     Int_t nselected=0;
 420     Int_t nselectedArray[6]={0};
 421     TVectorD vselectedArray(6);
 422     //
 423     // 1. First loop - select the tracks which fullfill the pt and theta criteria
 424     //
 425     for (Int_t itrack=0; itrack<entriesArray; itrack++){
 426       AliESDtrack * track = (AliESDtrack*)arrayTrack->At(itrack);
 427       if (!track) continue;
 428       if (!track->GetInnerParam()) continue;
 429       if (TMath::Abs(track->P()/mcp-1)>0.15) continue; // pt selection
 430       if (TMath::Abs(track->GetTgl()-mctheta)>0.15) continue; // theta selection
 431       Double_t mom = track->GetInnerParam()->P();
 432       Double_t dEdxTheor=50*AliMathBase::BetheBlochAleph(mom/mass);
 433       //
 434       vecPt0[nselected]=track->Pt();
 435       vecP[nselected]=mom;
 436       vecTheta[nselected]=track->GetTgl();
 437       vecdEdxTheor[nselected]=dEdxTheor;
 438       //
 439       if (track->GetTPCdEdxInfo()){
 440         vecdEdx0[nselectedArray[0]++]=track->GetTPCdEdxInfo()->GetSignal(0)/dEdxTheor;
 441         vecdEdx1[nselectedArray[1]++]=track->GetTPCdEdxInfo()->GetSignal(1)/dEdxTheor;
 442         vecdEdx2[nselectedArray[2]++]=track->GetTPCdEdxInfo()->GetSignal(2)/dEdxTheor;
 443         vecdEdxOROC[nselectedArray[4]++]=track->GetTPCdEdxInfo()->GetSignal(3)/dEdxTheor;
 444       }
 445       if (track->GetTRDsignal()>0 && track->GetTRDncls()>80){
 446         vecdEdxTRD[nselectedArray[3]++]=50.*track->GetTRDsignal()/dEdxTheor;
 447       }
 448       vecdEdxAll[nselectedArray[5]++]=track->GetTPCsignal()/dEdxTheor;
 449       nselected++;
 450     }
 451     for (Int_t idedx=0; idedx<6; idedx++) vselectedArray[idedx]=nselectedArray[idedx];
 452     if (nselected <20) continue;
 453     //
 454     Double_t meanPt=TMath::Mean(nselected,vecPt0.GetMatrixArray());
 455     Double_t meanP=TMath::Mean(nselected,vecP.GetMatrixArray());
 456     Double_t meanTheta=TMath::Mean(nselected,vecTheta.GetMatrixArray());
 457     Double_t meanTheor=TMath::Mean(nselected,vecdEdxTheor.GetMatrixArray());
 458     Double_t meanRob, rmsRob;
 459     //
 460     for (Int_t idet=0; idet<6; idet++){
 461       vecMeanTrunc60[idet]=0;
 462       vecRMSTrunc60[idet]=0;
 463       vecMeanTrunc80[idet]=0;
 464       vecRMSTrunc80[idet]=0;
 465       vecMeanTruncAll[idet]=0;
 466       vecRMSTruncAll[idet]=0;
 467       if  (nselectedArray[idet]<20) continue;
 468       AliMathBase::EvaluateUni( nselectedArray[idet], vecArray[idet]->GetMatrixArray(), meanRob, rmsRob, nselectedArray[idet]*0.6);
 469       vecMeanTrunc60[idet]=meanRob;
 470       vecRMSTrunc60[idet]=rmsRob;
 471       AliMathBase::EvaluateUni( nselectedArray[idet], vecArray[idet]->GetMatrixArray(), meanRob, rmsRob, nselectedArray[idet]*0.8);
 472       vecMeanTrunc80[idet]=meanRob;
 473       vecRMSTrunc80[idet]=rmsRob;
 474       vecMeanTruncAll[idet]= TMath::Mean(nselectedArray[idet],vecArray[idet]->GetMatrixArray());
 475       vecRMSTruncAll[idet] = TMath::RMS(nselectedArray[idet],vecArray[idet]->GetMatrixArray());
 476     }
 477     //
 478     if (isel%10==0) printf("%d\n",isel);
 479     Double_t dEdxTheorMeanP = AliMathBase::BetheBlochAleph(meanP/mass);
 480     Double_t dEdxTheorCenterP = AliMathBase::BetheBlochAleph(mcp/mass);
 481     (*pcstream)<<"pid"<<
 482       "ipid="<<ipid<<                     // pid type -0 pion, 1- proton, 2 electron
 483       "nAll="<<nselected<<                // number of primary points
 484       "year="<<year<<                     // year
 485       "period="<<period<<                 // period
 486       //
 487       "dEdxTheorMean="<<meanTheor<<                 // initial dEdx hypothesis - mean over given bin
 488       "dEdxTheorMeanP="<<dEdxTheorMeanP<<           // initial dEdx hypothesis - for mean momenta
 489       "dEdxTheorCenterP="<<dEdxTheorCenterP<<       // initial dEdx hypothesis - for central momenta
 490       // bin
 491       "p="<<mcp<<                                   // center of bin for particle momentum
 492       "theta="<<mctheta<<                           // particle theta
 493       //
 494       "mass="<<mass<<                               // mass of particle
 495       "meanPt="<<meanPt<<                           // mean pt in bin
 496       "meanP="<<meanP<<
 497       "meanTheta="<<meanTheta<<
 498       //
 499       "dedxMean60.="<<&vecMeanTrunc60<<
 500       "dedxRMS60.="<<&vecRMSTrunc60<<
 501       "dedxMean80.="<<&vecMeanTrunc80<<
 502       "dedxRMS80.="<<&vecRMSTrunc80<<
 503       "dedxMeanAll.="<<&vecMeanTruncAll<<
 504       "dedxRMSAll.="<<&vecRMSTruncAll<<
 505       "\n";
 506     //
 507     // Test
 508     //
 509   }
 510   delete pcstream;
 511
 512 }
 513
 514
 515
 516
 517 /*
 518
 519
 520
 521 void  AnalyzedEdx(){
 522   //
 523   pcstream = new TTreeSRedirector("trend.root","recreate");
 524   TString flist = gSystem->GetFromPipe("cat highpt.list");
 525   TObjArray * array = flist.Tokenize("\n");
 526   array->Print();
 527   Int_t entries = array->GetEntries();
 528   for (Int_t ifile=0; ifile<entries; ifile++){
 529     //TString name = array->At(ifile)->GetName();
 530     printf("\n\n\n");
 531     printf("Analyzing:\t%d\t%s\n", ifile,array->At(ifile)->GetName());
 532     AnalyzedEdxFile(Form("%s#FilterEvents_Trees.root",array->At(ifile)->GetName()));
 533     printf("\n\n\n");
 534   }
 535   delete pcstream;
 536 }
 537
 538 void AnalyzedEdxFile(const char * fname){
 539   //
 540   //scomment
 541     const char * fname="/hera/alice/local/filtered/alice/data/2013/LHC13e/000196201/vpass1/root_archive.zip#FilterEvents_Trees.root";
 542     const char * fname="/hera/alice/local/filtered/alice/data/2013/LHC13c/000195531/ESDs/pass1/root_archive.zip#FilterEvents_Trees.root";
 543   //ecomment
 544   //
 545   TFile * fin=TFile::Open(fname);
 546   if (!fin){
 547     printf("File\t%s  not existing\n",fname);
 548     return ;
 549   }
 550   TTree * treeIn = (TTree*)fin->Get("highPt");
 551   if (!fin){
 552     printf("Tree \t%s  not existing\n",fname);
 553     return ;
 554   }
 555   treeIn->SetCacheSize(1000000000.0);
 556   // TFile *fout = TFile::Open("aaa.root","recreate");
 557
 558   tree =treeIn->CopyTree("esdTrack.fITSncls>0||(esdTrack.fTOFsignal<100000-1)","",50000);
 559   if (!tree) return;
 560   if (tree->GetEntries()<2000) return;
 561   SetPionAliases();
 562   FitSlopeROC();
 563 }
 564
 565 void SetPionAliases(){
 566   //
 567   // Select the tracks with higher dEdx - at high momenta it is used in order to get clean sample of Pions
 568   // IROC selection used to select sample for OROC
 569   // OROC selection used to select sample for IROC
 570   //
 571   Int_t entries = 0;
 572   tree->SetAlias("pionOROC","(1+0)");
 573   tree->SetAlias("pionIROC","(1+0)");
 574
 575   Double_t mean=0, sigma=0;
 576   entries = tree->Draw("runNumber","","goff",1000);
 577   run = tree->GetV1()[0];
 578   {for (Int_t iter=0; iter<2; iter++){
 579     entries = tree->Draw("esdTrack.fTPCdEdxInfo.fTPCsignalRegion[3]/AliMathBase::BetheBlochAleph(esdTrack.P()/0.1)>>his","esdTrack.fTPCdEdxInfo.fTPCsignalRegion[3]>0&&pionIROC", "goff");
 580     AliMathBase::EvaluateUni(entries, tree->GetV1(), mean,sigma, 0.6*entries);
 581     tree->SetAlias("pionOROC",Form("esdTrack.fTPCdEdxInfo.fTPCsignalRegion[3]/AliMathBase::BetheBlochAleph(esdTrack.P()/0.1)-%f>%f",mean,0.0*sigma));
 582     entries = tree->Draw("esdTrack.fTPCdEdxInfo.fTPCsignalRegion[0]/AliMathBase::BetheBlochAleph(esdTrack.P()/0.1)>>his","esdTrack.fTPCdEdxInfo.fTPCsignalRegion[0]>0&&pionOROC", "goff");
 583     AliMathBase::EvaluateUni(entries, tree->GetV1(), mean,sigma, 0.6*entries);
 584     tree->SetAlias("pionIROC",Form("esdTrack.fTPCdEdxInfo.fTPCsignalRegion[0]/AliMathBase::BetheBlochAleph(esdTrack.P()/0.1)-%f>%f",mean,0.0*sigma));
 585     printf("%d\n",iter);
 586     printf("%s\n",tree->GetAlias("pionIROC"));
 587     printf("%s\n",tree->GetAlias("pionOROC"));
 588     }}
 589
 590 }
 591
 592
 593 void FitSlopeROC(){
 594   //
 595   // fit dEdx/eta for "identified high pt pions"
 596   //
 597   TF1 f1IROC("f1IROC","pol1");
 598   TF1 f1OROC("f1IROC","pol1");
 599   TVectorD vectorIROC(2);
 600   TVectorD vectorOROC(2);
 601   TH2 * hisOROC=0;
 602   TH2 * hisIROC=0;
 603   TObjArray fitIROC(3);
 604   TObjArray fitOROC(3);
 605   tree->Draw("esdTrack.fTPCdEdxInfo.fTPCsignalRegion[3]/AliMathBase::BetheBlochAleph(esdTrack.P()/0.1):abs(esdTrack.fP[3])>>hisOROC(5,0,1,100,20,80)","esdTrack.fTPCncls>70&&pionIROC","colz");
 606   hisOROC = (TH2*)tree->GetHistogram()->Clone();
 607   tree->Draw("esdTrack.fTPCdEdxInfo.fTPCsignalRegion[0]/AliMathBase::BetheBlochAleph(esdTrack.P()/0.1):abs(esdTrack.fP[3])>>hisIROC(5,0,1,100,20,80)","esdTrack.fTPCncls>70&&pionOROC","colz");
 608   hisIROC = (TH2*)tree->GetHistogram()->Clone();
 609   hisOROC->FitSlicesY(0,0,-1,0,"QNR",&fitOROC);
 610   hisIROC->FitSlicesY(0,0,-1,0,"QNR",&fitIROC);
 611   TH1 * phisOROC[3];
 612   TH1 * phisIROC[3];
 613   for (Int_t ihis=0; ihis<3; ihis++){
 614     phisOROC[ihis]=(TH1*)fitOROC.At(ihis);
 615     phisIROC[ihis]=(TH1*)fitIROC.At(ihis);
 616   }
 617   //
 618   phisOROC[1]->SetMarkerColor(2);
 619   phisOROC[1]->SetMarkerStyle(25);
 620   phisIROC[1]->SetMarkerColor(2);
 621   phisIROC[1]->SetMarkerStyle(25);
 622   phisOROC[1]->Fit(&f1OROC);
 623   f1OROC.GetParameters(vectorOROC.GetMatrixArray());
 624   phisOROC[1]->Draw();
 625   phisIROC[1]->Fit(&f1IROC);
 626   f1IROC.GetParameters(vectorIROC.GetMatrixArray());
 627   phisIROC[1]->Draw();
 628   //
 629   (*pcstream)<<"dedxFit"<<
 630     "run="<<run<<
 631     "hisIROC1.="<<phisIROC[1]<<  // IROC - dEdx/dEdxExp - mean as function of theta
 632     "hisIROC2.="<<phisIROC[2]<<  // IROC - dEdx/dEdxExp -        RMS as function of the theta
 633     "hisOROC1.="<<phisOROC[1]<<  // OROC - dEdx/dEdxExp - mean as function of theta
 634     "hisOROC2.="<<phisOROC[2]<<  // OROC - dEdx/dEdxExp -        RMS as function of theta
 635     "fIROC.="<<&vectorIROC<<     // IROC vector with fit parameters mean (at 0 theta) and splope IROC
 636     "fOROC.="<<&vectorOROC<<     // OROC vector
 637     "f1IROC.="<<&f1IROC<<
 638     "f1OROC.="<<&f1OROC<<
 639     "\n";
 640
 641   //scomment
 642   //resolution at given bin divided
 643   //
 644   Standard plots:
 645   TStatToolkit::MakeGraphSparse(chain,"hisOROC2.fArray[1]/hisOROC1.fArray[1]:run","",25,1,0)->Draw("alp");
 646   TStatToolkit::MakeGraphSparse(chain,"hisIROC2.fArray[1]/hisIROC1.fArray[1]:run","",25,1,0)->Draw("alp")
 647   //
 648   TStatToolkit::MakeGraphSparse(chain,"fIROC.fElements[0]+fIROC.fElements[1]*0.5:run","",25,1,0)->Draw("alp");
 649
 650   //ecomment
 651
 652 }
 653
 654
 655 void DrawdEdx(){
 656   //
 657   //pdg.GetParticle("Lambda0")->Mass
 658   //
 659   TDatabasePDG pdg;
 660   Double_t massLambda = pdg.GetParticle("Lambda0")->Mass();
 661   Double_t massK0 = pdg.GetParticle("K0")->Mass();
 662   Double_t massPion = pdg.GetParticle("pi+")->Mass();
 663   Double_t massProton = pdg.GetParticle("proton")->Mass();
 664
 665   chain->SetAlias("massPion",Form("(%f+0)",massPion));
 666   chain->SetAlias("massProton",Form("(%f+0)",massProton));
 667   chain->SetAlias("massK0",Form("(%f+0)",massK0));
 668   chain->SetAlias("massLambda",Form("(%f+0)",massLambda));
 669   chain->SetAlias("K0Pull","(v0.GetEffMass(2,2)-massK0)/v0.GetKFInfo(2,2,1)");
 670   chain->SetAlias("LPull","(v0.GetEffMass(4,2)-massLambda)/v0.GetKFInfo(4,2,1)");
 671   chain->SetAlias("K0sel", "abs(v0.GetEffMass(2,2)-massK0)<min(abs(v0.GetEffMass(4,2)-massLambda), abs(v0.GetEffMass(2,4)-massLambda))&&abs(K0Pull)<3");
 672   //
 673   //
 674   //
 675   chain->Draw("track0.fTPCsignal/AliMathBase::BetheBlochAleph(track0.fIp.P()/massPion)>>hisP(50,20,80)","K0sel","");
 676   chain->Draw("track1.fTPCsignal/AliMathBase::BetheBlochAleph(track1.fIp.P()/massPion)>>+hisN(50,20,80)","K0sel","");
 677
 678   chain->Draw("track0.fTPCsignal/AliMathBase::BetheBlochAleph(track0.fIp.P()/massPion):track0.fIp.P()>>his(20,0,10,100,30,70)","K0sel","colz");
 679
 680
 681   TH3F * hisBGPion3D = new TH3F("hisBGPion3D","hisBGPion3D",40,0.5,200,4,0,1,100,30,70);
 682   TH3F * hisBGProton3D = new TH3F("hisBGProton3D","hisBGProton3D",40,0.5,200,4,0,1,100,30,70);
 683   AliTPCcalibBase::BinLogX(hisBGPion3D->GetXaxis());
 684   AliTPCcalibBase::BinLogX(hisBGProton3D->GetXaxis());
 685
 686   chain->Draw("track0.fTPCsignal/AliMathBase::BetheBlochAleph(track0.fIp.P()/massPion):abs(track0.fIp.fP[3]):track0.fIp.P()/massPion>>hisBGPion3D","K0sel||(track0.fTOFr[2]+track0.fTOFr[1])>0.6","");
 687   chain->Draw("track0.fTPCsignal/AliMathBase::BetheBlochAleph(track0.fIp.P()/massProton):abs(track0.fIp.fP[3]):track0.fIp.P()/massProton>>hisBGProton3D","(!K0sel&&abs(LPull)<1.5)||track0.fTOFr[4]>0.6","");
 688
 689
 690   TTreeSRedirector *pcstream = new TTreeSRedirector("bb.root","recreate");
 691   TVectorD vecPionMean(4),vecPionRMS(4), vecPionEntries(4);
 692   TVectorD vecProtonMean(4),vecProtonRMS(4), vecProtonEntries(4);
 693   TVectorD vecTheta(4);
 694   {for (Int_t ibg=2; ibg<38; ibg++){
 695     for (Int_t itheta=0; itheta<4; itheta++){
 696       TH1 * hisPion = hisBGPion3D->ProjectionZ("hisPion",ibg-1, ibg+1, itheta+1, itheta+1);
 697       TH1 * hisProton = hisBGProton3D->ProjectionZ("hisProton",ibg-1, ibg+1, itheta+1, itheta+1);
 698       hisPion->Fit("gaus");
 699       vecPionMean[itheta]=gaus->GetParameter(1);
 700       vecPionRMS[itheta]=gaus->GetParameter(2);
 701       vecPionEntries[itheta]=hisPion->GetEntries();
 702       hisProton->Fit("gaus");
 703       vecProtonMean[itheta]=gaus->GetParameter(1);
 704       vecProtonRMS[itheta]=gaus->GetParameter(2);
 705       vecProtonEntries[itheta]=hisProton->GetEntries();
 706       vecTheta[itheta]=hisBGPion3D->GetXaxis()->GetBinCenter(itheta);
 707       delete hisPion;
 708       delete hisProton;
 709     }
 710     //
 711     Double_t bg = hisBGPion3D->GetXaxis()->GetBinCenter(ibg);
 712     Double_t dEdxRef = AliMathBase::BetheBlochAleph(bg);
 713     (*pcstream)<<"bb"<<
 714       "bg="<<bg<<
 715       "dedxRef="<<dEdxRef<<
 716       "ibg="<<ibg<<
 717       "pionMean.="<<&vecPionMean<<
 718       "pionRMS.="<<&vecPionRMS<<
 719       "pionEntries.="<<&vecPionEntries<<
 720       "protonMean.="<<&vecProtonMean<<
 721       "protonRMS.="<<&vecProtonRMS<<
 722       "protonEntries.="<<&vecProtonEntries<<
 723       "\n";
 724     }}
 725
 726   delete pcstream;
 727
 728
 729
 730
 731   hisBGPion->FitSlicesY();
 732   hisBGProton->FitSlicesY();
 733
 734   hisBGPion_1->SetMarkerStyle(25);
 735   hisBGProton_1->SetMarkerStyle(25);
 736   hisBGPion_1->SetMarkerColor(2);
 737   hisBGProton_1->SetMarkerColor(4);
 738
 739   hisBGPion_1->Draw();
 740   hisBGProton_1->Draw("same");
 741
 742
 743
 744 }
 745
 746
 747
 748 void aaa(){
 749
 750   TChain * chain  = AliXRDPROOFtoolkit::MakeChainRandom("/hera/alice/miranov/highptAll.list","V0s",0,1000);
 751   chain->SetCacheSize(1000000000);
 752
 753   TChain * chainHPT  = AliXRDPROOFtoolkit::MakeChainRandom("/hera/alice/miranov/highptAll.list","highPt",0,1000);
 754   chainHPT->SetCacheSize(1000000000);
 755   // for Kaons not good referenece data
 756   TFile ftofSelected("tofSelected.root","recreate");
 757   TTree * treeTOF = chainHPT->CopyTree("abs(esdTrack.fTOFsignalDz)<5&&abs(esdTrack.fTOFsignalDx)<5&&esdTrack.fITSncls>0&&abs(1/esdTrack.fP[4])<2");
 758   treeTOF->Write("treeTOF");
 759   //ftofSelected.Write();
 760
 761 }
 762
 763
 764
 765
 766
 767
 768 void DrawdEdxRatio(){
 769   //
 770   //pdg.GetParticle("Lambda0")->Mass
 771   //
 772   TDatabasePDG pdg;
 773   Double_t massLambda = pdg.GetParticle("Lambda0")->Mass();
 774   Double_t massK0 = pdg.GetParticle("K0")->Mass();
 775   Double_t massPion = pdg.GetParticle("pi+")->Mass();
 776   Double_t massProton = pdg.GetParticle("proton")->Mass();
 777
 778   chain->SetAlias("massPion",Form("(%f+0)",massPion));
 779   chain->SetAlias("massProton",Form("(%f+0)",massProton));
 780   chain->SetAlias("massK0",Form("(%f+0)",massK0));
 781   chain->SetAlias("massLambda",Form("(%f+0)",massLambda));
 782   chain->SetAlias("K0Pull","(v0.GetEffMass(2,2)-massK0)/v0.GetKFInfo(2,2,1)");
 783   chain->SetAlias("LPull","(v0.GetEffMass(4,2)-massLambda)/v0.GetKFInfo(4,2,1)");
 784   chain->SetAlias("K0sel", "abs(v0.GetEffMass(2,2)-massK0)<min(abs(v0.GetEffMass(4,2)-massLambda), abs(v0.GetEffMass(2,4)-massLambda))&&abs(K0Pull)<3");
 785
 786
 787   //
 788   //
 789   //
 790
 791   TH3F * hisBGPion3DR1 = new TH3F("hisBGPion3DR1","hisBGPion3DR1",40,0.5,200,4,0,1,100,0.5,1.5);
 792   TH3F * hisBGProton3DR1 = new TH3F("hisBGProton3DR1","hisBGProton3DR1",40,0.5,200,4,0,1,100,0.5,1.5);
 793   TH3F * hisBGPion3DR2 = new TH3F("hisBGPion3DR2","hisBGPion3DR2",40,0.5,200,4,0,1,100,0.5,1.5);
 794   TH3F * hisBGProton3DR2 = new TH3F("hisBGProton3DR2","hisBGProton3DR2",40,0.5,200,4,0,1,100,0.5,1.5);
 795   AliTPCcalibBase::BinLogX(hisBGPion3DR1->GetXaxis());
 796   AliTPCcalibBase::BinLogX(hisBGProton3DR1->GetXaxis());
 797   AliTPCcalibBase::BinLogX(hisBGPion3DR2->GetXaxis());
 798   AliTPCcalibBase::BinLogX(hisBGProton3DR2->GetXaxis());
 799
 800   chain->Draw("track0.fTPCdEdxInfo.fTPCsignalRegion[1]/track0.fTPCdEdxInfo.fTPCsignalRegion[0]:abs(track0.fIp.fP[3]):track0.fIp.P()/massPion>>hisBGPion3DR1","K0sel||(track0.fTOFr[2]+track0.fTOFr[1])>0.6","");
 801   chain->Draw("track0.fTPCdEdxInfo.fTPCsignalRegion[2]/track0.fTPCdEdxInfo.fTPCsignalRegion[0]:abs(track0.fIp.fP[3]):track0.fIp.P()/massPion>>hisBGPion3DR2","K0sel||(track0.fTOFr[2]+track0.fTOFr[1])>0.6","");
 802   //
 803   //
 804   //
 805   chain->Draw("track0.fTPCdEdxInfo.fTPCsignalRegion[1]/track0.fTPCdEdxInfo.fTPCsignalRegion[0]:abs(track0.fIp.fP[3]):track0.fIp.P()/massProton>>hisBGProton3DR1","(!K0sel&&abs(LPull)<1.5)||track0.fTOFr[4]>0.6","");
 806   chain->Draw("track0.fTPCdEdxInfo.fTPCsignalRegion[2]/track0.fTPCdEdxInfo.fTPCsignalRegion[0]:abs(track0.fIp.fP[3]):track0.fIp.P()/massProton>>hisBGProton3DR2","(!K0sel&&abs(LPull)<1.5)||track0.fTOFr[4]>0.6","");
 807
 808   TFile fratio("bbRatio.root","recreate");
 809   hisBGPion3DR1->Write();
 810   hisBGPion3DR2->Write();
 811   hisBGProton3DR1->Write();
 812   hisBGProton3DR2->Write();
 813   fratio.Close();
 814
 815   Double_t masses[5]={massPion, massProton,0,0,0};
 816   const Double_t kB2C=-0.299792458e-3;
 817
 818   AliTPCParamSR par;
 819   par.Update();
 820
 821   TTreeSRedirector *pcstream = new TTreeSRedirector("bbRatio.root","update");
 822   TVectorD vecMean(4),vecRMS(4), vecEntries(4);
 823   TVectorD vecTheta(4);
 824   TVectorD vecP(4);
 825   TVectorD vecPt(4);
 826   TVectorD vecPhi0(4);  // tracklet angle at 0
 827   TVectorD vecPhi1(4);  // tracklet angle at 1
 828   TVectorD vecPhi2(4);  // tracklet angle at 2
 829   TVectorD vecTheta0(4);  // tracklet angle at 0
 830   TVectorD vecTheta1(4);  // tracklet angle at 1
 831   TVectorD vecTheta2(4);  // tracklet angle at 2
 832   TVectorD vecL0(4);  // tracklet length at 0
 833   TVectorD vecL1(4);  // tracklet letgth at 1
 834   TVectorD vecL2(4);  // tracklet length at 2
 835   {
 836     for (Int_t ptype=0; ptype<2; ptype++){
 837       for (Int_t itype=0; itype<2; itype++){
 838         for (Int_t ibg=2; ibg<38; ibg++){
 839           Double_t bg = hisBGPion3DR1->GetXaxis()->GetBinCenter(ibg);
 840           Double_t dEdxRef = AliMathBase::BetheBlochAleph(bg);
 841           Double_t length0=0;
 842           for (Int_t itheta=0; itheta<4; itheta++){
 843             TH1 * hisProjection = 0;
 844             if (itype==0){
 845               if (ptype==0) hisProjection   = hisBGPion3DR1->ProjectionZ("hisPion",ibg-1, ibg+1, itheta+1, itheta+1);
 846               if (ptype==1) hisProjection   = hisBGProton3DR1->ProjectionZ("hisProton",ibg-1, ibg+1, itheta+1, itheta+1);
 847             }
 848             if (itype==1){
 849               if (ptype==0) hisProjection   = hisBGPion3DR2->ProjectionZ("hisPion",ibg-1, ibg+1, itheta+1, itheta+1);
 850               if (ptype==1) hisProjection   = hisBGProton3DR2->ProjectionZ("hisProton",ibg-1, ibg+1, itheta+1, itheta+1);
 851             }
 852             hisProjection->Fit("gaus");
 853             vecMean[itheta]=gaus->GetParameter(1);
 854             vecRMS[itheta]=gaus->GetParameter(2);
 855             vecEntries[itheta]=hisProjection->GetEntries();
 856             vecTheta[itheta]=hisBGPion3DR1->GetYaxis()->GetBinCenter(itheta+1);
 857             vecP[itheta]= bg*masses[ptype];
 858             vecPt[itheta]=bg*masses[ptype]/TMath::Sqrt(1+vecTheta[itheta]*vecTheta[itheta]);  //Pt
 859             //
 860             // Get angle and length
 861             //
 862             Double_t crv= 5*kB2C/vecPt[itheta];   //GetC(b); // bz*kB2C/pt;
 863             Double_t angleIROC= TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(0, par.GetNRow(0)/2.)*crv)*0.5,1.));
 864             Double_t angleOROC0= TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(36, par.GetNRowUp1()/2.)*crv)*0.5,1.));
 865             Double_t angleOROC1= TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(36, par.GetNRowUp1()+par.GetNRowUp2()/2.)*crv)*0.5,1.));
 866             vecPhi0[itheta]=angleIROC;
 867             vecPhi1[itheta]=angleOROC0;
 868             vecPhi2[itheta]=angleOROC1;
 869             //
 870             vecTheta0[itheta]=TMath::Sqrt(1+TMath::Tan(vecPhi0[itheta])*TMath::Tan(vecPhi0[itheta]))*vecTheta[itheta];
 871             vecTheta1[itheta]=TMath::Sqrt(1+TMath::Tan(vecPhi1[itheta])*TMath::Tan(vecPhi1[itheta]))*vecTheta[itheta];
 872             vecTheta2[itheta]=TMath::Sqrt(1+TMath::Tan(vecPhi2[itheta])*TMath::Tan(vecPhi2[itheta]))*vecTheta[itheta];
 873
 874             //
 875             vecL0[itheta]  = 0.75*TMath::Sqrt(1+TMath::Tan(vecPhi0[itheta])*TMath::Tan(vecPhi0[itheta]) + vecTheta0[itheta]*vecTheta0[itheta]);
 876             vecL1[itheta]  = 1.0*TMath::Sqrt(1+TMath::Tan(vecPhi1[itheta])*TMath::Tan(vecPhi1[itheta])  +  vecTheta1[itheta]*vecTheta1[itheta]);
 877             vecL2[itheta]  = 1.5*TMath::Sqrt(1+TMath::Tan(vecPhi2[itheta])*TMath::Tan(vecPhi2[itheta])  +  vecTheta2[itheta]*vecTheta2[itheta]);
 878             delete hisProjection;
 879           }
 880           //
 881           (*pcstream)<<"bbRatio"<<
 882             "ptype="<<ptype<<            // particle type
 883             "rtype="<<itype<<            // ratio type
 884             "mass="<<masses[ptype]<<     // mass of the
 885             "bg="<<bg<<
 886             "dedxRef="<<dEdxRef<<
 887             "ibg="<<ibg<<
 888             "vecTheta.="<<&vecTheta<<
 889             "vecTheta0.="<<&vecTheta0<<
 890             "vecTheta1.="<<&vecTheta1<<
 891             "vecTheta2.="<<&vecTheta2<<
 892             "ratioMean.="<<&vecMean<<
 893             "ratioRMS.="<<&vecRMS<<
 894             "ratioEntries.="<<&vecEntries<<
 895             "vecP.="<<&vecP<<             // momentum
 896             "vecPt.="<<&vecPt<<             // pt
 897             "vecPhi0.="<<&vecPhi0<<       // angle in the IROC
 898             "vecPhi1.="<<&vecPhi1<<       // angle in the OROC1
 899             "vecPhi2.="<<&vecPhi2<<       // angle in the OROC2
 900             "vecL0.="<<&vecL0<<       // angle in the IROC
 901             "vecL1.="<<&vecL1<<       // angle in the OROC1
 902             "vecL2.="<<&vecL2<<       // angle in the OROC2
 903             "\n";
 904         }}
 905     }
 906   }
 907
 908   delete pcstream;
 909
 910 }
 911
 912 void Fit(){
 913
 914   TFile ff("bbRatio.root");
 915   TTree * tree = ff.Get("bbRatio");
 916   tree->SetMarkerStyle(25);
 917   //
 918   Double_t chi2=0;
 919   Int_t    npoints=0;
 920   TVectorD param;
 921   TMatrixD covar;
 922   Int_t npointsMax=10000000;
 923
 924   TCut cutFit="ratioRMS.fElements<0.3&&ratioEntries.fElements>100&&vecL1.fElements<10&&vecL2.fElements<10";
 925   TString fstringFast="";
 926   fstringFast+="(1-(rtype==0)*2)++";
 927   fstringFast+="(rtype==0)*vecTheta.fElements++";
 928   fstringFast+="(rtype==1)*vecTheta.fElements++";
 929   //  fstringFast+="vecTheta.fElements++";
 930   //fstringFast+="vecPhi0.fElements++";
 931   fstringFast+="1/(vecL0.fElements*dedxRef)++";
 932   //fstringFast+="(rtype==0)*vecPhi1.fElements++";
 933   fstringFast+="(rtype==0)/(vecL1.fElements*dedxRef)++";
 934   //fstringFast+="(rtype==1)*vecPhi2.fElements++";
 935   fstringFast+="(rtype==1)/(vecL2.fElements*dedxRef)++";
 936   //
 937   TString *strDelta = TStatToolkit::FitPlaneConstrain(tree,"ratioMean.fElements:ratioRMS.fElements", fstringFast.Data(),cutFit, chi2,npoints,param,covar,0.9,0, npointsMax, 1);
 938   TObjArray* tokArr = strDelta->Tokenize("++");
 939   tokArr->Print();
 940   tree->SetAlias("corr",strDelta->Data());
 941 }
 942
 943
 944 void SelectTracks(){
 945   //
 946   //
 947   //
 948   TFile f("tofSelected.root");
 949   TChain * chain = f.Get("highPt");
 950   chain->SetAlias("ITSTOFrPion","sqrt((esdTrack.fTOFr[1]+esdTrack.fTOFr[2]+0.05)*(esdTrack.fITSr[1]+esdTrack.fITSr[2]+0.05))");
 951   chain->SetAlias("ITSTOFrProton","sqrt((esdTrack.fTOFr[4]+0.05)*(esdTrack.fITSr[4]+0.05))");
 952   chain->SetAlias("ITSTOFrKaon","sqrt((esdTrack.fTOFr[3]+0.05)*(esdTrack.fITSr[3]+0.05))");
 953   TDatabasePDG pdg;
 954   Double_t massLambda = pdg.GetParticle("Lambda0")->Mass();
 955   Double_t massK0 = pdg.GetParticle("K0")->Mass();
 956   Double_t massPion = pdg.GetParticle("pi+")->Mass();
 957   Double_t massProton = pdg.GetParticle("proton")->Mass();
 958
 959   chain->SetAlias("massPion",Form("(%f+0)",massPion));
 960   chain->SetAlias("massProton",Form("(%f+0)",massProton));
 961   chain->SetAlias("massK0",Form("(%f+0)",massK0));
 962
 963   //
 964   chain->SetAlias("isPion", "((ITSTOFrPion>1.5*ITSTOFrProton&&ITSTOFrPion>1.5*ITSTOFrKaon)&&esdTrack.fTRDsignal/AliMathBase::BetheBlochAleph(esdTrack.P()/0.13-1)<0.2)");
 965   //
 966   chain->SetAlias("isKaon", "((ITSTOFrKaon>1.5*ITSTOFrProton&&ITSTOFrKaon>1.5*ITSTOFrPion))&&abs(esdTrack.fTRDsignal/AliMathBase::BetheBlochAleph(esdTrack.P()/massK0)-1)<0.3");
 967
 968   chain->SetAlias("isProton", "((ITSTOFrProton>1.5*ITSTOFrKaon&&ITSTOFrProton>1.5*ITSTOFrPion))");
 969
 970 }
 971
 972 */