[u/mrichter/AliRoot.git] / TPC / Upgrade / macros / spaceChargeFluctuation.C

/*
  .x ~/rootlogon.C
  .L $ALICE_ROOT/TPC/Upgrade/macros/spaceChargeFluctuation.C+ 
//  GenerateMapRaw("/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 50);
  GenerateMapRawIons( "/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 25);
  
 */
#include "TMath.h"
#include "TRandom.h"
#include "TTreeStream.h"
#include "TVectorD.h"
#include "TCanvas.h"
#include "TStopwatch.h"
#include "AliTPCParam.h"
#include "AliTPCcalibDB.h"
#include "AliTPCAltroMapping.h"
#include "AliAltroRawStream.h"
#include "AliSysInfo.h"
#include "AliTPCRawStreamV3.h"
#include "AliCDBManager.h"
#include "TGeoGlobalMagField.h"
#include "AliMagF.h"
#include "AliRawReaderRoot.h"
#include "AliRawReader.h"
#include "TH3.h"
#include "TH2.h"


TH1 * GenerateMapRawIons(const char *fileName="raw.root", const char *outputName="histo.root", Int_t maxEvents=25);


void spaceChargeFluctuation(Int_t mode=0){
  //
  //
  if (mode==0) GenerateMapRawIons();
  
}

void spaceChargeFluctuationMC(Int_t nframes){
  //
  // Toy MC to generate space charge fluctuation
  //
  //
  TTreeSRedirector *pcstream = new TTreeSRedirector("spaceChargeFluctuation.root","recreate");
  Double_t interactionRate=100000;
  Double_t driftTime=0.1;

  for (Int_t iframe=0; iframe<nframes; iframe++){
    Int_t nevents=gRandom->Poisson(interactionRate*driftTime);
    Int_t ntracksAll=0;
    TVectorD vecAllPhi128(128);
    TVectorD vecAllPhi36(36);
    for (Int_t ievent=0; ievent<nevents; ievent++){
      //
      //
      Int_t ntracks=gRandom->Exp(500)+gRandom->Exp(500); // to be taken form the distribution          
      ntracksAll+=ntracks; 
      for (Int_t isec=0; isec<128; isec++){
	vecAllPhi128(isec)+=ntracks*gRandom->Rndm()/128;
      }
      for (Int_t isec=0; isec<36; isec++){
	vecAllPhi36(isec)+=ntracks*gRandom->Rndm()/36;
      }
    }
    (*pcstream)<<"ntracks"<<
      "nevents="<<nevents<<                  // number of events withing time frame
      "ntracksAll="<<ntracksAll<<            // number of tracks within  time frame
      "vecAllPhi36="<<&vecAllPhi36<<         // number of tracks in phi bin (36 bins)    within time frame
      "vecAllPhi128="<<&vecAllPhi128<<       // number of tracks in phi bin (128 bins)   within time frame
      "\n";
  }
  delete pcstream;
}


// void spaceChargeFluctuationDraw(){
//   //
//   //
//   //
//   TFile *f = TFile::Open("spaceChargeFluctuation.root");
//   TTree * tree = (TTree*)f->Get("ntracks");
//   TCanvas * canvasFluc = new TCanvas("cavasFluc","canvasFluc");  
// }


TH1 * GenerateMapRaw(const char *fileName, const char *outputName, Int_t maxEvents){
  //
  // Generate 3D maps of electrons 
  // different threshold considered
  // Paramaters:
  //    fileName      - name of input raw file
  //    outputName    - name of output file with the space charge histograms 
  //    maxEvents     - grouping of the events
  // 
  //  
  TTreeSRedirector * pcstream  = new TTreeSRedirector(outputName, "recreate");

  const char *  ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";  
  AliCDBManager * man = AliCDBManager::Instance();
  man->SetDefaultStorage(ocdbpath);
  man->SetRun(0);
  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG,       AliMagF::kBeamTypepp, 2.76/2.));
  AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
  AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
  //
  TStopwatch timer;
  timer.Start();
  TH3F * hisQ3D[3]={0};
  TH2F * hisQ2D[3]={0};
  TH1F * hisQ1D[3]={0};
  //
  for (Int_t ith=0; ith<3; ith++){
    char chname[100];
    snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
    hisQ3D[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,100,0,1100);
    snprintf(chname,100,"hisQ2D_Th%d",2*ith+2);
    hisQ2D[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
    //
    snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
    hisQ1D[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
    hisQ3D[ith]->SetDirectory(0);
    hisQ1D[ith]->SetDirectory(0);
    hisQ1D[ith]->SetDirectory(0);
  }
  //
  //  
  AliRawReader *reader = new AliRawReaderRoot(fileName);
  reader->Reset();
  AliAltroRawStream* stream = new AliAltroRawStream(reader);
  stream->SelectRawData("TPC");
  Int_t evtnr=0;
  Int_t chunkNr=0;
  //
  while (reader->NextEvent()) {
    if(evtnr>maxEvents) {
      chunkNr++;
      evtnr=0;
      pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
      pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
      for (Int_t ith=0; ith<3; ith++){	
	hisQ1D[ith]->Write(Form("His1D_%d",ith));
	hisQ2D[ith]->Write(Form("His2D_%d",ith));
	hisQ3D[ith]->Write(Form("His3D_%d",ith));
	if (ith<3) continue;
	hisQ1D[ith]->Reset();
	hisQ2D[ith]->Reset();
	hisQ3D[ith]->Reset();
      }
    }
    cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
    evtnr++;
    AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
    AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
    //
    while (input.NextDDL()){
      Int_t sector = input.GetSector(); 
      while ( input.NextChannel() ) {
	Int_t    row    = input.GetRow();
	Int_t rowAbs    = row+ ((sector<36) ? 0: param->GetNRow(0));
	Int_t    pad    = input.GetPad();
	Int_t    nPads   = param->GetNPads(sector,row);
	Double_t localX  = param->GetPadRowRadii(sector,row); 
	Double_t localY  = (pad-nPads)*param->GetPadPitchWidth(sector);
	Double_t localPhi= TMath::ATan2(localY,localX);
	Double_t phi     = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
	if (sector%36>=18) phi+=TMath::TwoPi();
	Double_t padLength=param->GetPadPitchLength(sector,row);
	//      Double_t padWidth=param->GetPadPitchWidth(sector);
	//    Double_t zLength=param->GetZLength();
	Double_t deltaPhi=hisQ3D[0]->GetXaxis()->GetBinWidth(0);
        Double_t deltaZ= hisQ3D[0]->GetZaxis()->GetBinWidth(0);
	//
	while ( input.NextBunch() ){
	  Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
	  Int_t  bunchlength  = (Int_t)input.GetBunchLength();
	  const UShort_t *sig = input.GetSignals();	  
	  Int_t aboveTh[3]={0};
	  for (Int_t i=0; i<bunchlength; i++){ 
	    for (Int_t ith=0; ith<3; ith++){
	      if (sig[i]>(ith*2)+2) aboveTh[ith]++; 
	    }
	  }
	  for (Int_t ith=0; ith<3; ith++){
	    if (aboveTh[ith%3]>1){
	      for (Int_t i=0; i<bunchlength; i++){
		//
		// normalization
		//
		Double_t volume3D = padLength*(localX*deltaPhi)*deltaZ;
		if (sector%36<18) hisQ1D[ith]->Fill(rowAbs, sig[i]/volume3D);
		hisQ2D[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
		hisQ3D[ith]->Fill(phi,rowAbs,startTbin+i,sig[i]/volume3D);
	      }
	    }
	  }
	}
      }
    }
  }
  timer.Print();
  delete pcstream;
  return 0;
}


TH1 * GenerateMapRawIons(const char *fileName, const char *outputName, Int_t maxEvents){
  //
  // Generate 3D maps of the space charge for the rad data maps
  // different threshold considered
  // Paramaters:
  //    fileName      - name of input raw file
  //    outputName    - name of output file with the space charge histograms 
  //    maxEvents     - grouping of the events
  // 
  //  
  TTreeSRedirector * pcstream  = new TTreeSRedirector(outputName, "recreate");
  const char *  ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";  
  AliCDBManager * man = AliCDBManager::Instance();
  man->SetDefaultStorage(ocdbpath);
  man->SetRun(0);
  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG,       AliMagF::kBeamTypepp, 2.76/2.));
  AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
  AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
  TStopwatch timer;
  timer.Start();
  //   arrays of space charges - different elements corresponds to different threshold to accumulate charge
  TH3F * hisQ3D[3]={0};                // 3D maps space charge from drift volume
  TH2F * hisQ2D[3]={0};                
  TH1F * hisQ1D[3]={0};
  TH3F * hisQ3DROC[3]={0};             // 3D maps space charge from ROC
  TH2F * hisQ2DROC[3]={0};
  TH1F * hisQ1DROC[3]={0};
  //
  for (Int_t ith=0; ith<3; ith++){
    char chname[100];
    snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
    hisQ3D[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
    snprintf(chname,100,"hisQ2D_Th%d",2*ith+2);
    hisQ2D[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
    snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
    hisQ1D[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
    //
    snprintf(chname,100,"hisQ3DROC_Th%d",2*ith+2);
    hisQ3DROC[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
    snprintf(chname,100,"hisQ2DROC_Th%d",2*ith+2);
    hisQ2DROC[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
    snprintf(chname,100,"hisQ1DROC_Th%d",2*ith+2);
    hisQ1DROC[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
    //
    hisQ3D[ith]->SetDirectory(0);
    hisQ1D[ith]->SetDirectory(0);
    hisQ1D[ith]->SetDirectory(0);
    hisQ3DROC[ith]->SetDirectory(0);
    hisQ1DROC[ith]->SetDirectory(0);
    hisQ1DROC[ith]->SetDirectory(0);
  }
  //
  //  
  AliRawReader *reader = new AliRawReaderRoot(fileName);
  reader->Reset();
  AliAltroRawStream* stream = new AliAltroRawStream(reader);
  stream->SelectRawData("TPC");
  Int_t evtnr=0;
  Int_t chunkNr=0;
  // 

  while (reader->NextEvent()) {
    Double_t shiftZ= gRandom->Rndm()*250;
    //
    if(evtnr>=maxEvents) {
      chunkNr++;
      evtnr=0;
      pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
      pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
      for (Int_t ith=0; ith<3; ith++){	
	hisQ1D[ith]->Write(Form("His1DDrift_%d",ith));
	hisQ2D[ith]->Write(Form("His2DDrift_%d",ith));
	hisQ3D[ith]->Write(Form("His3DDrift_%d",ith));
	hisQ1DROC[ith]->Write(Form("His1DROC_%d",ith));
	hisQ2DROC[ith]->Write(Form("His2DROC_%d",ith));
	hisQ3DROC[ith]->Write(Form("His3DROC_%d",ith));
	if (ith<3) continue;
	hisQ1D[ith]->Reset();
	hisQ2D[ith]->Reset();
	hisQ3D[ith]->Reset();
	hisQ1DROC[ith]->Reset();
	hisQ2DROC[ith]->Reset();
	hisQ3DROC[ith]->Reset();
      }
    }
    cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
    evtnr++;
    AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
    AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
    //
    while (input.NextDDL()){
      Int_t sector = input.GetSector(); 
      while ( input.NextChannel() ) {
	Int_t    row    = input.GetRow();
	Int_t rowAbs    = row+ ((sector<36) ? 0: param->GetNRow(0));
	Int_t    pad    = input.GetPad();
	Int_t    nPads   = param->GetNPads(sector,row);
	Double_t localX  = param->GetPadRowRadii(sector,row); 
	Double_t localY  = (pad-nPads)*param->GetPadPitchWidth(sector);
	Double_t localPhi= TMath::ATan2(localY,localX);
	Double_t phi     = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
	if (sector%36>=18) phi+=TMath::TwoPi();
	Double_t padLength=param->GetPadPitchLength(sector,row);
	//  Double_t padWidth=param->GetPadPitchWidth(sector);
	//  Double_t zLength=param->GetZLength();
	Double_t deltaPhi=hisQ3D[0]->GetXaxis()->GetBinWidth(0);
        Double_t deltaZ= hisQ3D[0]->GetZaxis()->GetBinWidth(0);
	//
	while ( input.NextBunch() ){
	  Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
	  Int_t  bunchlength  = (Int_t)input.GetBunchLength();
	  const UShort_t *sig = input.GetSignals();	  
	  Int_t aboveTh[3]={0};
	  for (Int_t i=0; i<bunchlength; i++){ 
	    for (Int_t ith=0; ith<3; ith++){
	      if (sig[i]>(ith*2)+2) aboveTh[ith]++; 
	    }
	  }
	  for (Int_t ith=0; ith<3; ith++){
	    if (aboveTh[ith%3]>1){
	      for (Int_t i=0; i<bunchlength; i++){
		//
		// normalization
		//
		Double_t volume3D  = padLength*(localX*deltaPhi)*deltaZ;
		Double_t zIonDrift = (param->GetZLength()-startTbin*param->GetZWidth());
		zIonDrift+=shiftZ;
                if (TMath::Abs(zIonDrift)<param->GetZLength()){
		  if ((sector%36)>=18) zIonDrift*=-1;   // c side has opposite sign
		  if (sector%36<18) hisQ1D[ith]->Fill(rowAbs, sig[i]/volume3D);
		  hisQ2D[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
		  hisQ3D[ith]->Fill(phi,rowAbs,zIonDrift,sig[i]/volume3D);
		}
		//
		Double_t zIonROC = ((sector%36)<18)? shiftZ: -shiftZ;  // z position of the "ion disc" -  A side C side opposite sign
		if (sector%36<18) hisQ1DROC[ith]->Fill(rowAbs, sig[i]/volume3D);
		hisQ2DROC[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
		hisQ3DROC[ith]->Fill(phi,rowAbs,zIonROC,sig[i]/volume3D);
	      }
	    }
	  }
	}
      }
    }
  }
  timer.Print();
  delete pcstream;
  return 0;
}
Commit	Line	Data
c0172f82	1	/*
	2	.x ~/rootlogon.C
	3	.L $ALICE_ROOT/TPC/Upgrade/macros/spaceChargeFluctuation.C+
	4	// GenerateMapRaw("/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 50);
	5	GenerateMapRawIons( "/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 25);
	6
	7	*/
	8	#include "TMath.h"
	9	#include "TRandom.h"
	10	#include "TTreeStream.h"
	11	#include "TVectorD.h"
	12	#include "TCanvas.h"
	13	#include "TStopwatch.h"
	14	#include "AliTPCParam.h"
	15	#include "AliTPCcalibDB.h"
	16	#include "AliTPCAltroMapping.h"
	17	#include "AliAltroRawStream.h"
	18	#include "AliSysInfo.h"
	19	#include "AliTPCRawStreamV3.h"
	20	#include "AliCDBManager.h"
	21	#include "TGeoGlobalMagField.h"
	22	#include "AliMagF.h"
	23	#include "AliRawReaderRoot.h"
	24	#include "AliRawReader.h"
	25	#include "TH3.h"
	26	#include "TH2.h"
	27
	28
	29	TH1 * GenerateMapRawIons(const char fileName="raw.root", const char outputName="histo.root", Int_t maxEvents=25);
	30
	31
	32	void spaceChargeFluctuation(Int_t mode=0){
	33	//
	34	//
	35	if (mode==0) GenerateMapRawIons();
	36
	37	}
	38
	39	void spaceChargeFluctuationMC(Int_t nframes){
	40	//
	41	// Toy MC to generate space charge fluctuation
	42	//
	43	//
	44	TTreeSRedirector *pcstream = new TTreeSRedirector("spaceChargeFluctuation.root","recreate");
	45	Double_t interactionRate=100000;
	46	Double_t driftTime=0.1;
	47
	48	for (Int_t iframe=0; iframe<nframes; iframe++){
	49	Int_t nevents=gRandom->Poisson(interactionRate*driftTime);
	50	Int_t ntracksAll=0;
	51	TVectorD vecAllPhi128(128);
	52	TVectorD vecAllPhi36(36);
	53	for (Int_t ievent=0; ievent<nevents; ievent++){
	54	//
	55	//
	56	Int_t ntracks=gRandom->Exp(500)+gRandom->Exp(500); // to be taken form the distribution
	57	ntracksAll+=ntracks;
	58	for (Int_t isec=0; isec<128; isec++){
	59	vecAllPhi128(isec)+=ntracks*gRandom->Rndm()/128;
	60	}
	61	for (Int_t isec=0; isec<36; isec++){
	62	vecAllPhi36(isec)+=ntracks*gRandom->Rndm()/36;
	63	}
	64	}
65	(*pcstream)<<"ntracks"<<
66	"nevents="<<nevents<< // number of events withing time frame
67	"ntracksAll="<<ntracksAll<< // number of tracks within time frame
68	"vecAllPhi36="<<&vecAllPhi36<< // number of tracks in phi bin (36 bins) within time frame
69	"vecAllPhi128="<<&vecAllPhi128<< // number of tracks in phi bin (128 bins) within time frame
70	"\n";
71	}
72	delete pcstream;
73	}
74
75
76	// void spaceChargeFluctuationDraw(){
77	// //
78	// //
79	// //
80	// TFile *f = TFile::Open("spaceChargeFluctuation.root");
81	// TTree * tree = (TTree*)f->Get("ntracks");
82	// TCanvas * canvasFluc = new TCanvas("cavasFluc","canvasFluc");
83	// }
84
85
86	TH1 * GenerateMapRaw(const char fileName, const char outputName, Int_t maxEvents){
87	//
88	// Generate 3D maps of electrons
89	// different threshold considered
90	// Paramaters:
91	// fileName - name of input raw file
92	// outputName - name of output file with the space charge histograms
93	// maxEvents - grouping of the events
94	//
95	//
96	TTreeSRedirector * pcstream = new TTreeSRedirector(outputName, "recreate");
97
98	const char * ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";
99	AliCDBManager * man = AliCDBManager::Instance();
100	man->SetDefaultStorage(ocdbpath);
101	man->SetRun(0);
102	TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 2.76/2.));
103	AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
104	AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
105	//
106	TStopwatch timer;
107	timer.Start();
108	TH3F * hisQ3D[3]={0};
109	TH2F * hisQ2D[3]={0};
110	TH1F * hisQ1D[3]={0};
111	//
112	for (Int_t ith=0; ith<3; ith++){
113	char chname[100];
114	snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
115	hisQ3D[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,100,0,1100);
116	snprintf(chname,100,"hisQ2D_Th%d",2*ith+2);
117	hisQ2D[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
118	//
119	snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
120	hisQ1D[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
121	hisQ3D[ith]->SetDirectory(0);
122	hisQ1D[ith]->SetDirectory(0);
123	hisQ1D[ith]->SetDirectory(0);
124	}
125	//
126	//
127	AliRawReader *reader = new AliRawReaderRoot(fileName);
128	reader->Reset();
129	AliAltroRawStream* stream = new AliAltroRawStream(reader);
130	stream->SelectRawData("TPC");
131	Int_t evtnr=0;
132	Int_t chunkNr=0;
133	//
134	while (reader->NextEvent()) {
135	if(evtnr>maxEvents) {
136	chunkNr++;
137	evtnr=0;
138	pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
139	pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
140	for (Int_t ith=0; ith<3; ith++){
141	hisQ1D[ith]->Write(Form("His1D_%d",ith));
142	hisQ2D[ith]->Write(Form("His2D_%d",ith));
143	hisQ3D[ith]->Write(Form("His3D_%d",ith));
144	if (ith<3) continue;
145	hisQ1D[ith]->Reset();
146	hisQ2D[ith]->Reset();
147	hisQ3D[ith]->Reset();
148	}
149	}
150	cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
151	evtnr++;
152	AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
153	AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
154	//
155	while (input.NextDDL()){
156	Int_t sector = input.GetSector();
157	while ( input.NextChannel() ) {
158	Int_t row = input.GetRow();
159	Int_t rowAbs = row+ ((sector<36) ? 0: param->GetNRow(0));
160	Int_t pad = input.GetPad();
161	Int_t nPads = param->GetNPads(sector,row);
162	Double_t localX = param->GetPadRowRadii(sector,row);
163	Double_t localY = (pad-nPads)*param->GetPadPitchWidth(sector);
164	Double_t localPhi= TMath::ATan2(localY,localX);
165	Double_t phi = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
166	if (sector%36>=18) phi+=TMath::TwoPi();
167	Double_t padLength=param->GetPadPitchLength(sector,row);
168	// Double_t padWidth=param->GetPadPitchWidth(sector);
169	// Double_t zLength=param->GetZLength();
170	Double_t deltaPhi=hisQ3D[0]->GetXaxis()->GetBinWidth(0);
171	Double_t deltaZ= hisQ3D[0]->GetZaxis()->GetBinWidth(0);
172	//
173	while ( input.NextBunch() ){
174	Int_t startTbin = (Int_t)input.GetStartTimeBin();
175	Int_t bunchlength = (Int_t)input.GetBunchLength();
176	const UShort_t *sig = input.GetSignals();
177	Int_t aboveTh[3]={0};
178	for (Int_t i=0; i<bunchlength; i++){
179	for (Int_t ith=0; ith<3; ith++){
180	if (sig[i]>(ith*2)+2) aboveTh[ith]++;
181	}
182	}
183	for (Int_t ith=0; ith<3; ith++){
184	if (aboveTh[ith%3]>1){
185	for (Int_t i=0; i<bunchlength; i++){
186	//
187	// normalization
188	//
189	Double_t volume3D = padLength(localXdeltaPhi)*deltaZ;
190	if (sector%36<18) hisQ1D[ith]->Fill(rowAbs, sig[i]/volume3D);
191	hisQ2D[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
192	hisQ3D[ith]->Fill(phi,rowAbs,startTbin+i,sig[i]/volume3D);
193	}
194	}
195	}
196	}
197	}
198	}
199	}
200	timer.Print();
201	delete pcstream;
202	return 0;
203	}
204
205
206
207
208
209	TH1 * GenerateMapRawIons(const char fileName, const char outputName, Int_t maxEvents){
210	//
211	// Generate 3D maps of the space charge for the rad data maps
212	// different threshold considered
213	// Paramaters:
214	// fileName - name of input raw file
215	// outputName - name of output file with the space charge histograms
216	// maxEvents - grouping of the events
217	//
218	//
219	TTreeSRedirector * pcstream = new TTreeSRedirector(outputName, "recreate");
220	const char * ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";
221	AliCDBManager * man = AliCDBManager::Instance();
222	man->SetDefaultStorage(ocdbpath);
223	man->SetRun(0);
224	TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 2.76/2.));
225	AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
226	AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
227	TStopwatch timer;
228	timer.Start();
229	// arrays of space charges - different elements corresponds to different threshold to accumulate charge
230	TH3F * hisQ3D[3]={0}; // 3D maps space charge from drift volume
231	TH2F * hisQ2D[3]={0};
232	TH1F * hisQ1D[3]={0};
233	TH3F * hisQ3DROC[3]={0}; // 3D maps space charge from ROC
234	TH2F * hisQ2DROC[3]={0};
235	TH1F * hisQ1DROC[3]={0};
236	//
237	for (Int_t ith=0; ith<3; ith++){
238	char chname[100];
239	snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
240	hisQ3D[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
241	snprintf(chname,100,"hisQ2D_Th%d",2*ith+2);
242	hisQ2D[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
243	snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
244	hisQ1D[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
245	//
246	snprintf(chname,100,"hisQ3DROC_Th%d",2*ith+2);
247	hisQ3DROC[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
248	snprintf(chname,100,"hisQ2DROC_Th%d",2*ith+2);
249	hisQ2DROC[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
250	snprintf(chname,100,"hisQ1DROC_Th%d",2*ith+2);
251	hisQ1DROC[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
252	//
253	hisQ3D[ith]->SetDirectory(0);
254	hisQ1D[ith]->SetDirectory(0);
255	hisQ1D[ith]->SetDirectory(0);
256	hisQ3DROC[ith]->SetDirectory(0);
257	hisQ1DROC[ith]->SetDirectory(0);
258	hisQ1DROC[ith]->SetDirectory(0);
259	}
260	//
261	//
262	AliRawReader *reader = new AliRawReaderRoot(fileName);
263	reader->Reset();
264	AliAltroRawStream* stream = new AliAltroRawStream(reader);
265	stream->SelectRawData("TPC");
266	Int_t evtnr=0;
267	Int_t chunkNr=0;
268	//
269
270	while (reader->NextEvent()) {
271	Double_t shiftZ= gRandom->Rndm()*250;
272	//
273	if(evtnr>=maxEvents) {
274	chunkNr++;
275	evtnr=0;
276	pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
277	pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
278	for (Int_t ith=0; ith<3; ith++){
279	hisQ1D[ith]->Write(Form("His1DDrift_%d",ith));
280	hisQ2D[ith]->Write(Form("His2DDrift_%d",ith));
281	hisQ3D[ith]->Write(Form("His3DDrift_%d",ith));
282	hisQ1DROC[ith]->Write(Form("His1DROC_%d",ith));
283	hisQ2DROC[ith]->Write(Form("His2DROC_%d",ith));
284	hisQ3DROC[ith]->Write(Form("His3DROC_%d",ith));
285	if (ith<3) continue;
286	hisQ1D[ith]->Reset();
287	hisQ2D[ith]->Reset();
288	hisQ3D[ith]->Reset();
289	hisQ1DROC[ith]->Reset();
290	hisQ2DROC[ith]->Reset();
291	hisQ3DROC[ith]->Reset();
292	}
293	}
294	cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
295	evtnr++;
296	AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
297	AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
298	//
299	while (input.NextDDL()){
300	Int_t sector = input.GetSector();
301	while ( input.NextChannel() ) {
302	Int_t row = input.GetRow();
303	Int_t rowAbs = row+ ((sector<36) ? 0: param->GetNRow(0));
304	Int_t pad = input.GetPad();
305	Int_t nPads = param->GetNPads(sector,row);
306	Double_t localX = param->GetPadRowRadii(sector,row);
307	Double_t localY = (pad-nPads)*param->GetPadPitchWidth(sector);
308	Double_t localPhi= TMath::ATan2(localY,localX);
309	Double_t phi = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
310	if (sector%36>=18) phi+=TMath::TwoPi();
311	Double_t padLength=param->GetPadPitchLength(sector,row);
312	// Double_t padWidth=param->GetPadPitchWidth(sector);
313	// Double_t zLength=param->GetZLength();
314	Double_t deltaPhi=hisQ3D[0]->GetXaxis()->GetBinWidth(0);
315	Double_t deltaZ= hisQ3D[0]->GetZaxis()->GetBinWidth(0);
316	//
317	while ( input.NextBunch() ){
318	Int_t startTbin = (Int_t)input.GetStartTimeBin();
319	Int_t bunchlength = (Int_t)input.GetBunchLength();
320	const UShort_t *sig = input.GetSignals();
321	Int_t aboveTh[3]={0};
322	for (Int_t i=0; i<bunchlength; i++){
323	for (Int_t ith=0; ith<3; ith++){
324	if (sig[i]>(ith*2)+2) aboveTh[ith]++;
325	}
326	}
327	for (Int_t ith=0; ith<3; ith++){
328	if (aboveTh[ith%3]>1){
329	for (Int_t i=0; i<bunchlength; i++){
330	//
331	// normalization
332	//
333	Double_t volume3D = padLength(localXdeltaPhi)*deltaZ;
334	Double_t zIonDrift = (param->GetZLength()-startTbin*param->GetZWidth());
335	zIonDrift+=shiftZ;
336	if (TMath::Abs(zIonDrift)<param->GetZLength()){
337	if ((sector%36)>=18) zIonDrift*=-1; // c side has opposite sign
338	if (sector%36<18) hisQ1D[ith]->Fill(rowAbs, sig[i]/volume3D);
339	hisQ2D[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
340	hisQ3D[ith]->Fill(phi,rowAbs,zIonDrift,sig[i]/volume3D);
341	}
342	//
343	Double_t zIonROC = ((sector%36)<18)? shiftZ: -shiftZ; // z position of the "ion disc" - A side C side opposite sign
344	if (sector%36<18) hisQ1DROC[ith]->Fill(rowAbs, sig[i]/volume3D);
345	hisQ2DROC[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
346	hisQ3DROC[ith]->Fill(phi,rowAbs,zIonROC,sig[i]/volume3D);
347	}
348	}
349	}
350	}
351	}
352	}
353	}
354	timer.Print();
355	delete pcstream;
356	return 0;
357	}