[u/mrichter/AliRoot.git] / TRD / AliTRDclusterErrors.C

//#ifndef __CINT__
  #include <iostream.h>

  #include "AliTRDtracker.h"
  #include "AliTRDclusterMI.h"
  #include "AliTRDhit.h"
  #include "AliTRDv1.h"
  #include "AliTRDgeometry.h"
  #include "AliTRDgeometryDetail.h"
  #include "AliTRDparameter.h"

  #include "alles.h"
  #include "TFile.h"
  #include "TStopwatch.h"
  #include "Rtypes.h"
  #include "TTree.h"
#include "AliRunLoader.h"
#include "AliStack.h"
#include "AliTRDclusterErrors.h"
#include "TF1.h"
#include "AliTrackReference.h"
//#endif    


ClassImp(AliTRDCI)
ClassImp(AliTRDExactPoint)
ClassImp(AliTRDClusterErrAnal)
ClassImp(AliTRDClusterErrDraw)


AliTRDClusterErrAnal::AliTRDClusterErrAnal(Char_t *chloader )
{
  //
  //SET Input loaders
  if (gAlice){
    delete gAlice->GetRunLoader();
    delete gAlice;
    gAlice = 0x0;
  }    
  fRunLoader = AliRunLoader::Open(chloader);
  if (fRunLoader == 0x0){
    cerr<<"Can not open session"<<endl;
    return;
  }
  fTRDLoader = fRunLoader->GetLoader("TRDLoader");  
  if (fTRDLoader == 0x0){
    cerr<<"Can not get TRD Loader"<<endl;
    return ;
  }   
  if (fRunLoader->LoadgAlice()){
    cerr<<"Error occured while l"<<endl;
    return;
  }
  fRunLoader->CdGAFile();
  fTracker = new AliTRDtracker(gFile);
  fParam = (AliTRDparameter*) gFile->Get("TRDparameter");
  fGeometry = new AliTRDgeometryDetail();   
  fTRD      = (AliTRDv1*) gAlice->GetDetector("TRD");
  //
  AliTRDCI * clinfo = new AliTRDCI();
  fFileA  = new TFile("trdclusteranal.root","recreate");
  fFileA->cd();
  fTreeA  = new TTree("trdcl","trdcl");
  fTreeA->Branch("trdcl","AliTRDCI",&clinfo);
  delete clinfo;
}

void AliTRDClusterErrAnal::SetIO(Int_t event)
{
  //
  //set input output for given event
  fRunLoader->SetEventNumber(event);
  fRunLoader->LoadHeader();
  fRunLoader->LoadKinematics();
  fRunLoader->LoadTrackRefs();
  fTRDLoader->LoadHits();
  fTRDLoader->LoadRecPoints("read");
  //
  fStack = fRunLoader->Stack();
  fHitTree = fTRDLoader->TreeH();
  fClusterTree = fTRDLoader->TreeR();
  fReferenceTree = fRunLoader->TreeTR();
  //
}

void AliTRDClusterErrAnal::LoadClusters()
{
  //
  //
  // Load clusters  
  TObjArray *ClusterArray = new TObjArray(400);   
  TObjArray carray(2000);
  TBranch *branch=fClusterTree->GetBranch("TRDcluster");
  if (!branch) {
    Error("ReadClusters","Can't get the branch !");
    return;
  }
  branch->SetAddress(&ClusterArray);
  Int_t nentries = (Int_t)fClusterTree->GetEntries();
  for (Int_t i=0;i<nentries;i++){
    fClusterTree->GetEvent(i);
    Int_t nCluster = ClusterArray->GetEntriesFast();
    for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
      AliTRDcluster* c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
      carray.AddLast(c);
      ClusterArray->RemoveAt(iCluster);
    }
  }
  Int_t nClusters = carray.GetEntriesFast();
  printf("Total number of clusters %d \n", nClusters);
  //
  //
  //SORT clusters
  //
  Int_t all=0;  
  for (Int_t i=0;i<nClusters;i++){
    AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i);  
    Int_t plane = fGeometry->GetPlane(cl->GetDetector());
    if (plane>=12) continue;
    Int_t time = cl->GetLocalTimeBin(); 
    if (time>=100) continue;
    Int_t sector = fGeometry->GetSector(cl->GetDetector());
    if (sector>=18){
      printf("problem1\n");
      continue;
    }    
    fClusters[plane][time][sector].AddLast(cl);
    all++;
  }
}

void AliTRDClusterErrAnal::SortReferences()
{
  //
  //
  //
  printf("Sorting references\n");
  fReferences = new TObjArray;
  Int_t ntracks = fStack->GetNtrack();
  fReferences->Expand(ntracks);
  Int_t nentries = (Int_t)fReferenceTree->GetEntries();
  TClonesArray * arr = new TClonesArray("AliTrackReference");
  TBranch * br = fReferenceTree->GetBranch("TRD");
  br->SetAddress(&arr);
  //
  TClonesArray *labarr=0;
  for (Int_t iprim=0;iprim<nentries;iprim++){
    if (br->GetEntry(iprim)){
      for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){
	AliTrackReference *ref =(AliTrackReference*)arr->At(iref);
	if (!ref) continue;
	Int_t lab = ref->GetTrack();
	if ( (lab<0) || (lab>ntracks)) continue;
	//
	if (fReferences->At(lab)==0) {
	  labarr = new TClonesArray("AliTrackReference"); 
	  fReferences->AddAt(labarr,lab);
	}
	TClonesArray &larr = *labarr;
	new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
      }
    }
  }
  printf("end - Sorting references\n");
}

AliTrackReference * AliTRDClusterErrAnal::FindNearestReference(Int_t lab, Float_t pos[3], Float_t dmax)
{
  //
  //
  //
  if (fReferences->At(lab)==0) return 0;
  AliTrackReference *nearest=0;
  TClonesArray * arr = (TClonesArray *)fReferences->At(lab);
  for (Int_t iref =0;iref<arr->GetEntriesFast();iref++){
    AliTrackReference * ref = ( AliTrackReference *)arr->UncheckedAt(iref);
    if (!ref) continue;
    Float_t delta = (pos[0]-ref->X())*(pos[0]-ref->X());
    delta += (pos[1]-ref->Y())*(pos[1]-ref->Y());
    delta += (pos[2]-ref->Z())*(pos[2]-ref->Z());
    delta = TMath::Sqrt(delta);
    if (delta<dmax){
      dmax=delta;
      nearest = ref;
    }
  }
  return nearest;
}

void AliTRDClusterErrAnal::MakeExactPoints(Int_t trackmax)
{
  //
  //make exact points:-)	
  //
  // 

  fExactPoints.Delete();
  fExactPoints.Expand(fStack->GetNtrack());
  //
  Double_t fSum=0;
  Double_t fSumQ =0;
  Double_t fSumX=0;
  Double_t fSumX2=0;
  Double_t fSumXY=0;
  Double_t fSumXZ=0;
  Double_t fSumY=0;
  Double_t fSumZ=0;
  //
  Int_t entries = Int_t(fHitTree->GetEntries());
  printf("Number of primary  entries\t%d\n",entries);
  entries = TMath::Min(trackmax,entries);
  Int_t nallpoints = 0;

  Int_t nalltracks =0;
  Int_t pointspertrack =0;

  for (Int_t entry=0;entry<entries; entry++){
    gAlice->ResetHits();
    fHitTree->GetEvent(entry);
    Int_t lastlabel    = -1;
    Int_t lastdetector = -1;
    Int_t lasttimebin   = -1;
    Float_t lastpos[3];
    //
    for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1); hit; 
	hit = (AliTRDhit *) fTRD->NextHit()) {
      //
      Int_t label    = hit->Track();
      TParticle * particle = fStack->Particle(label);
      if (!particle) continue;
      if (particle->Pt()<0.05) continue;
      Int_t detector = hit->GetDetector();
      Int_t plane    = fGeometry->GetPlane(detector);
      //
      //
      if (hit->GetCharge()==0) continue;
      Float_t pos[3] = {hit->X(),hit->Y(),hit->Z()};
      Int_t indexes[3];
      fGeometry->Global2Detector(pos,indexes,fParam);
      //
      Float_t rot[3];
      fGeometry->Rotate(detector,pos,rot);
      //
      //
      Float_t  time0    = fParam->GetTime0(plane);
      Int_t    timebin  = Int_t(TMath::Nint(((time0 - rot[0])/fParam->GetTimeBinSize())+ fParam->GetTimeBefore())+0.1); 
      if (timebin<0) continue;
      //
      //
      if (label!=lastlabel || detector != lastdetector || lasttimebin !=timebin){
	//
	if (label!=lastlabel){
	  fExactPoints.AddAt(new TClonesArray("AliTRDExactPoint",0),label);
	  //printf("new particle\t%d\n",label);
	  nalltracks++;
	  //	  printf("particle\t%d- hits\t%d\n",lastlabel, pointspertrack);
	  pointspertrack=0;
	}
	
	if ( (fSum>1) && lasttimebin>=0 && lasttimebin<fParam->GetTimeMax() ){
	  //if we have enough info for given layer time bin - store it
	  AliTrackReference * ref = FindNearestReference(lastlabel,lastpos,4.);
	  Float_t rotmom[3];
	  Float_t rotpos[3];
	  Float_t refangle=0;
	  if (ref){
	    Float_t mom[3] = {ref->Px(),ref->Py(),ref->Pz()};
	    Float_t refpos[3] = {ref->X(),ref->Y(),ref->Z()};
	    fGeometry->Rotate(detector,mom,rotmom);
	    fGeometry->Rotate(detector,refpos,rotpos);
	    refangle = rotmom[1]/rotmom[0];

	  }

	  Double_t ay,by,az,bz;
	  Double_t det = fSum*fSumX2-fSumX*fSumX;
	  if (TMath::Abs(det)> 0.000000000000001) { 
	    by = (fSum*fSumXY-fSumX*fSumY)/det;
	    ay = (fSumX2*fSumY-fSumX*fSumXY)/det;
	    
	  }else{
	    ay =fSumXY/fSumX;
	    by =0;	   
	  }
	  if (TMath::Abs(det)> 0.000000000000001) { 
	    bz = (fSum*fSumXZ-fSumX*fSumZ)/det;
	    az = (fSumX2*fSumZ-fSumX*fSumXZ)/det;	  
	  }else{
	    az =fSumXZ/fSumX;
	    bz =0;	   
	  }
	  //
	  Float_t lastplane = fGeometry->GetPlane(lastdetector);
	  Float_t time0    = fParam->GetTime0(lastplane);
	  Float_t xcenter = time0 - (lasttimebin - fParam->GetTimeBefore())*fParam->GetTimeBinSize();
	  Float_t ty = ay + by * xcenter;
	  Float_t tz = az + bz * xcenter;
	  //

	  TClonesArray * arr = (TClonesArray *) fExactPoints.At(label);
	  TClonesArray & larr= *arr;
	  Int_t arrent = arr->GetEntriesFast();
	  AliTRDExactPoint * point = new (larr[arrent]) AliTRDExactPoint;
	  nallpoints++;

	  if (ref){
	    point->SetReference(ref);
	    point->fTRefAngleY = rotmom[1]/rotmom[0];
	  }
	  point->fTX  = xcenter;
	  point->fTY  = ty;
	  point->fTZ  = tz;
	  point->fTAY = by;
	  point->fTAZ = bz;
	  //
	  point->fGx  = lastpos[0];
	  point->fGy  = lastpos[1];
	  point->fGz  = lastpos[2];

	  //
	  point->fDetector     = lastdetector;
	  point->fLocalTimeBin = lasttimebin;
	  point->fPlane        = fGeometry->GetPlane(lastdetector); 
	  point->fSector      = fGeometry->GetSector(lastdetector); 
	  point->fPlaneMI      = indexes[0]; 
	  //
	  point->fTPrim = fSum;
	  point->fTQ    = fSumQ;	    
	  //
	}
	lastdetector = detector;
	lastlabel    = label;
	lasttimebin  = timebin;
	fSum=fSumQ=fSumX=fSumX2=fSumXY=fSumXZ=fSumY=fSumZ=0.;
      }
      //
      lastpos[0] = hit->X();
      lastpos[1] = hit->Y();
      lastpos[2] = hit->Z();
      fSum++;
      fSumQ = hit->GetCharge();
      fSumX  +=rot[0];
      fSumX2 +=rot[0]*rot[0];
      fSumXY +=rot[0]*rot[1];      
      fSumXZ +=rot[0]*rot[2];
      fSumY  +=rot[1];      
      fSumZ  +=rot[2];     
      pointspertrack++;
    }
  }
  //
  printf("Found %d exact points\n",nallpoints); 
}


Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) {  

  //
  // comparison works with both cluster types MI and old also
  //dummy cluster to be fill if not cluster info
  AliTRDclusterMI clmi;  
  TClass * classmi =  clmi.IsA();
  //
  //SetOutput
  AliTRDCI * clinfo = new AliTRDCI();
  TBranch * clbr = fTreeA->GetBranch("trdcl");
  clbr->SetAddress(&clinfo);

  SetIO(0);
  SortReferences();
  MakeExactPoints(trackmax);
  LoadClusters();
  //
  trackmax =  fStack->GetNtrack();
  //
  // Get the number of entries in the hit tree
  // (Number of primary particles creating a hit somewhere)
  Int_t nTrack = (Int_t)fExactPoints.GetEntries();
  printf("Found %d charged in TRD in first %d primaries", nTrack, trackmax);
  //

  for (Int_t itrack = 0; itrack<trackmax; itrack++){
    TClonesArray *arrpoints = (TClonesArray*)fExactPoints.At(itrack);
    
    if (!arrpoints) continue;
    //printf("new particle\t%d\n",itrack);
    TParticle * particle = fStack->Particle(itrack);
    if (!particle) continue;
    //printf("founded in kine tree \t%d\n",itrack);
    Int_t npoints = arrpoints->GetEntriesFast();
    if (npoints<10) continue;
    //printf("have enough points \t%d\t%d\n",itrack,npoints);

    for (Int_t ipoint=0;ipoint<npoints;ipoint++){
      AliTRDExactPoint * point = (AliTRDExactPoint*)arrpoints->UncheckedAt(ipoint);
      if (!point) continue;
      //
      Int_t sec = fGeometry->GetSector(point->fDetector);
      if (sec>18){
	printf("problem2\n");
      }
      TObjArray & cllocal = fClusters[point->fPlane][point->fLocalTimeBin][sec];       
      Int_t nclusters = cllocal.GetEntriesFast();
      Float_t maxdist = 10;
      AliTRDcluster * nearestcluster =0;
      //find nearest cluster to hit with given label
      for (Int_t icluster =0; icluster<nclusters; icluster++){
	AliTRDcluster * cluster = (AliTRDcluster*)cllocal.UncheckedAt(icluster);
	if (!cluster) continue;
	if ( (cluster->GetLabel(0)!= itrack) &&  (cluster->GetLabel(1)!= itrack)&&(cluster->GetLabel(2)!= itrack))
	  continue;
	Float_t dist = TMath::Abs(cluster->GetY()-point->fTY);
	if (TMath::Abs(cluster->GetZ()-point->fTZ)>5.5) continue; 
	if (dist<maxdist){
	  maxdist = dist;
	  nearestcluster = cluster;
	}
      }      
      //
      clinfo->fEp  = *point;
      clinfo->fP   = *particle;
      if (!nearestcluster) {
	clinfo->fStatus=1;
	clinfo->fCl = clmi;
      }
      else{
	clinfo->fStatus=0;
	if (nearestcluster->IsA()==classmi){
	  clinfo->fCl    =*((AliTRDclusterMI*)nearestcluster);    
	}
	else{
	  clinfo->fCl    = *nearestcluster;
	}
	//     
	Float_t dz      = clinfo->fCl.GetZ()-point->fTZ;
	Double_t h01    = sin(TMath::Pi() / 180.0 * fParam->GetTiltingAngle());
	clinfo->fTDistZ = dz;
	clinfo->fDYtilt = h01*dz*((point->fPlane%2)*2.-1.); 
	//
	clinfo->fNTracks =1;
	if (nearestcluster->GetLabel(1)>=0)  clinfo->fNTracks++;
	if (nearestcluster->GetLabel(2)>=0)  clinfo->fNTracks++;  
	clinfo->Update();
      }
      //
      fTreeA->Fill();
    }
  }
  
  
  fFileA->cd();
  fTreeA->Write();
  fFileA->Close();
  return 0;
}


TH1F * AliTRDClusterErrDraw::ResDyVsAmp(TTree* tree, const char* selection, Float_t t0, Float_t ampmin, Float_t ampmax)
{
  //
  //
  TH2F hisdy("resy","resy",10,ampmin,ampmax,30,-0.3,0.3);
  char expression[1000];
  sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fQ>>resy",t0);
  char selectionall[1000];
  sprintf(selectionall,"fStatus==0&&%s",selection);
  printf("%s\n",expression);
  printf("%s\n",selectionall);
  tree->Draw(expression,selectionall);
  return CreateResHisto(&hisdy);
}


TH1F * AliTRDClusterErrDraw::ResDyVsRelPos(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
{
  //
  //
  min *=128;
  max *=128;
  TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
  char expression[1000];
  sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fRelPos>>resy",t0);
  char selectionall[1000];
  sprintf(selectionall,"fStatus==0&&%s",selection);
  printf("%s\n",expression);
  printf("%s\n",selectionall);
  tree->Draw(expression,selectionall);
  return CreateResHisto(&hisdy);
}


TH1F * AliTRDClusterErrDraw::ResDyVsAngleY(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
{
  //
  //
  TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);

  char expression[1000];
  sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%f*fEp.fTAY:fEp.fTAY>>resy",t0);
  char selectionall[1000];
  sprintf(selectionall,"fStatus==0&&%s",selection);

  tree->Draw(expression,selectionall);
  return CreateResHisto(&hisdy);
}

void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
{
  histo->GetXaxis()->SetTitle(xAxisTitle);
  histo->GetYaxis()->SetTitle(yAxisTitle);
}


TH1F* AliTRDClusterErrDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
{
  Int_t nBins = hGen->GetNbinsX();
  TH1F* hEff = (TH1F*) hGen->Clone("hEff");
  hEff->SetTitle("");
  hEff->SetStats(kFALSE);
  hEff->SetMinimum(0.);
  hEff->SetMaximum(110.);

  for (Int_t iBin = 0; iBin <= nBins; iBin++) {
    Double_t nGen = hGen->GetBinContent(iBin);
    Double_t nRec = hRec->GetBinContent(iBin);
    if (nGen > 0) {
      Double_t eff = nRec/nGen;
      hEff->SetBinContent(iBin, 100. * eff);
      Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);      
      //      if (error == 0) error = sqrt(0.1/nGen);
      //
      if (error == 0) error = 0.0001;
      hEff->SetBinError(iBin, 100. * error);
    } else {
      hEff->SetBinContent(iBin, 100. * 0.5);
      hEff->SetBinError(iBin, 100. * 0.5);
    }
  }

  return hEff;
}


TH1F* AliTRDClusterErrDraw::CreateResHisto(TH2F* hRes2, Bool_t draw,  Bool_t drawBinFits, 
		     Bool_t overflowBinFits)
{
  TVirtualPad* currentPad = gPad;
  TAxis* axis = hRes2->GetXaxis();
  Int_t nBins = axis->GetNbins();
  TH1F* hRes, *hMean;
  if (axis->GetXbins()->GetSize()){
    hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
    hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
  }
  else{
    hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
    hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());

  }
  hRes->SetStats(false);
  hRes->SetOption("E");
  hRes->SetMinimum(0.);
  //
  hMean->SetStats(false);
  hMean->SetOption("E");
 
  // create the fit function
  //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
  TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
  
  fitFunc->SetLineWidth(2);
  fitFunc->SetFillStyle(0);
  // create canvas for fits
  TCanvas* canBinFits = NULL;
  Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
  Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
  Int_t ny = (nPads-1) / nx + 1;
  if (drawBinFits) {
    canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
    if (canBinFits) delete canBinFits;
    canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
    canBinFits->Divide(nx, ny);
  }

  // loop over x bins and fit projection
  Int_t dBin = ((overflowBinFits) ? 1 : 0);
  for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
    if (drawBinFits) canBinFits->cd(bin + dBin);
    TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
    //    
    if (hBin->GetEntries() > 10) {
      fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS(),0.02*hBin->GetMaximum());
      hBin->Fit(fitFunc,"s");
      Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));

      if (sigma > 0.){
	hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
	hMean->SetBinContent(bin, fitFunc->GetParameter(1));	
      }
      else{
	hRes->SetBinContent(bin, 0.);
	hMean->SetBinContent(bin,0);
      }
      hRes->SetBinError(bin, fitFunc->GetParError(2));
      hMean->SetBinError(bin, fitFunc->GetParError(1));
      
      //
      //

    } else {
      hRes->SetBinContent(bin, 0.);
      hRes->SetBinError(bin, 0.);
      hMean->SetBinContent(bin, 0.);
      hMean->SetBinError(bin, 0.);
    }
    

    if (drawBinFits) {
      char name[256];
      if (bin == 0) {
	sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
      } else if (bin == nBins+1) {
	sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
      } else {
	sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
		axis->GetTitle(), axis->GetBinUpEdge(bin));
      }
      canBinFits->cd(bin + dBin);
      hBin->SetTitle(name);
      hBin->SetStats(kTRUE);
      hBin->DrawCopy("E");
      canBinFits->Update();
      canBinFits->Modified();
      canBinFits->Update();
    }
    
    delete hBin;
  }

  delete fitFunc;
  currentPad->cd();
  if (draw){
    currentPad->Divide(1,2);
    currentPad->cd(1);
    hRes->Draw();
    currentPad->cd(2);
    hMean->Draw();
  }
  
  return hRes;
}
Commit	Line	Data
289478c7	1	//#ifndef __CINT__
59dfcadd	2	#include <iostream.h>
	3
	4	#include "AliTRDtracker.h"
289478c7	5	#include "AliTRDclusterMI.h"
59dfcadd	6	#include "AliTRDhit.h"
	7	#include "AliTRDv1.h"
	8	#include "AliTRDgeometry.h"
289478c7	9	#include "AliTRDgeometryDetail.h"
59dfcadd	10	#include "AliTRDparameter.h"
	11
	12	#include "alles.h"
	13	#include "TFile.h"
	14	#include "TStopwatch.h"
289478c7	15	#include "Rtypes.h"
	16	#include "TTree.h"
	17	#include "AliRunLoader.h"
	18	#include "AliStack.h"
	19	#include "AliTRDclusterErrors.h"
	20	#include "TF1.h"
	21	#include "AliTrackReference.h"
	22	//#endif
	23
	24
	25	ClassImp(AliTRDCI)
	26	ClassImp(AliTRDExactPoint)
	27	ClassImp(AliTRDClusterErrAnal)
	28	ClassImp(AliTRDClusterErrDraw)
	29
	30
	31
	32	AliTRDClusterErrAnal::AliTRDClusterErrAnal(Char_t *chloader )
	33	{
	34	//
	35	//SET Input loaders
	36	if (gAlice){
	37	delete gAlice->GetRunLoader();
	38	delete gAlice;
	39	gAlice = 0x0;
	40	}
	41	fRunLoader = AliRunLoader::Open(chloader);
	42	if (fRunLoader == 0x0){
	43	cerr<<"Can not open session"<<endl;
	44	return;
	45	}
	46	fTRDLoader = fRunLoader->GetLoader("TRDLoader");
	47	if (fTRDLoader == 0x0){
	48	cerr<<"Can not get TRD Loader"<<endl;
	49	return ;
	50	}
	51	if (fRunLoader->LoadgAlice()){
	52	cerr<<"Error occured while l"<<endl;
	53	return;
	54	}
	55	fRunLoader->CdGAFile();
	56	fTracker = new AliTRDtracker(gFile);
	57	fParam = (AliTRDparameter*) gFile->Get("TRDparameter");
	58	fGeometry = new AliTRDgeometryDetail();
	59	fTRD = (AliTRDv1*) gAlice->GetDetector("TRD");
	60	//
	61	AliTRDCI * clinfo = new AliTRDCI();
	62	fFileA = new TFile("trdclusteranal.root","recreate");
	63	fFileA->cd();
	64	fTreeA = new TTree("trdcl","trdcl");
	65	fTreeA->Branch("trdcl","AliTRDCI",&clinfo);
	66	delete clinfo;
	67	}
59dfcadd	68
289478c7	69	void AliTRDClusterErrAnal::SetIO(Int_t event)
	70	{
	71	//
	72	//set input output for given event
	73	fRunLoader->SetEventNumber(event);
	74	fRunLoader->LoadHeader();
	75	fRunLoader->LoadKinematics();
	76	fRunLoader->LoadTrackRefs();
	77	fTRDLoader->LoadHits();
	78	fTRDLoader->LoadRecPoints("read");
	79	//
	80	fStack = fRunLoader->Stack();
	81	fHitTree = fTRDLoader->TreeH();
	82	fClusterTree = fTRDLoader->TreeR();
	83	fReferenceTree = fRunLoader->TreeTR();
	84	//
	85	}
59dfcadd	86
289478c7	87	void AliTRDClusterErrAnal::LoadClusters()
	88	{
	89	//
	90	//
	91	// Load clusters
	92	TObjArray *ClusterArray = new TObjArray(400);
59dfcadd	93	TObjArray carray(2000);
289478c7	94	TBranch *branch=fClusterTree->GetBranch("TRDcluster");
	95	if (!branch) {
	96	Error("ReadClusters","Can't get the branch !");
	97	return;
	98	}
	99	branch->SetAddress(&ClusterArray);
	100	Int_t nentries = (Int_t)fClusterTree->GetEntries();
	101	for (Int_t i=0;i<nentries;i++){
	102	fClusterTree->GetEvent(i);
	103	Int_t nCluster = ClusterArray->GetEntriesFast();
	104	for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
	105	AliTRDcluster* c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
	106	carray.AddLast(c);
	107	ClusterArray->RemoveAt(iCluster);
	108	}
	109	}
59dfcadd	110	Int_t nClusters = carray.GetEntriesFast();
59dfcadd	111	printf("Total number of clusters %d \n", nClusters);
289478c7	112	//
	113	//
	114	//SORT clusters
	115	//
	116	Int_t all=0;
	117	for (Int_t i=0;i<nClusters;i++){
	118	AliTRDcluster cl = (AliTRDcluster ) carray.UncheckedAt(i);
	119	Int_t plane = fGeometry->GetPlane(cl->GetDetector());
	120	if (plane>=12) continue;
	121	Int_t time = cl->GetLocalTimeBin();
	122	if (time>=100) continue;
	123	Int_t sector = fGeometry->GetSector(cl->GetDetector());
	124	if (sector>=18){
	125	printf("problem1\n");
	126	continue;
	127	}
	128	fClusters[plane][time][sector].AddLast(cl);
	129	all++;
59dfcadd	130	}
289478c7	131	}
59dfcadd	132
289478c7	133	void AliTRDClusterErrAnal::SortReferences()
	134	{
	135	//
	136	//
	137	//
	138	printf("Sorting references\n");
	139	fReferences = new TObjArray;
	140	Int_t ntracks = fStack->GetNtrack();
	141	fReferences->Expand(ntracks);
	142	Int_t nentries = (Int_t)fReferenceTree->GetEntries();
	143	TClonesArray * arr = new TClonesArray("AliTrackReference");
	144	TBranch * br = fReferenceTree->GetBranch("TRD");
	145	br->SetAddress(&arr);
	146	//
	147	TClonesArray *labarr=0;
	148	for (Int_t iprim=0;iprim<nentries;iprim++){
	149	if (br->GetEntry(iprim)){
	150	for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){
	151	AliTrackReference ref =(AliTrackReference)arr->At(iref);
	152	if (!ref) continue;
	153	Int_t lab = ref->GetTrack();
	154	if ( (lab<0) \|\| (lab>ntracks)) continue;
	155	//
	156	if (fReferences->At(lab)==0) {
	157	labarr = new TClonesArray("AliTrackReference");
	158	fReferences->AddAt(labarr,lab);
	159	}
	160	TClonesArray &larr = *labarr;
	161	new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
	162	}
	163	}
	164	}
	165	printf("end - Sorting references\n");
	166	}
59dfcadd	167
289478c7	168	AliTrackReference * AliTRDClusterErrAnal::FindNearestReference(Int_t lab, Float_t pos[3], Float_t dmax)
	169	{
	170	//
	171	//
	172	//
	173	if (fReferences->At(lab)==0) return 0;
	174	AliTrackReference *nearest=0;
	175	TClonesArray * arr = (TClonesArray *)fReferences->At(lab);
	176	for (Int_t iref =0;iref<arr->GetEntriesFast();iref++){
	177	AliTrackReference * ref = ( AliTrackReference *)arr->UncheckedAt(iref);
	178	if (!ref) continue;
	179	Float_t delta = (pos[0]-ref->X())*(pos[0]-ref->X());
	180	delta += (pos[1]-ref->Y())*(pos[1]-ref->Y());
	181	delta += (pos[2]-ref->Z())*(pos[2]-ref->Z());
	182	delta = TMath::Sqrt(delta);
	183	if (delta<dmax){
	184	dmax=delta;
	185	nearest = ref;
	186	}
	187	}
	188	return nearest;
	189	}
59dfcadd	190
289478c7	191	void AliTRDClusterErrAnal::MakeExactPoints(Int_t trackmax)
	192	{
	193	//
	194	//make exact points:-)
	195	//
	196	//
	197
	198	fExactPoints.Delete();
	199	fExactPoints.Expand(fStack->GetNtrack());
	200	//
	201	Double_t fSum=0;
	202	Double_t fSumQ =0;
	203	Double_t fSumX=0;
	204	Double_t fSumX2=0;
	205	Double_t fSumXY=0;
	206	Double_t fSumXZ=0;
	207	Double_t fSumY=0;
	208	Double_t fSumZ=0;
	209	//
	210	Int_t entries = Int_t(fHitTree->GetEntries());
	211	printf("Number of primary entries\t%d\n",entries);
	212	entries = TMath::Min(trackmax,entries);
	213	Int_t nallpoints = 0;
	214
	215	Int_t nalltracks =0;
	216	Int_t pointspertrack =0;
	217
	218	for (Int_t entry=0;entry<entries; entry++){
	219	gAlice->ResetHits();
	220	fHitTree->GetEvent(entry);
	221	Int_t lastlabel = -1;
	222	Int_t lastdetector = -1;
	223	Int_t lasttimebin = -1;
	224	Float_t lastpos[3];
	225	//
	226	for(AliTRDhit hit = (AliTRDhit ) fTRD->FirstHit(-1); hit;
	227	hit = (AliTRDhit *) fTRD->NextHit()) {
	228	//
	229	Int_t label = hit->Track();
	230	TParticle * particle = fStack->Particle(label);
	231	if (!particle) continue;
	232	if (particle->Pt()<0.05) continue;
	233	Int_t detector = hit->GetDetector();
	234	Int_t plane = fGeometry->GetPlane(detector);
	235	//
	236	//
	237	if (hit->GetCharge()==0) continue;
	238	Float_t pos[3] = {hit->X(),hit->Y(),hit->Z()};
	239	Int_t indexes[3];
	240	fGeometry->Global2Detector(pos,indexes,fParam);
	241	//
	242	Float_t rot[3];
	243	fGeometry->Rotate(detector,pos,rot);
	244	//
	245	//
	246	Float_t time0 = fParam->GetTime0(plane);
	247	Int_t timebin = Int_t(TMath::Nint(((time0 - rot[0])/fParam->GetTimeBinSize())+ fParam->GetTimeBefore())+0.1);
	248	if (timebin<0) continue;
	249	//
	250	//
	251	if (label!=lastlabel \|\| detector != lastdetector \|\| lasttimebin !=timebin){
	252	//
	253	if (label!=lastlabel){
	254	fExactPoints.AddAt(new TClonesArray("AliTRDExactPoint",0),label);
255	//printf("new particle\t%d\n",label);
256	nalltracks++;
257	// printf("particle\t%d- hits\t%d\n",lastlabel, pointspertrack);
258	pointspertrack=0;
259	}
260
261	if ( (fSum>1) && lasttimebin>=0 && lasttimebin<fParam->GetTimeMax() ){
262	//if we have enough info for given layer time bin - store it
263	AliTrackReference * ref = FindNearestReference(lastlabel,lastpos,4.);
264	Float_t rotmom[3];
265	Float_t rotpos[3];
266	Float_t refangle=0;
267	if (ref){
268	Float_t mom[3] = {ref->Px(),ref->Py(),ref->Pz()};
269	Float_t refpos[3] = {ref->X(),ref->Y(),ref->Z()};
270	fGeometry->Rotate(detector,mom,rotmom);
271	fGeometry->Rotate(detector,refpos,rotpos);
272	refangle = rotmom[1]/rotmom[0];
273
274	}
275
276	Double_t ay,by,az,bz;
277	Double_t det = fSumfSumX2-fSumXfSumX;
278	if (TMath::Abs(det)> 0.000000000000001) {
279	by = (fSumfSumXY-fSumXfSumY)/det;
280	ay = (fSumX2fSumY-fSumXfSumXY)/det;
281
282	}else{
283	ay =fSumXY/fSumX;
284	by =0;
285	}
286	if (TMath::Abs(det)> 0.000000000000001) {
287	bz = (fSumfSumXZ-fSumXfSumZ)/det;
288	az = (fSumX2fSumZ-fSumXfSumXZ)/det;
289	}else{
290	az =fSumXZ/fSumX;
291	bz =0;
292	}
293	//
294	Float_t lastplane = fGeometry->GetPlane(lastdetector);
295	Float_t time0 = fParam->GetTime0(lastplane);
296	Float_t xcenter = time0 - (lasttimebin - fParam->GetTimeBefore())*fParam->GetTimeBinSize();
297	Float_t ty = ay + by * xcenter;
298	Float_t tz = az + bz * xcenter;
299	//
300
301	TClonesArray * arr = (TClonesArray *) fExactPoints.At(label);
302	TClonesArray & larr= *arr;
303	Int_t arrent = arr->GetEntriesFast();
304	AliTRDExactPoint * point = new (larr[arrent]) AliTRDExactPoint;
305	nallpoints++;
306
307	if (ref){
308	point->SetReference(ref);
309	point->fTRefAngleY = rotmom[1]/rotmom[0];
310	}
311	point->fTX = xcenter;
312	point->fTY = ty;
313	point->fTZ = tz;
314	point->fTAY = by;
315	point->fTAZ = bz;
316	//
317	point->fGx = lastpos[0];
318	point->fGy = lastpos[1];
319	point->fGz = lastpos[2];
320
321	//
322	point->fDetector = lastdetector;
323	point->fLocalTimeBin = lasttimebin;
324	point->fPlane = fGeometry->GetPlane(lastdetector);
325	point->fSector = fGeometry->GetSector(lastdetector);
326	point->fPlaneMI = indexes[0];
327	//
328	point->fTPrim = fSum;
329	point->fTQ = fSumQ;
330	//
331	}
332	lastdetector = detector;
333	lastlabel = label;
334	lasttimebin = timebin;
335	fSum=fSumQ=fSumX=fSumX2=fSumXY=fSumXZ=fSumY=fSumZ=0.;
336	}
337	//
338	lastpos[0] = hit->X();
339	lastpos[1] = hit->Y();
340	lastpos[2] = hit->Z();
341	fSum++;
342	fSumQ = hit->GetCharge();
343	fSumX +=rot[0];
344	fSumX2 +=rot[0]*rot[0];
345	fSumXY +=rot[0]*rot[1];
346	fSumXZ +=rot[0]*rot[2];
347	fSumY +=rot[1];
348	fSumZ +=rot[2];
349	pointspertrack++;
350	}
351	}
352	//
353	printf("Found %d exact points\n",nallpoints);
354	}
59dfcadd	355
59dfcadd	356
59dfcadd	357
59dfcadd	358
59dfcadd	359
59dfcadd	360
289478c7	361	Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) {
59dfcadd	362
289478c7	363	//
	364	// comparison works with both cluster types MI and old also
	365	//dummy cluster to be fill if not cluster info
	366	AliTRDclusterMI clmi;
	367	TClass * classmi = clmi.IsA();
	368	//
	369	//SetOutput
	370	AliTRDCI * clinfo = new AliTRDCI();
	371	TBranch * clbr = fTreeA->GetBranch("trdcl");
	372	clbr->SetAddress(&clinfo);
59dfcadd	373
289478c7	374	SetIO(0);
	375	SortReferences();
	376	MakeExactPoints(trackmax);
	377	LoadClusters();
	378	//
	379	trackmax = fStack->GetNtrack();
	380	//
	381	// Get the number of entries in the hit tree
	382	// (Number of primary particles creating a hit somewhere)
	383	Int_t nTrack = (Int_t)fExactPoints.GetEntries();
	384	printf("Found %d charged in TRD in first %d primaries", nTrack, trackmax);
	385	//
	386
	387	for (Int_t itrack = 0; itrack<trackmax; itrack++){
	388	TClonesArray arrpoints = (TClonesArray)fExactPoints.At(itrack);
	389
	390	if (!arrpoints) continue;
	391	//printf("new particle\t%d\n",itrack);
	392	TParticle * particle = fStack->Particle(itrack);
	393	if (!particle) continue;
	394	//printf("founded in kine tree \t%d\n",itrack);
	395	Int_t npoints = arrpoints->GetEntriesFast();
	396	if (npoints<10) continue;
	397	//printf("have enough points \t%d\t%d\n",itrack,npoints);
	398
	399	for (Int_t ipoint=0;ipoint<npoints;ipoint++){
	400	AliTRDExactPoint * point = (AliTRDExactPoint*)arrpoints->UncheckedAt(ipoint);
	401	if (!point) continue;
	402	//
	403	Int_t sec = fGeometry->GetSector(point->fDetector);
	404	if (sec>18){
	405	printf("problem2\n");
	406	}
	407	TObjArray & cllocal = fClusters[point->fPlane][point->fLocalTimeBin][sec];
	408	Int_t nclusters = cllocal.GetEntriesFast();
	409	Float_t maxdist = 10;
	410	AliTRDcluster * nearestcluster =0;
	411	//find nearest cluster to hit with given label
	412	for (Int_t icluster =0; icluster<nclusters; icluster++){
	413	AliTRDcluster * cluster = (AliTRDcluster*)cllocal.UncheckedAt(icluster);
	414	if (!cluster) continue;
	415	if ( (cluster->GetLabel(0)!= itrack) && (cluster->GetLabel(1)!= itrack)&&(cluster->GetLabel(2)!= itrack))
	416	continue;
	417	Float_t dist = TMath::Abs(cluster->GetY()-point->fTY);
	418	if (TMath::Abs(cluster->GetZ()-point->fTZ)>5.5) continue;
	419	if (dist<maxdist){
	420	maxdist = dist;
	421	nearestcluster = cluster;
	422	}
	423	}
	424	//
	425	clinfo->fEp = *point;
	426	clinfo->fP = *particle;
	427	if (!nearestcluster) {
	428	clinfo->fStatus=1;
	429	clinfo->fCl = clmi;
	430	}
	431	else{
	432	clinfo->fStatus=0;
	433	if (nearestcluster->IsA()==classmi){
	434	clinfo->fCl =((AliTRDclusterMI)nearestcluster);
	435	}
	436	else{
	437	clinfo->fCl = *nearestcluster;
59dfcadd	438	}
289478c7	439	//
	440	Float_t dz = clinfo->fCl.GetZ()-point->fTZ;
	441	Double_t h01 = sin(TMath::Pi() / 180.0 * fParam->GetTiltingAngle());
	442	clinfo->fTDistZ = dz;
	443	clinfo->fDYtilt = h01dz((point->fPlane%2)*2.-1.);
	444	//
	445	clinfo->fNTracks =1;
	446	if (nearestcluster->GetLabel(1)>=0) clinfo->fNTracks++;
	447	if (nearestcluster->GetLabel(2)>=0) clinfo->fNTracks++;
	448	clinfo->Update();
	449	}
	450	//
	451	fTreeA->Fill();
	452	}
	453	}
	454
	455
	456	fFileA->cd();
	457	fTreeA->Write();
	458	fFileA->Close();
	459	return 0;
	460	}
59dfcadd	461
59dfcadd	462
289478c7	463	TH1F * AliTRDClusterErrDraw::ResDyVsAmp(TTree* tree, const char* selection, Float_t t0, Float_t ampmin, Float_t ampmax)
	464	{
	465	//
	466	//
	467	TH2F hisdy("resy","resy",10,ampmin,ampmax,30,-0.3,0.3);
	468	char expression[1000];
	469	sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fQ>>resy",t0);
	470	char selectionall[1000];
	471	sprintf(selectionall,"fStatus==0&&%s",selection);
	472	printf("%s\n",expression);
	473	printf("%s\n",selectionall);
	474	tree->Draw(expression,selectionall);
	475	return CreateResHisto(&hisdy);
	476	}
59dfcadd	477
59dfcadd	478
289478c7	479	TH1F * AliTRDClusterErrDraw::ResDyVsRelPos(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
	480	{
	481	//
	482	//
	483	min *=128;
	484	max *=128;
	485	TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
	486	char expression[1000];
	487	sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fRelPos>>resy",t0);
	488	char selectionall[1000];
	489	sprintf(selectionall,"fStatus==0&&%s",selection);
	490	printf("%s\n",expression);
	491	printf("%s\n",selectionall);
	492	tree->Draw(expression,selectionall);
	493	return CreateResHisto(&hisdy);
	494	}
59dfcadd	495
59dfcadd	496
289478c7	497	TH1F * AliTRDClusterErrDraw::ResDyVsAngleY(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
	498	{
	499	//
	500	//
	501	TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
59dfcadd	502
289478c7	503	char expression[1000];
	504	sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%f*fEp.fTAY:fEp.fTAY>>resy",t0);
	505	char selectionall[1000];
	506	sprintf(selectionall,"fStatus==0&&%s",selection);
59dfcadd	507
289478c7	508	tree->Draw(expression,selectionall);
	509	return CreateResHisto(&hisdy);
	510	}
59dfcadd	511
289478c7	512	void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
	513	{
	514	histo->GetXaxis()->SetTitle(xAxisTitle);
	515	histo->GetYaxis()->SetTitle(yAxisTitle);
	516	}
59dfcadd	517
59dfcadd	518
289478c7	519	TH1F* AliTRDClusterErrDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
	520	{
	521	Int_t nBins = hGen->GetNbinsX();
	522	TH1F* hEff = (TH1F*) hGen->Clone("hEff");
	523	hEff->SetTitle("");
	524	hEff->SetStats(kFALSE);
	525	hEff->SetMinimum(0.);
	526	hEff->SetMaximum(110.);
	527
	528	for (Int_t iBin = 0; iBin <= nBins; iBin++) {
	529	Double_t nGen = hGen->GetBinContent(iBin);
	530	Double_t nRec = hRec->GetBinContent(iBin);
	531	if (nGen > 0) {
	532	Double_t eff = nRec/nGen;
	533	hEff->SetBinContent(iBin, 100. * eff);
	534	Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
	535	// if (error == 0) error = sqrt(0.1/nGen);
	536	//
	537	if (error == 0) error = 0.0001;
	538	hEff->SetBinError(iBin, 100. * error);
	539	} else {
	540	hEff->SetBinContent(iBin, 100. * 0.5);
	541	hEff->SetBinError(iBin, 100. * 0.5);
	542	}
	543	}
59dfcadd	544
289478c7	545	return hEff;
289478c7	546	}
59dfcadd	547
59dfcadd	548
59dfcadd	549
289478c7	550	TH1F* AliTRDClusterErrDraw::CreateResHisto(TH2F* hRes2, Bool_t draw, Bool_t drawBinFits,
	551	Bool_t overflowBinFits)
	552	{
	553	TVirtualPad* currentPad = gPad;
	554	TAxis* axis = hRes2->GetXaxis();
	555	Int_t nBins = axis->GetNbins();
	556	TH1F* hRes, *hMean;
	557	if (axis->GetXbins()->GetSize()){
	558	hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
	559	hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
	560	}
	561	else{
	562	hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
	563	hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
59dfcadd	564
289478c7	565	}
	566	hRes->SetStats(false);
	567	hRes->SetOption("E");
	568	hRes->SetMinimum(0.);
	569	//
	570	hMean->SetStats(false);
	571	hMean->SetOption("E");
	572
	573	// create the fit function
	574	//TKFitGaus* fitFunc = new TKFitGaus("resFunc");
	575	TF1 * fitFunc = new TF1("G","[3]+[0]exp(-(x-[1])(x-[1])/(2.[2][2]))",-3,3);
	576
	577	fitFunc->SetLineWidth(2);
	578	fitFunc->SetFillStyle(0);
	579	// create canvas for fits
	580	TCanvas* canBinFits = NULL;
	581	Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
	582	Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
	583	Int_t ny = (nPads-1) / nx + 1;
	584	if (drawBinFits) {
	585	canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
	586	if (canBinFits) delete canBinFits;
	587	canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
	588	canBinFits->Divide(nx, ny);
	589	}
59dfcadd	590
289478c7	591	// loop over x bins and fit projection
	592	Int_t dBin = ((overflowBinFits) ? 1 : 0);
	593	for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
	594	if (drawBinFits) canBinFits->cd(bin + dBin);
	595	TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
	596	//
	597	if (hBin->GetEntries() > 10) {
	598	fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS(),0.02*hBin->GetMaximum());
	599	hBin->Fit(fitFunc,"s");
	600	Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
	601
	602	if (sigma > 0.){
	603	hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
	604	hMean->SetBinContent(bin, fitFunc->GetParameter(1));
59dfcadd	605	}
289478c7	606	else{
	607	hRes->SetBinContent(bin, 0.);
	608	hMean->SetBinContent(bin,0);
	609	}
	610	hRes->SetBinError(bin, fitFunc->GetParError(2));
	611	hMean->SetBinError(bin, fitFunc->GetParError(1));
	612
	613	//
	614	//
	615
	616	} else {
	617	hRes->SetBinContent(bin, 0.);
	618	hRes->SetBinError(bin, 0.);
	619	hMean->SetBinContent(bin, 0.);
	620	hMean->SetBinError(bin, 0.);
59dfcadd	621	}
289478c7	622
	623
	624	if (drawBinFits) {
	625	char name[256];
	626	if (bin == 0) {
	627	sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
	628	} else if (bin == nBins+1) {
	629	sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
	630	} else {
	631	sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
	632	axis->GetTitle(), axis->GetBinUpEdge(bin));
	633	}
	634	canBinFits->cd(bin + dBin);
	635	hBin->SetTitle(name);
	636	hBin->SetStats(kTRUE);
	637	hBin->DrawCopy("E");
	638	canBinFits->Update();
	639	canBinFits->Modified();
	640	canBinFits->Update();
	641	}
	642
	643	delete hBin;
	644	}
	645
	646	delete fitFunc;
	647	currentPad->cd();
	648	if (draw){
	649	currentPad->Divide(1,2);
	650	currentPad->cd(1);
	651	hRes->Draw();
	652	currentPad->cd(2);
	653	hMean->Draw();
59dfcadd	654	}
59dfcadd	655
289478c7	656	return hRes;
59dfcadd	657	}