PWGPP/ITS/AliAnalysisTaskdEdxSSDQA.cxx

   1 //SSD dEdx gain calibration task
   2 //to run in QA train
   3 //autor Marek Chojnacki
   4 //Marek.Chojnacki@cern.ch
   5 // last line
   6 #include "TChain.h"
   7 #include "TTree.h"
   8 #include "TH1F.h"
   9 #include "TH2F.h"
  10 #include "TH3F.h"
  11 #include "TFile.h"
  12
  13 #include "AliAnalysisTaskdEdxSSDQA.h"
  14 #include "AliAnalysisManager.h"
  15
  16 #include "AliESDEvent.h"
  17 #include "AliESDInputHandler.h"
  18
  19
  20 #include "AliESDtrack.h"
  21 #include "AliExternalTrackParam.h"
  22 #include "Riostream.h"
  23
  24 #include "AliESDfriend.h"
  25 #include "AliESDfriendTrack.h"
  26 #include "AliTrackPointArray.h"
  27 #include "AliCDBManager.h"
  28 #include "AliGeomManager.h"
  29 #include "TMath.h"
  30 #include "TGeoMatrix.h"
  31
  32 using namespace std;
  33 ClassImp(AliAnalysisTaskdEdxSSDQA)
  34
  35 //________________________________________________________________________
  36 AliAnalysisTaskdEdxSSDQA::AliAnalysisTaskdEdxSSDQA(const char *name)
  37 : AliAnalysisTaskSE(name), fHist1(0),fHist2(0),fHist3(0),fHist4(0),fHist5(0),fHist6(0),fListOfHistos(0),fPcut(0.0),fdothecorrection(0)
  38 {
  39   // Constructor
  40         // Define input and output slots here
  41   // Input slot #1 works with a TChain
  42  // fcorrections=new Float_t[1698*12];
  43   for(int i=0;i<1698*12;i++)
  44         fcorrections[i]=1.0;
  45   DefineOutput(1, TList::Class());
  46 }
  47
  48 //________________________________________________________________________
  49 void AliAnalysisTaskdEdxSSDQA::UserCreateOutputObjects()
  50 {
  51   // Create histograms
  52   // Called once
  53
  54         fListOfHistos = new TList();
  55         fListOfHistos->SetOwner();
  56         fHist1 =new TH2F("QAChargeRatio","QAChargeRatio;Module;CR",1698,-0.5,1697.5,80,-1.0,1.0);
  57         fListOfHistos->Add(fHist1);
  58         fHist2=new TH2F("QACharge","QACharge;Module;Q",1698,-0.5,1697.5,150,0,300);
  59         fListOfHistos->Add(fHist2);
  60         fHist3=new TH3F("ChargeRatiovCharge","CRvQ;module;Q;CR",1698,-0.5,1697.5,80,0,1600,50,-1.0,1.0);
  61         fListOfHistos->Add(fHist3);
  62         fHist4=new TH2F("QAChargedinchips","chargedinchips;chip;Q",1698*12,-0.5,1698*12-0.5,150,0,300);
  63         fListOfHistos->Add(fHist4);
  64         if(fdothecorrection)
  65         {
  66                 fHist5=new TH2F("QAChargeCorrected","chargedCorrected;Module;QCorrected",1698,-0.5,1697.5,150,0,300);
  67                 fListOfHistos->Add(fHist5);
  68         }
  69         fHist6=new TH2F("QNvQP",";QN;QP",160,0,1600,1600,0,1600);
  70         fListOfHistos->Add(fHist6);
  71
  72         PostData(1,  fListOfHistos);
  73 }
  74 //________________________________________________________________________
  75 void AliAnalysisTaskdEdxSSDQA::LocalInit()
  76 {
  77
  78         Printf("end of LocalInit");
  79 }
  80 //______________________________________________________________________
  81 AliAnalysisTaskdEdxSSDQA::~AliAnalysisTaskdEdxSSDQA()
  82 {
  83 }
  84
  85 //________________________________________________________________________
  86
  87 void AliAnalysisTaskdEdxSSDQA::UserExec(Option_t *)
  88 {
  89     // Main loop
  90     // Called for each event
  91     AliESDEvent* esd = dynamic_cast<AliESDEvent*> (InputEvent());
  92     if(!esd)
  93     {
  94         Printf("ERROR: Input ESD Event not available");
  95         PostData(1,  fListOfHistos);
  96         return;
  97     }
  98
  99     Int_t nTracks=esd->GetNumberOfTracks();
 100
 101
 102
 103         if(!ESDfriend())
 104         {
 105                 Printf("problem with friend");
 106                 PostData(1,  fListOfHistos);
 107                 return;
 108         }
 109         //Printf("Event nTracks %d", nTracks);
 110         AliTrackPointArray*  trackar=0x0;
 111         Bool_t l5;
 112         Bool_t l6;
 113         Int_t npoints;
 114         AliTrackPoint point;
 115         Int_t nlayer=0;
 116         Int_t id=0;
 117         Float_t chargeratio=0.0;
 118         Float_t charge=0.0;
 119         for(int itr=0;itr<nTracks;itr++)
 120         {
 121                 AliESDtrack* track= esd->GetTrack(itr);
 122
 123                 if(TMath::Abs(track->Eta())>0.9)
 124                         continue;
 125                 if(track->GetP()>10.0)
 126                         continue;
 127                 if (track->IsOn(AliESDtrack::kITSrefit)&&track->IsOn(AliESDtrack::kTPCrefit))
 128                 {
 129                          l5=track->HasPointOnITSLayer(4);
 130                          l6=track->HasPointOnITSLayer(5);
 131                         if (!(l5||l6))//only tracks with SSD point
 132                                 continue;
 133                         Double_t tmpQESD[4]={-1.0,-1.0,-1.0,-1.0};
 134                         track->GetITSdEdxSamples(tmpQESD);
 135                         trackar=(AliTrackPointArray*)track->GetTrackPointArray();
 136                         if(!trackar)
 137                                 continue;
 138                         npoints=trackar->GetNPoints();
 139                         for(int itnp=0;itnp<npoints;itnp++)
 140                         {
 141                                 if(trackar->GetPoint(point,itnp))
 142                                 {
 143                                         nlayer=AliGeomManager::VolUIDToLayer(point.GetVolumeID(), id);//layer number
 144                                 }
 145                                 else
 146                                         continue;
 147                                 if(nlayer==5||nlayer==6)
 148                                 {
 149                                         TGeoHMatrix* geomatrix=AliGeomManager::GetMatrix(point.GetVolumeID());
 150                                         //geomatrix->Print();
 151                                         if(point.GetCharge()>0.0&&point.IsExtra()==kFALSE)//do not use additional clusters
 152                                         {
 153                                                 Double_t local[3]={0.0,0.0,0.0};
 154                                                 Double_t global[3]={0.0,0.0,0.0};
 155                                                 global[0]=point.GetX();
 156                                                 global[1]=point.GetY();
 157                                                 global[2]=point.GetZ();
 158                                                 geomatrix->MasterToLocal(global,local);
 159                                                  chargeratio=point.GetChargeRatio();
 160                                                  charge=point.GetCharge();
 161                                                  Float_t fQNnotcorr=charge*(1.0+chargeratio);
 162                                                  Float_t fQPnotcorr=charge*(1.0-chargeratio);
 163                                                  fHist6->Fill(fQNnotcorr,fQPnotcorr);
 164                                                 // cout<<point.GetCharge()<<" "<<point.GetChargeRatio()<<" "<<nlayer<<" "<<id<<endl;
 165                                                 if(nlayer==5&&tmpQESD[2]>0.0)
 166                                                 {
 167                                                         fHist1->Fill(id,chargeratio);
 168                                                         if(track->GetP()>fPcut)
 169                                                         {
 170                                                                 Float_t fQP=tmpQESD[2]*(1.0-chargeratio);
 171                                                                 Float_t fQN=tmpQESD[2]*(1.0+chargeratio);
 172                                                                 Int_t fPchip=Pstrip5(10000.0*local[0],10000.0*local[2])/128;
 173                                                                 Int_t fNchip=Nstrip5(10000.0*local[0],10000.0*local[2])/128;
 174
 175                                                                 fHist2->Fill(id,tmpQESD[2]);
 176                                                                 fHist4->Fill(id*12+fPchip,fQP);
 177                                                                 fHist4->Fill(id*12+6+fNchip,fQN);
 178                                                                 if(fdothecorrection)
 179                                                                         fHist5->Fill(id,0.5*(fQP*fcorrections[12*id+fPchip]+fQN*fcorrections[12*id+6+fNchip]));
 180                                                         }
 181                                                         fHist3->Fill(id,charge,chargeratio);
 182                                                 }
 183                                                 if(nlayer==6&&tmpQESD[3]>0.0)
 184                                                 {
 185                                                         fHist1->Fill(id+748,chargeratio);
 186                                                         if(track->GetP()>fPcut)
 187                                                         {
 188
 189                                                                 Float_t fQP=tmpQESD[3]*(1.0-chargeratio);
 190                                                                 Float_t fQN=tmpQESD[3]*(1.0+chargeratio);
 191                                                                 Int_t fPchip=Pstrip6(10000.0*local[0],10000.0*local[2])/128;
 192                                                                 Int_t fNchip=Nstrip6(10000.0*local[0],10000.0*local[2])/128;
 193
 194                                                                 fHist2->Fill(id+748,tmpQESD[3]);
 195                                                                 fHist4->Fill((id+748)*12+fPchip,fQP);
 196                                                                 fHist4->Fill((id+748)*12+6+fNchip,fQN);
 197                                                                 if(fdothecorrection)
 198                                                                         fHist5->Fill((id+748),0.5*(fQP*fcorrections[12*(id+748)+fPchip]+fQN*fcorrections[12*(id+748)+6+fNchip]));
 199
 200                                                         }
 201                                                         fHist3->Fill(id+748,charge,chargeratio);
 202                                                 }
 203
 204                                         }
 205                                 }
 206                         }
 207                 }
 208         }
 209         //Printf("Event nTracks end %d", nTracks);
 210         // Post output data.
 211         PostData(1,  fListOfHistos);
 212 }
 213
 214 //________________________________________________________________________
 215 void AliAnalysisTaskdEdxSSDQA::Terminate(Option_t *)
 216 {
 217         //terminate
 218         Printf("end of Terminate");
 219 }
 220 //_____________________________________________________________________________
 221 Int_t  AliAnalysisTaskdEdxSSDQA::Pstrip5(Float_t x,Float_t z) const
 222 {
 223         // P strip for layer 5 from local cooridnates
 224         Float_t value;
 225         x=x+36432.5-223.9;
 226         z=z-(76000.0/7.0);
 227         value=-(3.0/38000.0)*z+x/95.0;
 228         return TMath::Nint(value);
 229 }
 230 //___________________________________________________________________________
 231 Int_t  AliAnalysisTaskdEdxSSDQA::Nstrip5(Float_t x,Float_t z)  const
 232 {
 233         // N strip for layer 5 from local cooridnates
 234         Float_t value;
 235         x=x+36432.5-223.9;
 236         z=z-(76000.0/7.0);
 237         value=x/95.0+(11.0/38000.0)*z;
 238         return TMath::Nint(value);
 239 }
 240 //___________________________________________________________________________
 241 Int_t  AliAnalysisTaskdEdxSSDQA::Pstrip6(Float_t x,Float_t z)  const
 242 {
 243         // P strip for layer 6 from local cooridnates
 244         Float_t value;
 245         x=x-36432.5-223.9;
 246         z=z+(76000.0/7.0);
 247         value=-(3.0/38000.0)*z-x/95.0;
 248         return TMath::Nint(value);
 249 }
 250 //__________________________________________________________________________________
 251 Int_t  AliAnalysisTaskdEdxSSDQA::Nstrip6(Float_t x,Float_t z)  const
 252 {
 253         // N strip for layer 6 from local cooridnates
 254         Float_t value;
 255         x=x-36432.5-223.9;
 256         z=z+(76000.0/7.0);
 257         value=x*(-1.0/95.0)+(11.0/38000.0)*z;
 258         return TMath::Nint(value);
 259 }
 260 //___________________________________________________________________________________
 261 void  AliAnalysisTaskdEdxSSDQA::SetDoChipCorretions(const char* filename)
 262 {
 263 //Upload corrections for each chip only for test
 264         cout<<filename<<endl;
 265         ifstream infile(filename);
 266         if(!infile.good())
 267                 return;
 268         Float_t value=0.0;
 269         fdothecorrection=1;
 270         for (int i=0;i<1698;i++)
 271         {
 272                 for (int j=0;j<6;j++)
 273                 {
 274                         infile>>value;
 275                         cout<<value<<" ";
 276                         if(value>0.0)
 277                                 fcorrections[i*12+j]=84.0/value;
 278                         else
 279                                 fcorrections[i*12+j]=1.0;
 280                         infile>>value;
 281                         infile>>value;
 282                         cout<<value<<" ";
 283                         if(value>0.0)
 284                                 fcorrections[i*12+j+6]=84.0/value;
 285                         else
 286                                 fcorrections[i*12+j+6]=1.0;
 287                         infile>>value;
 288                 }
 289                 cout<<endl;
 290         }
 291 }
 292